Supermassive black holes (BHs) are believed to reside in the centers of all massive galaxies that are not pure disks. They encode information about the formation and evolution of their host galaxies through a number of scaling relations such as the M BH -σ å , M BH -L bul , and M BH -M bul relations (Kormendy & Richstone 1995;Gebhardt et al. 2000;Ferrarese & Merritt 2000;Kormendy & Ho 2013), which relate BH mass to spheroid stellar velocity dispersion, luminosity, and mass. These scaling relations are often used to estimate BH masses in galaxies over broad ranges in both galaxy type and distance (Kormendy & Ho 2013;McConnell & Ma 2013;Saglia et al. 2016;van den Bosch et al. 2016).
Increasing the sample of measured BH masses that define these scaling relations is necessary to improve our understanding of BH-host galaxy coevolution. Questions regarding when these scaling relations came to be, the ranges of BH masses and galaxy types they apply to, and the physical mechanisms that allowed them to form remain unanswered. To complicate matters, M BH predictions from the M BH -σ å and M BH -L bul relations can strongly disagree, such as in the case of the most luminous early-type galaxies (ETGs), where the discrepancies reach an order of magnitude at M BH ∼ 10 10 M e (Bernardi et al. 2007;Lauer et al. 2007). A larger sample of reliable BH mass measurements is needed to address these questions and issues.
About 100 supermassive BH mass measurements have been obtained through modeling the motions of either stars or gas (Kormendy & Ho 2013). Accurately determining these masses requires modeling the motions of kinematic tracers that extend within the BH's radius of influence, r GM g BH 2 s » , where the BH as opposed to the host galaxy is the main contributor to their combined gravitational potential. Among current measurements, the most robust are of the Milky Way's own supermassive BH (Boehle et al. 2016;Gravity Collaboration et al. 2019), which use observations of individual stellar orbits, and of BHs in galaxies with rotating H 2 O megamaser disks whose emission originates deep within r g (Miyoshi et al. 1995;Kuo et al. 2011Kuo et al. , 2018)). For other galaxies, BH masses have been measured primarily through stellar-dynamical and ionized gas-dynamical modeling. Both methods are prone to a variety of challenges in modeling and interpreting the data. Stellar-dynamical modeling of massive ETGs is sensitive to the treatment of galaxy triaxiality, dark matter halo structure, and stellar orbital anisotropy (van den Bosch et al. 2008;Gebhardt & Thomas 2009;Shen & Gebhardt 2010), while ionized gas-dynamical modeling can be affected by noncircular motions and substantial gas turbulence (Barth et al. 2001;Shapiro et al. 2006).
Molecular gas has emerged as a dynamical tracer capable of circumventing the aforementioned issues with stellar-dynamical and ionized gas-dynamical modeling. Tracers such as H 2 , HCN, HCO + , and CO emission lines have been used to constrain BH masses in late-type galaxies (Neumayer et al. 2007;Scharwächter et al. 2013;Onishi et al. 2015;den Brok et al. 2015). In addition, a number of CO surveys have shown that a fraction of ETGs have dynamically cold and regularly rotating molecular gas disks at their centers (Combes et al. 2007;Young et al. 2011;Alatalo et al. 2013;Bolatto et al. 2017). These molecular gas disks are ideal targets for precision BH mass measurements. As with ionized gas, modeling the dynamics of molecular gas on scales comparable to the BH's sphere of influence is insensitive to factors such as the distribution of dark matter and triaxial structure, which affect stellar-dynamical models. Molecular gas also has the added benefit that it is much less turbulent than ionized gas (Davis et al. 2013a;Utomo et al. 2015;Boizelle et al. 2017). Davis et al. (2013b) demonstrated the potential of molecular gas as an effective kinematic tracer by measuring the mass of the BH in NGC 4526 with the Combined Array for Millimeterwave Astronomy. Since then, the Atacama Large Millimeter/ submillimeter Array (ALMA) has emerged as the premier radio interferometer for these types of measurements. There are now several ALMA-based BH mass measurements derived from observations of molecular gas emission on scales comparable to and even within the BHs' spheres of influence for nearby galaxies (Barth et al. 2016a;Davis et al. 2017;Onishi et al. 2017;Boizelle et al. 2019;Smith et al. 2019;North et al. 2019;Boizelle et al. 2021;Cohn et al. 2021;Nguyen et al. 2021).
In this paper, we analyze ALMA observations of NGC 1380 and NGC 6861, two galaxies that have been shown to host a central circumnuclear gas disk. Boizelle et al. (2017) mapped the distribution of CO(2-1) emission within these disks and found that they exhibited dynamically cold rotation.
NGC 1380 is classified as an SA0 galaxy in both the Third Reference Catalogue of Bright Galaxies (RC3; de Vaucouleurs et al. 1991) and in the Hyperleda database (Paturel et al. 2003).
It is located at a luminosity distance of 17.1 Mpc in the Fornax cluster based on surface brightness fluctuations from Tonry et al. (2001) after applying the Cepheid zero-point correction from Mei et al. (2007). With this assumed luminosity distance, and using an observed redshift of z = 0.00618 obtained from initial dynamical modeling results, the corresponding angular scale is 81.9 pc arcsec 1 -. We adopt the Hyperleda average stellar velocity dispersion of σ å = 215 km s -1 (Makarov et al. 2014), a total apparent K-band magnitude of m K = 6.87 mag from the Two Micron All Sky Survey (2MASS; Jarrett et al. 2003) and a bulge-to-total ratio of B/T = 0.34 from the Carnegie-Irvine Galaxy Survey (Gao et al. 2019).
NGC 6861 is located at a luminosity distance of 27.3 Mpc in the Telescopium galaxy group, which corresponds to an angular scale of 129.9 pc arcsec 1 -when using a redshift of z = 0.00944 from our initial dynamical models. This galaxy is classified as an E/S0 in Hyperleda and as an S0A-(s) in RC3. Kormendy & Ho (2013) noted that the main body of NGC 6861 does not deviate significantly from an n ; 2 Sérsicfunction profile, although the galaxy has extra central light; we adopt that paper's classification as an extra-light elliptical. In addition, we also adopt the stellar velocity dispersion of σ å = 389 km s -1 measured within the effective radius by Rusli et al. (2013), and the 2MASS total apparent K-band magnitude of m K = 7.75 mag (Vaddi et al. 2016).
For each galaxy, we obtained a BH mass by constructing gas-dynamical models and fitting directly to the ALMA data cubes. This work provides the first dynamical mass measurement of the supermassive BH in NGC 1380, and a second, independent dynamical mass measurement of the supermassive BH in NGC 6861, which was previously measured through stellar-dynamical modeling to be (2.0 ± 0.2) × 10 9 M e by Rusli et al. (2013).
Our paper is organized as follows. In Section 2, we describe both the Hubble Space Telescope (HST) and ALMA observations and the data reduction process. In Section 3, we decompose the stellar surface brightness distribution of each galaxy with Multi-Gaussian expansions and construct extinction-corrected host galaxy models. A description of our dynamical modeling formalism, including the thin disk model, is presented in Section 4. We present and compare the results of our dynamical models in Section 5. In Section 6, we compare our measurements of M BH to the local BH-host galaxy scaling relations and discuss limiting factors in the measurements. We conclude and summarize our findings in Section 7.
For NGC 1380, we retrieved and used archival HST F160W (H-band) images from HST program 11712. The observation was subdivided into four separate exposures of 299 s each that were taken with the Wide Field Camera 3 (WFC3). We processed the images with the calwf3 pipeline and subsequently combined them in AstroDrizzle to produce a cleaned and distortioncorrected image with a pixel scale of 0 08 pixel -1 . We followed a similar procedure for NGC 6861. We used archival HST data for NGC 6861 from Program 15226, which was designed to obtain host galaxy imaging to complement our ALMA program. The observation consisted of four separate exposures of 249 s each taken with the F160W filter on WFC3. We processed and combined the images with calwf3 and AstroDrizzle to produce a composite image with a 0 08 pixel -1 scale. To identify regions of substantial dust attenuation, we also obtained archival F110W (J-band) images for each galaxy. Archival F110W observations of NGC 1380 were obtained from HST program 11712 and consisted of four separate exposures of 299 s each, while F110W observations of NGC 6861 were obtained from HST program 15226 and consisted of two separate exposures of 249 s each. These were processed using the calwf3 pipeline and AstroDrizzle as above.
We obtained ALMA imaging of NGC 1380 and NGC 6861 as part of Program 2013.1.00229.S. The data were studied by Boizelle et al. (2017) as part of a larger sample of galaxies to map CO(2-1) emission in nearby ETGs, and we review their data reduction process and findings below.
NGC 1380 was observed in ALMA Band 6 for 23 minutes on both 2015 June 11 and 2015 September 18 with maximum baselines of 783 m and 2125 m, respectively. For the redshifted 12 CO(2-1) line, the observation covered a 1.875 GHz bandwidth from 228.199 to 230.074 GHz, centered at an estimated redshifted line frequency of 229.136 GHz. The frequency channel widths were 488.281 kHz, corresponding to a velocity channel resolution of 0.64 km s -1 at the redshifted frequency. For continuum emission, two separate 2 GHz spectral windows were centered at 227.210 and 244.902 GHz with 15.625 MHz channel widths, equating to velocity resolutions of 20.6 and 19.1 km s -1 , respectively. The data were initially processed through the ALMA pipeline with version 4.3.1 of the Common Astronomy Software Applications package (CASA; McMullin et al. 2007) and then imaged into data cubes using Briggs weighting with a robust parameter of 0.5 following continuum phase self-calibration and continuum subtraction in the uv plane. We reimaged the cube to have 10 km s -1 velocity channel widths (with respect to the rest frequency of the 12 CO(2-1) line) to isolate narrower line features in spatial regions close to the disk center and chose a pixel size of 0 03 to sufficiently sample the synthesized beam's minor axis. The beam's FWHM is 0 24 and 0 18 along the major and minor axis, respectively, and the beam has a position angle of 86°.9 measured east of north.
NGC 6861 was observed on 2014 September 1 in ALMA Band 6 for 24 minutes with a maximum baseline of 1091 m. The line observation was centered at an estimated redshifted 12 CO(2-1) line frequency of 228.390 GHz, while the continuum windows were centered at 226.466 and 244.098 GHz, with the same bandwidth and channel spacing properties as the NGC 1380 observation. The data were processed using CASA version 4.2.2 and imaged into a data cube with 20 km s -1 velocity channel widths following the continuum phase selfcalibration and continuum subtraction processes described for the NGC 1380 data. The NGC 6861 data cube has the following properties: a synthesized beam size of 0 32 × 0 23 with a position angle of 58°.2 and a pixel size of 0 065. Boizelle et al. (2017) determined several properties of the circumnuclear disks in NGC 1380 and NGC 6861, which we summarize here. The gas in each disk is cospatial with the dust, as seen in the HST optical images and ALMA integrated intensity maps in Figure 1. Both disks are very inclined (i ≈ 75°) and exhibit orderly rotation around their respective centers, with projected line-of-sight (LOS) velocities of ∼300 km s -1 and ∼500 km s -1 for NGC 1380 and NGC 6861, respectively. LOS velocity and dispersion maps indicate nearly circular and dynamically cold rotation about the disk centers. Independent stellar kinematic observations with the Multi Unit Spectroscopic Explorer (MUSE) have also revealed the presence of a large-scale cold disk component in NGC 1380 (Sarzi et al. 2018). The radial extents of the CO emission were measured to be 5 2 (426 pc) and 6″ (784 pc) for NGC 1380 and NGC 6861, respectively. A major axis position-velocity diagram (PVD) extracted from the NGC 1380 data cube shows a slight rise in velocity within the innermost ∼0 1, although this does not extend past the velocities observed in the outer parts of the PVD. This central upturn in gas velocity indicates the presence of a massive and compact object at the disk center. For NGC 6861, the PVD and moment maps reveal a central ∼1″ radius hole in CO emission. The gas mass of each disk was determined by summing the CO flux and assuming an α CO factor of 3.1 M e pc -2 (K km s -1 ) -1 (Sandstrom et al. 2013) as the extragalactic mass-to-luminosity ratio, a CO(2-1)/CO(1-0) ≈ 0.7 line ratio in brightness temperature units (Lavezzi et al. 1999), and a correction factor of 1.36 for helium. Given these assumptions, the gas masses were estimated to be (8.4 ± 1.6) × 10 7 M e and (25.6 ± 8.9) × 10 7 M e for NGC 1380 and NGC 6861, respectively. For NGC 1380, our assumption about CO excitation can be tested: Zabel et al. (2019) measured a CO(1-0) line flux for NGC 1380 that, in combination with the CO(2-1) line flux from Boizelle et al. (2017), implies a CO(2-1)/CO (1-0) ratio of 1.08 0.20 0.24 -+ in brightness temperature units. This is higher than our assumed value of 0.7, and implies a lower gas mass. Because a ratio >1 is unphysically high if both CO lines are tracing the same material, we consider a value ≈0.9 (still lying within the measurement uncertainties) to be more appropriate and in Section 5.2 explore the implications of the correspondingly lower gas mass for our dynamical models.
A key input to the dynamical modeling program is the stellar mass profile, M å (r), which is determined by measuring and deprojecting the host galaxy's observed surface brightness profile. One approach to obtaining M å (r) is the Multi-Gaussian expansion (MGE) method, which fits the observed brightness in galaxy images with a series expansion of 2D Gaussian functions (Emsellem et al. 1994;Cappellari 2002). For galaxies such as NGC 1380 and NGC 6861, which possess optically thick dust disks, the impact of dust attenuation on the host galaxy light is mitigated by using observations in the nearinfrared (NIR) regime, but the impact is not completely negligible. To assess the variation of H-band extinction across the dust disk, we performed a simple correction to the observed H-band major axis surface brightness profile for each galaxy as described in Section 3.1. This was done as an attempt to quantify the possible impact of dust on the host galaxy's surface brightness. We then proceeded to fit dust-masked and dust-corrected MGE models to the drizzled H-band images of NGC 1380 and NGC 6861, as described in Sections 3.2 and 3.3.
The central dust disk in each galaxy is clearly visible in Figures 1 and 2, due to its dimming and reddening of the observed stellar light. We attempted to estimate the amount of dust extinction with a color-based correction method, which was complicated by the fact that the dust disk is embedded within the galaxy and cannot be treated as a simple foreground screen. To estimate the amount of dust extinction, we extracted and corrected the observed H-band major axis surface brightness profile of each galaxy. Using the sectors_photometry routine from the MgeFit package in Python (Cappellari 2002), we plotted the surface brightness profile of each galaxy in Figure 3. In NGC 1380, a slight dip in the profile can be seen at around r = 4″, which marks the outer edge of the dust disk, while a more noticeable dip is seen between ∼1″ and 5 5 in NGC 6861.
To determine the pixels that are most affected by dust, we created J -H color maps as seen in Figure 2. We note that all H-and J-band magnitudes in this work are in the Vega magnitude system. The color maps revealed that the dust extends from the nuclei, and the pixels most affected by dust are about 0.25 mag redder than the median J -H color of ∼0.80 mag outside the disk. Furthermore, each nucleus has a bluer color that is nearly identical to the median color outside the disk. A blue nucleus suggests the presence of star formation or a weak active galactic nucleus (AGN); we discuss these possibilities below.
We attempted to correct the major axis H-band surface brightness profiles by examining Δ(J -H), the observed color excess relative to the median J -H color outside the disk, along the major axis. Using Equation (1) from Boizelle et al. (2019), which predicts the ratio of observed to intrinsic integrated stellar light based on the embedded-screen model described by Viaene et al. (2017), we generated a model Δ(J -H) curve as a function of intrinsic V-band extinction, A V , to compare with the observations.
The embedded-screen model assumes that the obscuring dust lies in a thin, inclined disk that bisects the galaxy, and that the fraction of stellar light originating behind the disk, b, is obscured by simple screen extinction, while the fraction of stellar light in front of the disk, f, is unaffected. In addition, the model assumes that there is no scattering of stellar light back into the LOS, and that the J -H color outside of the disk is the intrinsic color of the host galaxy. One important aspect of this color excess model is that Δ(J -H) is not a strictly monotonically increasing function of intrinsic extinction. As seen in Figure 4, the color excess increases approximately linearly with increasing A V up to a turnover point, after which it will begin to decrease to zero as the light originating behind the disk becomes completely obscured. As a result, there are two possible A V values for a given Δ(J -H). To maintain a one-toone correspondence between Δ(J -H) and A V , we considered only the low-extinction branch of the curve and therefore adopted the lesser of the two possible A V values. We inverted the relationship to derive A V as a function of observed Δ (J -H) by fitting a third-order polynomial up to the turnover point and determining its inverse. Using a standard Galactic (R V = 3.1) extinction curve (Rieke & Lebofsky 1985) where A H /A V = 0.175, and assuming that the fractions of stellar light originating in front of and behind the disk are equal ( f = b = 0.5) along the major axis, we associated our observed major axis Δ(J -H) values with corresponding values of A H .
We generated point-by-point corrections to our major axis surface brightness profiles using the fact that our modeled A H values only applied to the fraction of light originating behind the disk. The corrected values are shown in blue in Figure 3 for points within the dust disk for each galaxy. For NGC 1380, the corrected surface brightness profile still exhibits a slight dip near the edge of the dust disk. For NGC 6861, where we can anchor the surface brightness profile inside and outside the dust disk, it is clear there is still some remaining extinction, as the decrease in observed surface brightness is still visible in the corrected profile.
Since the method described above appears to undercorrect for extinction, we also tried applying the method using the high-extinction branch of the Δ(J -H) versus A V curve to both galaxies; however, this approach led to overcorrection all along the major axis of each galaxy, as each point's H-band surface brightness was raised by nearly the theoretical maximum of 0.75 mag arcsec -2 . While our method provides some insight on how extinction varies across the disk, based on the results for these two galaxies using both branches of the color excess curve, it is clear that a simple extinction correction for a thin embedded disk does not fully correct for dust extinction or give us accurate host galaxy profiles. Thus, we opted to create dustmasked and dust-corrected MGEs to model the host galaxy's light following methods used for similar galaxies by Boizelle et al. (2019Boizelle et al. ( , 2021) ) and Cohn et al. (2021).
Before we constructed an MGE model for NGC 1380, we created a mask that isolated the host galaxy light from the light of foreground stars and background galaxies and identified pixels affected by dust. For the NGC 1380 drizzled H-band image, we masked out foreground stars and background galaxies in the image and corrected for a foreground Galactic reddening in the H band of A H = 0.009 mag based on reddening measurements from Sloan Digital Sky Survey data by Schlafly & Finkbeiner (2011). Using the J -H color map we constructed earlier, we also masked pixels that had J -H > 1.05 mag, which were on the disk's near side. This step prevented pixels with the most apparent dust obscuration from being used in the MGE fit, but there is still clear evidence of extinction in other regions of the disk.
We modeled the observed H-band surface brightness within the inner 10″ × 10″ with an MGE created in GALFIT (Peng et al. 2002) and required each Gaussian component to have the HST major axis H-band surface brightness profiles to their respective dust-corrected models. The red points are the observed values from the H-band images, while the blue points are the dust-corrected values described in Section 3.1. The different lines in each panel correspond to extracted major axis surface brightness profiles for our 2D MGE models, which are described in Sections 3.2 and 3.3. For NGC 1380, we mark the outer edge of the dust disk with a vertical dashed line and indicate that the dust extends down to the nucleus with an arrow, while for NGC 6861, we indicate the inner and outer boundaries of the dusty region with the two dashed lines. (2019) and a standard Galactic (R V = 3.1) extinction curve (Rieke & Lebofsky 1985), under the assumptions that along the major axis, the fractions of light in front of and behind the disk are equal ( f = b = 0.50), and that only the starlight originating from behind the disk is subject to dust extinction. The gray region left of the dotted line indicates the part of the curve we used to define a one-to-one correspondence between Δ(J -H) and A V and A H for the low-extinction branch of the curve. same center and position angle. To account for the HST PSF, we generated a model H-band PSF using Tiny Tim (Krist & Hook 2004). This PSF was drizzled and dithered in the same manner as the H-band image and, along with our mask, was used during the GALFIT optimization to create the MGE. We refer to this initial, dust-masked MGE as our A H = 0.00 mag model, since it does not attempt to correct for the impact of extinction at locations that were not masked out.
A robust pixel-by-pixel dust correction model would require radiative transfer modeling to account for factors such as disk geometry, thickness, scattering from dust, and extinction within the disk itself (De Geyter et al. 2013;Camps & Baes 2015). Additionally, light originating from recent star formation or a weakly active nucleus would add further complications to a dust correction model. In our J -H color map, the nucleus of NGC 1380 is bluer than the most reddened pixels in our mask by about ∼0.2 mag, suggesting the presence of star formation and/or a weak AGN. Zabel et al. (2020) used combined MUSE and ALMA data to study the relationship between molecular gas surface density and star formation rate in NGC 1380. They concluded that there was no Hα emission from star formation, and that the presence of Hα in NGC 1380 was primarily due to what they defined as composite regions such as shocks or an AGN. Indeed, through integral field spectroscopy, Ricci et al. (2014) determined that NGC 1380 contains a low ionization nuclear emission-line region. Viaene et al. (2019) also studied the dust mix and gas properties in NGC 1380 with MUSE observations and detected low-level star formation within the inner portion of the disk. They construct 2D A V maps of the dust lane area by comparing MGE model fits (after having masked out the dust lane) to MUSE V-band images and estimated a maximum A V value of 1.00 mag, corresponding to A H ≈ 0.18 mag for a standard Galactic (R V = 3.1) extinction curve. In addition, they use 3D radiative transfer models to reproduce the observed V-band attenuation (defined as the combination of extinction and scattering of light back into the LOS) curve. However, their methods assume that only the near side of the dust disk experiences any V-band extinction, whereas our J -H map shows that pixels on the far side are redder relative to the median color outside the disk, indicating that light from the far side is still affected by extinction.
We used the simpler method described by Boizelle et al. (2019), which assumed an analytic surface brightness profile model to correct for dust extinction. Their method examined the impact of extinction on the inferred host galaxy circular velocity profile by adjusting the central H-band surface brightness profile to correct for three fiducial values of dust extinction. We chose the same values of H-band extinction as Boizelle et al. (2019), which were A H = 0.31, 0.75, and 1.50 mag. These values correspond to fractions of one-fourth, onehalf, and three-fourths of the stellar light originating behind the dust disk. We emphasize that these values were chosen to explore the impact of dust on the inferred host galaxy models (and subsequently, the values of M BH derived from our dynamical models) over a range in extinction and note that among the major axis surface brightness profiles shown in Figure 3, the A H = 0.31 mag MGE model most closely resembles the observed profile after correction using the lowextinction branch of the reddening curve.
We created three dust-corrected MGE models based on the three fiducial H-band extinction values mentioned above. We followed the steps outlined by Boizelle et al. (2021) and describe our process below. To start, we fit a 2D Nuker model (Faber et al. 1997) in GALFIT to the central 10″ × 10″ region of the H-band image, using the same dust mask and PSF model as before. Nuker profiles are known to effectively model the central surface brightness profiles of ETGs, and we can easily adjust their parameters to produce dust-corrected models matching the H-band image. The Nuker model's surface brightness profile is characterized by inner and outer power-law profiles, with γ and β representing the inner and outer profile logarithmic slopes. The transition between these two regimes occurs at a break radius, r b , and the transition sharpness is controlled by the parameter α. We allowed all free parameters of the Nuker model to vary in this initial fit. The Nuker model parameters converged to α = 0.42, β = 1.49, γ = 0.31, and r b = 2 5 (≈ 200 pc). These parameters characterized the Nuker model fit to the H-band image prior to any dust correction. We then manually corrected the central surface brightness values of the H-band image for extinction levels of A H = 0.31, 0.75, and 1.50 mag, and fit three separate Nuker models to these three dust-corrected images, keeping all parameters other than γ fixed. This approach allowed the Nuker model to adjust its inner slope to the dust-corrected values of the central pixels, but retain its outer slope shape from the initial fit. For extinctions of A H = 0.31, 0.75, and 1.50 mag, the value of γ converged to values of 0.39, 0.44, and 0.47, respectively. Finally, to create dust-corrected MGE models, we replaced the H-band data within the disk region with the corresponding pixels in the Nuker models, and fit MGE models in GALFIT to these dust-corrected H-band images without using a mask. The major axis surface brightness profiles of these three MGE models are shown in Figure 3, and a plot of their enclosed mass profiles is shown in Figure 5. As we show in Section 5.2, our best-fit dynamical model uses the A H = 0.31 mag MGE model. We display and compare this model's isophotal contours to those of the data in Figure 6. While there is good agreement between data and model outside of the central dust lane, there is some deviation between the two within the dusty region. The observed isophotes become nonelliptical toward the center, which we attribute to the presence of the dust disk. The A H = 0.31 mag MGE model's components are listed in Table 1, while the other MGE models' components are listed in Table 3 in the Appendix.
We created our MGE models for NGC 6861 in a different fashion from our models for NGC 1380 given the differences in their dust disks. We started again by creating a J -H color map using the drizzled J-and H-band HST images of NGC 6861 to identify the pixels most affected by dust and corrected for Galactic reddening based on a foreground A H = 0.028 (Schlafly & Finkbeiner 2011). The color map revealed the presence of a ring-like structure within the disk with a 1″ radius hole at its center, as seen in Figure 2. A measurement of the surface brightness along the major axis of the disk (shown in Figure 3) revealed a clear decrease of stellar light due to dust. This decrease is most noticeable between 1″ (the outer radius of the central hole) and 5 5 (the outer edge of the dust disk). The central hole is also visible in the absence of CO emission within the inner 1″ of the PVD shown in Figure 7.
Given the lack of dynamical tracers in the central region, we wanted to test how our choice of host galaxy model impacted the inferred value of M BH . We explored systematic effects due to the choice of surface brightness model used to interpolate over the dusty region, by constructing MGEs using two different models in the correction of the H-band image for extinction: (1) a 2D Nuker model, and (2) a 2D Core-Sérsic model (Graham et al. 2003). Similar to the Nuker model, the Core-Sérsic model was designed to characterize the surface brightness profiles of ETGs. Unlike the Nuker model, however, it characterizes the outer structure of ETGs with a Sérsic profile (Sérsic 1963) and the inner structure as a power law.
To fit a Nuker model to the H-band image, we followed a process similar to that for NGC 1380, but we did not make any adjustments to the central pixels in the drizzled H-band image. We first created H-band PSF models in Tiny Tim that were dithered and drizzled in an identical fashion to the H-band image, and built a mask for foreground stars, background galaxies, and the dust disk itself. Because the J -H color map indicated a lack of color excess in the central hole, we masked the entire dust disk as seen in the J -H image, but kept the pixels within the hole to anchor the model. We fit the inner 10″ × 10″ of the H-band image with a Nuker model in GALFIT and allowed all free parameters to vary. The Nuker model parameters converged to α = 0.65, β = 1.29, γ = 0.0002, and r b = 0 31(≈ 40 pc). Finally, we replaced pixels in the original H-band image located in the dust disk with the corresponding pixels in the Nuker model, and proceeded to fit this image with an MGE model in GALFIT without using a mask. We measured and compared this MGE model's major axis surface brightness profile and isophotal contours to those of the H-band image in Figure 3 and Figure 6; we refer to this model as our Nuker interpolation model.
We constructed a 2D Core-Sérsic model for NGC 6861ʼs Hband image using the imfit program (Erwin 2015) and used the same mask, model PSF, and fitting region as for our Nuker model. The parameters that characterize the Core-Sérsic model include the Sérsic index, n, break radius, r b , effective half-light radius, r e , inner slope parameter, γ, and transition sharpness, α. The optimization in imfit converged on n = 7.1, r b = 3 5(≈455 pc), r e = 7 9(≈1 kpc), γ = 0.61, and α = 2.13. We replaced the pixels in the dust disk region in the original Hband image with the corresponding pixels in the Core-Sérsic model, and proceeded to fit this new image with an MGE in GALFIT. We deprojected this MGE in an identical fashion to the MGE created with the 2D Nuker model, and refer to it as our Core-Sérsic interpolation model. We extracted its major axis surface brightness profile and compared it with both the Hband image and the Nuker interpolation in Figure 3. Over the extent of the dust disk, the Core-Sérsic interpolation produces higher corrected surface brightness values than the Nuker interpolation. At the nucleus, the Nuker interpolation matches the observed central surface brightness better than the Core-Sérsic interpolation, whose innermost point slightly exceeds the observed value. A comparison with the observed H-band isophotes is shown in Figure 6. As in the case of NGC 1380, the observed isophotes became noticeably nonelliptical toward the center, although there appears to be reasonable agreement between the data and models within the central hole region. The MGE components of the Nuker interpolation are listed in Table 1, while the components of the Core-Sérsic interpolation are listed in Table 4 in the Appendix.
Our dynamical modeling formalism is a Python-based adaptation of the ALMA gas-dynamical modeling framework described by Barth et al. (2016a) and Boizelle et al. (2019), which was written in IDL and was used by those authors to measure the BH masses in NGC 1332 and NGC 3258. We describe the methods used in the Python version and the modifications that differentiate it from its IDL antecedent.
Modeling the observed gas kinematics in an ALMA data cube relies on a few key steps and assumptions. First, we assume that the gas is distributed in a thin disk and is in circular rotation. A model velocity field is built on a grid that is oversampled relative to an ALMA spatial pixel by a factor of s = 3, such that each pixel is subdivided into s × s = 9 subpixels in order to model steep velocity gradients near the disk's center. The disk's velocity field is determined by the enclosed mass at a given radius, which consists of a central BH, the stellar mass profile of the host galaxy and a corresponding mass-to-light (M/L) ratio ϒ, and the mass profile of the gas disk. For a given disk inclination i and a major axis position angle Γ (both of which are free parameters), and an assumed (fixed) distance to the galaxy, D, we calculate the LOS projection of this velocity field as seen on the plane of the sky. The construction and geometry of the model disk are as described by Macchetto et al. (1997) and Barth et al. (2001). The LOS velocity projections are used to generate a model cube that we can compare directly to the ALMA data. For each subpixel with CO emission, we assume that the emergent line profile along the spectral dimension is intrinsically Gaussian. , where the v å,
H ¡ values in Table 2. The resolution of the ALMA observation is denoted by the vertical dotted line and is comparable to the BH's expected radius of influence. The BH mass determined from our fiducial model is represented by the horizontal dashed line, and the range of BH masses determined from models A-D is indicated by the gray shaded region. The Gaussian's line centroid and line width can be calculated at each subpixel by transforming both the LOS velocity projections and a spatially uniform turbulent velocity dispersion term, σ(r) = σ 0 , into observed frequency units using the redshift z obs (related to the systemic velocity through v sys = cz obs ). The model cube must have its line profiles scaled by a model CO flux map, have each of its frequency slices convolved with the ALMA synthesized beam, and be downsampled to an appropriate resolution before being fitted to the ALMA data cube. We discuss these steps in further detail in the subsequent paragraphs and in Section 4.2.
In total, our dynamical models use a minimum of nine free parameters: the BH mass M BH , the stellar H-band M/L ratio ϒ H , the disk's dynamical center in pixels (x c , y c ), the disk's inclination and major axis position angle i and Γ, the turbulent velocity dispersion σ, the observed redshift z obs , and a fluxscaling factor F 0 , that correctly normalizes the model to the data.
The circular velocity v c (relative to the disk's systemic velocity, v sys ) as a function of radius is calculated as
where v å,MGE and v gas are the circular velocities due to the gravitational potential of the stars and the gas disk, respectively. The BH is modeled as a point mass, while the stellar and gas mass distributions are radially extended and are constructed using different methods. We modeled the stellar mass distribution using the MGE method described in Section 3.2 and Section 3.3. We deprojected the MGE under the assumptions that NGC 1380 and NGC 6861 are oblate and axisymmetric and have inclination angles of 77 • and 73 • , respectively, based on initial gas-dynamical modeling runs. We calculated the contribution to the circular velocity from the stars in the midplane of each disk by using the mge_vcirc routine from the JamPy package in Python (Cappellari 2008) to derive a fiducial velocity profile from our MGEs. Ideally, one should match the stellar inclination angle of the MGEs to the inclination angle found for the gas disk, as mismatches between the two lead to nonequilibrium configurations for the disk. However, this matching process is difficult to implement within our framework, and we found that the differences between the stellar and gas inclination angles were small (<2°) in both NGC 1380 and NGC 6861. We will explore this aspect of the modeling process in a future work. We use ϒ MGE = 1 Table 1 H-band MGE Parameters k
NGC 1380 NGC 6861 A H = 0.31 mag A H = 0.00 mag (Nuker Model) 1 5.560 0.059 0.739 4.960 0.091 0.934 2 4.996 0.197 0.710 4.874 0.179 0.563 3 4.662 0.391 0.726 4.455 0.369 0.624 4 4.513 0.743 0.722 4.175 0.370 0.624 5 4.356 1.461 0.807 4.334 0.372 0.625 6 3.973 3.414 0.613 4.228 0.630 0.788 7 3.386 3.880 0.997 4.061 0.787 0.686 8 3.757 6.018 0.723 4.391 1.120 0.512 9 3.382 12.932 0.732 4.165 1.824 0.879 10 3.043 18.800 0.400 4.085 4.227 0.542 11 2.689 42.731 0.400 3.658 7.591 0.506 12 2.101 55.077 0.851 3.299 13.064 0.562 13 0.926 92.556 0.948 2.508 27.828 1.000 14 L L L 1.662 33.017 0.641 15 L L L 1.684 52.142 0.883 16 L L L 1.700 62.263 0.743 Note. NGC 1380 and NGC 6861 MGE solutions created from the combination of HST H-band images and best-fitting GALFIT Nuker models. These correspond to the A H = 0.31 mag MGE model for NGC 1380 and the Nuker interpolation MGE model for NGC 6861. As described in Section 5, these two MGEs are used in the dynamical models with the lowest χ 2 . The first column is the component number, the second is the central surface brightness corrected for Galactic extinction and assuming an absolute solar magnitude of M e,H = 3.37 mag (Willmer 2018), the third is the Gaussian standard deviation along the major axis, and the fourth is the axial ratio, which was constrained to have a minimum value of 0.400 to allow for a broader range in the inclination angle during the deprojection process. Primes indicate projected quantities. when deriving v å,MGE for each galaxy. At each model iteration, v ,MGE 2
is scaled by the ratio ϒ H /ϒ MGE , which scales the stellar mass profile by the free parameter ϒ H .
As stated earlier, Boizelle et al. (2017) created mass profiles for both gas disks by averaging the CO flux in elliptical annuli centered on the continuum peaks. They determined M gas to be (8.4 ± 1.6) × 10 7 M e and (25.6 ± 8.9) × 10 7 M e for NGC 1380 and NGC 6861, respectively, but their gas mass profiles did not assume specific shapes for the mass distributions. We assumed the mass was distributed in a thin disk and numerically integrated the projected surface mass densities to determine each gas disk's contribution to the circular velocity (v gas ) using Equation (2).157 from Binney & Tremaine (2008).
We disregard the mass contribution of dark matter in our models, as the stars are expected to dominate the mass budget across the length of the circumnuclear disk. For NGC 6861, we estimated the dark matter mass within the region we fit our models (r ≈ 400 pc) by integrating the spherical cored logarithmic density profile used by Rusli et al. (2013) in their stellar-dynamical BH mass measurement. Their model suggests an enclosed dark matter mass between 10 6 -10 7 M e , which is lower than our estimated stellar mass by two to three orders of magnitude.
The Gaussian line profiles at each point on the disk must be weighted by an observed CO flux map obtained from the ALMA observation. To create this flux map, we first visually identified channels in the data cube that contained CO emission. In these channels, we created a unique mask that separated pixels with visible emission from those without any. Spatially, each mask has the same size as a single frequency slice, with pixel values set to unity if the corresponding data cube pixel displays CO emission and zero if not. A channel that displays no emission would have a corresponding mask with all of its elements set to zero. An entire mask is 3D, with the same dimensions as the ALMA data cube. We then multiplied each slice of this mask by the corresponding slice in the ALMA data cube, and summed the products along the spectral axis. This approach produced a less noisy image of the CO flux than if we had simply summed the data cube across channels with visible emission without any masking. To deconvolve this image, we applied five iterations of the Richardson-Lucy algorithm (Richardson 1972;Lucy 1974). The deconvolution is performed with an elliptical Gaussian PSF that matches the specifications of the ALMA synthesized beam and uses the Richardson-Lucy algorithm implemented in the scikitimage package in Python (Van der Walt et al. 2014). The deconvolved image is initially constructed on the original ALMA pixel scale, with each pixel then subdivided into an s × s grid of subpixels that matches the dimensions of the oversampled model grid. The scaled and deconvolved CO flux map is normalized so that the line profiles at each subpixel element for a given original ALMA pixel have equal fluxes. Thus, if the total flux in an original ALMA pixel is F, each subpixel in the s × s grid has a flux of F/s 2 .
The next steps consist of rescaling the oversampled model back to the original ALMA pixel scale, convolving each frequency channel within it with the ALMA synthesized beam, and minimizing χ 2 between data and model. Ideally, beam convolution should occur on the oversampled spatial grid for the highest model fidelity. However, beam convolution is the most time-consuming part of the entire modeling process and becomes prohibitively slow for oversampling factors of s > 3. Both Barth et al. (2016b) and Boizelle et al. (2019) found that modeling results do not change appreciably if the convolution step is done on the original ALMA pixel scale, so we followed the same approach. We summed each s × s group of subpixels in our oversampled model to form a single pixel on the original ALMA scale, and then convolved each frequency slice of our model with the ALMA synthesized beam, using the convolution implementation in the astropy package for Python (Astropy Collaboration et al. 2013, 2018).
For a given model parameter set, we create a simulated data cube with the same spatial and spectral dimensions as the ALMA data. Therefore, our models can be fitted directly to the ALMA data cubes and can be optimized by χ 2 minimization. We optimized models with the Levenberg-Marquardt algorithm (Levenberg 1944;Marquardt 1963) within the LMFIT framework (Newville et al. 2016) in Python and fitting to pixels that lay within the elliptical regions illustrated in the data moment 0 maps in Figures 8 and 9, and within the frequency channels that span the full width of the CO emission line for each pixel. We describe the fitting regions in detail in Section 4.2.2.
In order to calculate χ 2 and assess the goodness-of-fit of our models to the ALMA data, we require an estimate of the flux uncertainty at each data pixel. The most straightforward approach here would be to calculate the standard deviation of pixel values in emission-line-free regions of the data cube. However, the background noise in ALMA data cubes is correlated on scales comparable to the synthesized beam in each frequency channel. This correlation prevents the determination of a meaningful χ 2 value without appropriate adjustments. Ideally, one would calculate a covariance matrix accounting for these correlations to compute χ 2 , but such an approach would be computationally expensive and challenging to implement. Barth et al. (2016b) and Boizelle et al. (2019) adopted the simpler approach of rebinning the data by blockaveraging over m × m pixel blocks within each frequency channel, where the value of m was the approximate number of pixels across the width of the synthesized beam. Their method creates a data cube with a scale of approximately one rebinned pixel per synthesized beam, and mitigates the noise correlation among neighboring pixels. They then measured the standard deviation of emission-free pixels in the rebinned data cube to produce a unique value of flux uncertainty for each frequency channel, and similarly rebinned their models to compute χ 2 on the block-averaged scale.
We also incorporated the effects of the ALMA primary beam on the background noise level. Prior to primary beam correction, the noise level in an ALMA data cube is spatially uniform, but post-correction it increases with distance from the phase center. Dynamical models are created and fitted to data cubes that have been corrected for the primary beam attenuation, so we incorporated this spatial modification into our noise model. As part of the ALMA data reduction process, a primary beam cube is generated along with the beamcorrected data cube. Multiplying the corresponding slices of these cubes together generates an uncorrected version of the data in which the background noise is spatially uniform. At this step, we block-average the data to roughly the size of the synthesized beam, using 7 × 7 pixel blocks for the NGC 1380 data cube and 4 × 4 pixel blocks for the NGC 6861 data cube. Once the data have been rebinned, we measure the standard deviations of pixel values in blank regions of each frequency channel, and populate an array having the spatial dimensions of a block-averaged image with the value of the standard deviation at each element. To replicate the spatial modification of the noise in each channel, we block-averaged the primary beam cubes over the same pixel blocks as was done for the data and divided the block-averaged array of standard deviations by the block-averaged primary beam cube at the same frequency. In essence, we create a block-averaged noise cube that captures both the spatial and frequency dependence of the noise, which we use to compute χ 2 . This approach differs from previous methods where the given background noise is assumed to be spatially uniform across a given frequency slice (Barth et al. 2016b;Boizelle et al. 2019).
Although our noise model is designed to represent the rms noise in emission-line-free regions of each frequency channel of the data cube, an additional complication is that the mean background level can be slightly offset from zero (e.g., as a residual of imperfect passband calibration or continuum subtraction). If this is the case, the line-free regions of a cube will indirectly contribute to elevated χ 2 values for model fits.
For both galaxies, we find roughly equal number of channels having positive and negative mean background levels, with typical magnitudes ∼10% of the respective per-channel rms noise levels. As a simple test, we empirically measured the mean background level in each frequency channel included in our fit and added this value into the corresponding channels in our synthetic model cubes. We found that our values of reduced χ 2 were smaller by about ∼1% with this adjustment, and hence the impact of these background levels can be regarded as minimal.
We computed χ 2 over elliptical fitting regions to assess the goodness-of-fit of our dynamical models. These elliptical fitting regions were centered on the disk centers and used the axial ratios and major axis position angles of the gas disks, as measured by Boizelle et al. (2017). The size of the fitting region can influence the inferred value of M BH . While fitting models to the entire disk uses all of the available data, the majority of pixels in the fit are at radii much greater than the BH's radius of influence. In this regime, the uncertainty in the stellar mass profile accounts for a large portion of the error budget, and model fits can lead to tight statistical constraints on both ϒ H and M BH . If the assumed and intrinsic shapes of the stellar mass profile are discrepant, full-disk fits can force M BH to an inaccurate, but highly precise value. Alternatively, fitting to smaller regions can mitigate effects from discrepancies in the stellar mass profile because the BH mass represents a larger fraction of the total enclosed mass. Smaller fit regions can also limit systematic effects due to the structural mismatch of a thin disk model with a mildly warped disk. We discuss our selected fitting regions for NGC 1380 and NGC 6861 below.
For NGC 1380, we initially created an ellipse centered on the disk center with an axial ratio of q = 0.27, and a position angle of Γ = 187°based on results from Boizelle et al. (2017). For our fiducial dynamical model, we chose to fit within an ellipse that encompassed the inner half of the CO disk in order to limit the sensitivity of our dynamical models to the shape of the stellar mass profile and the disk's slightly warped structure. We also modified the size of this ellipse to see how the choice of fit region affected the inferred value of M BH in Section 5.2. Our final fitting ellipse has a semimajor axis of a = 2 05 and a semiminor axis of b = 0 55. This ellipse was used across 62 consecutive frequency channels that spanned the full width of the visible CO emission in the data and can be seen in Figure 8. On the final rebinned scale, this choice of spatial and spectral regions resulted in 61 block-averaged pixels over 62 frequency channels for a total of 3782 data points used to calculate χ 2 .
For NGC 6861, we initially followed the same procedure, starting with the values of q = 0.32 and Γ = 141°found by Boizelle et al. (2017). However, the disk structure in NGC 6861 is more complicated than in NGC 1380. The NGC 6861 gas disk contains a central hole that is ∼1″ in radius along the major axis. Thus, the innermost CO emission is at a radius that is three times larger than the BH's estimated radius of influence. Additionally, the presence of rings and spiral-like substructure can be seen toward the edge of the disk. Fitting models to the entire disk led to reduced χ 2 values between 2.5 and 3, as the thin disk models struggled to reproduce kinematic features in the outer disk. The inner half of the gas disk shows the most regularity in its structure, and we found that fitting dynamical models in this region led to lower reduced χ 2 values and better overall fits to the data. Therefore, we created an elliptical fitting region with dimensions a = 3″ and b = 0 96, as seen in Figure 9. In order to prevent pixels within the hole from contributing to the fit, we masked out a 1″ ellipse with the same axial ratio (q = 0.32) at the center of our fitting region, which yielded our final annular fitting region. Along the spectral axis, we fit across 52 frequency channels that extended slightly beyond the channels with visible emission. On the final rebinned scale, with 75 rebinned pixels per channel, we included a total of 3900 data points in the fit.
We present results for four models for NGC 1380, which we refer to as models A, B, C, and D in Table 2. The key difference among them is the input host galaxy circular velocity profile, based on one of the four MGE models described in Section 3.2, which accounted for four fiducial values of central dust extinction (A H = 0.00, 0.31, 0.75, and 1.50 mag for models A, B, C, and D, respectively). These dynamical models use a uniform turbulent velocity dispersion across the entire disk and are optimized over the elliptical region described in Section 4.2.2.
Models A-D yield best-fit values of M BH in the range of (1.02-1.85) × 10 8 M e , ϒ H between 1.30 and 1.42, and a range in reduced χ 2 2 ( ) c n between 1.525 and 1.563 over 3773 degrees of freedom. Our measured range of ϒ H is slightly above the predictions made from single stellar population (SSP) models (Vazdekis et al. 2010) that assume either a Kroupa (2001) or Chabrier (2003) initial-mass function (IMF) and is lower than predictions made with a Salpeter (1955) IMF. These SSP models assume an old stellar population (10-14 Gyr) and solar metallicity, which are consistent with 2D IMF analyses of NGC 1380 by Martín-Navarro et al. (2019). As seen in Table 2, all other free parameters remained virtually unchanged among the different dynamical models. We discuss the interplay between M BH and ϒ H in Section 6.1.
The major axis PVD, moment maps, and example line profiles for our best-fit model (model B) for NGC 1380 are presented and compared with the ALMA data in Figures 7, 8, and 10. The moment maps and PVD reveal that the data and model are in good agreement over a majority of the disk, although a mismatch in the observed surface brightness is seen in both the moment 0 (surface brightness) map and in the structure of the PVD, especially within the innermost ∼0 5. Moment 1 (v LOS ) maps show that data and model velocity fields are also in good agreement, although differences of ∼30 km s -1 are noticeable toward the disk edge, outside the fit region, and at the disk center, where the impact of beamsmearing is most severe. Most likely, these differences are a result of the differences between the model and intrinsic stellar velocity profiles. The extracted line profiles of the blockaveraged data and best-fit model highlight our models' ability to reproduce the observed shapes of the line profiles, even when they display asymmetric structure.
To determine the statistical uncertainties of the free parameters, we performed a Monte Carlo simulation for model B. We created 150 realizations of the best-fit model by adding noise to each pixel of the model cube. The value of the noise at each pixel was determined by choosing a random value drawn from a Gaussian distribution with a standard deviation equal to the value of the corresponding pixel in the noise cube described in Section 4.2.1. We re-fit to each noise-added model realization using the values found in Table 2 as initial guesses; the standard deviation of each recovered parameter was identified as the 1σ uncertainty. For M BH , we found a tight distribution centered at our initial guess of M BH = 1.47 × 10 8 M e with a standard deviation of 2 × 10 6 M e , or 1.4% of the mean. The statistical uncertainties for the free parameters are listed under model B's best-fit values in Table 2. Based on other Monte Carlo simulations we ran, these statistical uncertainties are representative of those for models A, C, and D.
While the statistical uncertainties from our Monte Carlo simulation are small, there are several other sources of uncertainty that stem from the choices we made when building our dynamical models. Thus, we conducted numerous tests to determine the impact these choices had on the value of M BH .
Dust extinction: The value of M BH in our model optimizations is highly sensitive to the choice of MGE model. Using the initial, dust-masked MGE, our models converge on M BH = 1.85 × 10 8 M e . As we increase the central extinction from A H = 0.00 mag to A H = 0.31 mag, corresponding to a loss of 25% of the total stellar light behind the dusty disk, M BH decreased to 1.47 × 10 8 M e , representing a ∼20% decrease from the initial fit. This particular model also shows a decrease in the resultant 2 c n , as model A with the initial MGE returns Table 2 Dynamical Modeling Results Model MGE M BH ϒ H i Γ σ 0 x c y c v sys F 0 2 c n (M e /L e ) ( °) ( °) ( km s -1 ) (arcsec) (arcsec) (km s -1 ) NGC 1380 (A H mag) ( 10 8 M e ) A 0.00 1.85 1.33 76.9 187.1 10.8 0.016 0.013 1853.83 0.99 1.544 B 0.31 1.47 1.42 76.9 187.2 10.5 0.017 0.013 1853.86 0.99 1.525 (0.02) ( 0.003) ( 0.04) (0.04) ( 0.13) ( 0.001) (0.001) ( 0.15) ( 0.003) C 0.75 1.27 1.36 76.9 187.2 10.5 0.017 0.013 1853.86 0.99 1.545 D 1.50 1.02 1.30 76.8 187.2 10.5 0.017 0.013 1853.88 0.99 1.563 NGC 6861 (Nuclear Profile) (10 9 M e ) E Nuker 1.13 2.52 72.7 142.5 7.2 0.057 0.038 2795.63 1.03 1.987 (0.04) ( 0.01) ( 0.07) (0.05) ( 0.29) ( 0.002) (0.002) ( 0.29) ( 0.006) F Core-Sérsic 2.89 2.14 73.6 142.6 7.4 0.049 0.046 2795.65 1.04 2.004 Note. Best-fit parameter values obtained by fitting thin disk dynamical models to the NGC 1380 and NGC 6861 CO(2-1) data cubes. We derive 1σ statistical uncertainties for the parameters of fiducial models B and E, based on a Monte Carlo resampling procedure described in Section 5.1 and list them under the results for models B and E. These models have 3773 and 3891 degrees of freedom, respectively. The major axis PA, Γ, is measured east of north for the receding side of the disk. The disk dynamical center, (x c , y c ) is measured in arcsecond offsets from the nuclear continuum centroids for NGC 1380 and NGC 6861 determined in Boizelle et al. (2017). The observed redshift, z obs , is used in our dynamical models as a proxy for the systemic velocity of the disk, v sys , in the barycentric frame via the relation: v sys = cz obs and is used to translate the model velocities to observed frequency units.
respectively. However, both models C and D result in higher 2 c n values. We chose the A H = 0.31 mag MGE model (model B) to use as our fiducial model for a number of reasons. First, this MGE model accounts for the impact of dust extinction, whereas the initial MGE simply masks dust out. Next, the dynamical model that uses the A H = 0.31 mag MGE has the lowest value of χ 2 out of all of our models, signifying the best overall match to the ALMA data. Lastly, extracting the major axis surface brightness profile from this MGE model reveals that it most closely resembles our surface brightness profile after correction for the low-extinction branch of the reddening curve shown in Figure 3. Therefore, for all of the remaining systematic tests, we use the A H = 0.31 mag MGE as our host galaxy model. The model line profiles are assumed to be intrinsically Gaussian at the subpixel level, with the asymmetric structure coming from beamsmearing. The noise band represents values in the range of data ± 1σ. These line profiles were extracted from the cubes on the final scale of 1 block-averaged pixel per synthesized beam. The x and y labels indicate the pixel locations in terms of offsets in arcseconds from the disk dynamical centers, which are given in Table 2.
Radial motion: Although there is no evidence of strong deviations from circular motion in the NGC 1380 gas disk, we constructed a simple model that allows for radial motion in our dynamical models. We followed an approach similar to Boizelle et al. (2019) and Cohn et al. (2021) and added a radially inward velocity term to our dynamical models. We included an additional free parameter, α, which lies in the range [0, 1] and controls the balance between pure rotational (α = 1) and radially inflowing (α = 0) motion. Mathematically, we defined α to be the ratio between the rotational velocity, v rot , and the ideal circular velocity in our model grid (i.e., α = v rot /v c ). We defined the relationship between α, radial inflow velocity, and the ideal circular velocity as
. Thus, when α = 1, our model velocities are circular, and when α = 0, the velocities are radially inflowing at the ideal freefall speed of a test particle falling in from infinity. We projected the radial velocity component along the LOS, and added it linearly to the projected LOS rotation velocity at each pixel in our model grid. Upon optimizing, we found that the models strongly favored pure rotation, with a best-fitting value of α = 1. This model converged on M BH = 1.47 × 10 8 M e , leaving the results found for the fiducial model unchanged.
Turbulent velocity dispersion: The ratio of the turbulent velocity dispersion to the rotational velocity (σ/v rot ) determines if the disk can be treated as dynamically cold (where
, or if dynamical pressure effects from turbulence must be accounted for. By using our A H = 0.31 mag MGE model, our radial gas mass profile, the fiducial BH mass, and our best-fit value of σ 0 = 10.5 km s -1 , we find a maximum value of σ 0 /v c = 0.05 (where we have set v rot = v c , the ideal circular velocity) at about r = 50 pc, indicating that treating the disk as dynamically cold and neglecting dynamical pressure effects are justified. Nevertheless, a spatially uniform turbulent velocity dispersion term might be insufficient to characterize possible variations in turbulence across the entire disk. Therefore, in addition to using a spatially uniform gas turbulent velocity dispersion term, σ(r) = σ 0 , we also tried incorporating a Gaussian turbulent velocity dispersion profile r r r exp 2
-into our fiducial model. This profile adds three free parameters and allows for more flexibility in characterizing the overall velocity dispersion. However, the model is not physically motivated, and is used here as a simple tool for exploring possible variations in σ(r). The preferred turbulent velocity dispersion parameters of σ 0 = 10.6 km s -1 , σ 1 = 0.21 km s -1 , r 0 = 0.06 pc, and μ = 0.02 pc yield a turbulent velocity dispersion profile that is dominated by the spatially uniform term of 10.6 km s -1 , and is nearly identical to the fiducial model's spatially uniform σ 0 = 10.5 km s -1 . The modified model yields M BH = 1.47 × 10 8 M e and 1.526
, demonstrating almost no effect on the results of the fiducial model.
Fit region: We tested the sensitivity of M BH to the fit regions of the dynamical models by making two separate adjustments to the spatial fitting ellipse described in Section 4.2.2 and used to calculate χ 2 . The first adjustment was expanding the fitting ellipse to cover the entirety of the gas disk, with semimajor axis length 4 1 and semiminor axis length 1″. For this fitting ellipse, the models converged on M BH = 1.63 × 10 8 M e , an increase of 10.9% from the fiducial model, and 1.451 2 c = n . This larger fitting ellipse includes nearly four times as many data points, but a majority of points are at radii where the extended stellar mass distribution dominates the total enclosed mass.
Our second adjustment was reducing the fitting ellipse to fit only the inner third of the disk, with semimajor and semiminor axis lengths of 1 37 and 0 33. The resultant value of M BH was 1.52 × 10 8 M e , an increase of 3.4% relative to our fiducial model, with 1.337 2 c = n . On this scale, the fit contains a higher fraction of data points that display unresolved gas kinematics, particularly along the disk's minor axis, where beam-smearing effects are more severe.
Pixel oversampling: Other molecular gas-dynamical studies such as Barth et al. (2016a) have found that M BH is relatively insensitive to the choice of pixel oversampling factor, s. We tested our own dynamical models with oversampling factors of s = 1 and s = 4. For s = 1, the result was M BH = 1.45 × 10 8 M e , a decrease of 1.4% relative to model B, with a higher 2 c n of 1.543, as expected for no pixel oversampling. At s = 4, the resulting M BH was 1.47 × 10 8 M e , identical to our fiducial model result, with a slightly improved 2 c n of 1.524. These results show that M BH has little sensitivity to the choice of s, even in the no-oversampling case of s = 1.
Gas mass: We performed two tests to observe the dependence of M BH on the inclusion or exclusion of the gas disk's contribution to the mass model. First, we optimized a dynamical model that did not include the circular velocity contribution from the gas disk, which was derived from the mass surface densities in Boizelle et al. (2017), but otherwise used the same inputs as model B. Without this contribution, the model converged upon M BH = 1.43 × 10 8 M e , a decrease of 2.7%, and the same 1.525 2 c = n as our fiducial model. In addition, we rescaled the gas mass surface densities to the lower value implied by a CO(2-1)/CO(1-0) intensity ratio ≈0.9 (see Section 2.2.2) and reoptimized our dynamical model with this adjusted circular velocity contribution. With this adjustment, our model converged upon M BH = 1.45 × 10 8 M e , and an identical 1.525 2 c = n . These results suggests the inclusion or exclusion of the gas component in models can be important in percent-level precision BH mass measurements. However, for our measurements, it is a relatively minor contribution to the error budget in comparison to that of the dust extinction.
Unresolved active calactic nucleus emission: Given that there is evidence of a weak AGN in NGC 1380, we explored the possibility that our MGEs could be incorporating the light from this AGN in addition to the stars. As a test, we removed the innermost component (FWHM = 0 14 = 11.5 pc) of the A H = 0.31 mag MGE model, and deprojected the remaining components. Using this altered MGE model, the BH mass rose to M BH = 2.04 × 10 8 M e , an increase of 38.7% from the fiducial BH mass (1.47 × 10 8 M e ) found when we included the innermost component, but the reduced χ 2 value also rose to 1.528, indicating a poorer fit to the data. If we deproject this innermost component individually, and assume it is composed entirely of starlight, the corresponding stellar mass is 0.59 × 10 8 M e (assuming ϒ H = 1.42), which is slightly higher (by 2 × 10 6 M e ) than the difference in BH mass. While it is difficult to determine the amount of AGN light in the innermost MGE component, given that the increase in BH mass between the two dynamical models is commensurate with the decrease in assumed stellar masses, this test shows that the BH and stellar mass of the innermost component are degenerate.
The large range of M BH found in Table 2 and the sensitivity tests performed above show that the systematic error is dominated by the uncertainties associated with the host galaxy models. Specifically, these uncertainties are associated with the amount of stellar mass, which changes depending on the assumed dust extinction and/or the presence or absence of an AGN contributing to the central host galaxy light, as shown when we removed the innermost component of our MGE. We chose model B as our fiducial model for the reasons described in Section 5.2 and treat the ΔM BH between its M BH value of 1.47 × 10 8 M e and the maximum (1.85 × 10 8 M e ) and minimum M BH (1.02 × 10 8 M e ) values found in Table 2 as a rough estimate of the uncertainty due to the dust correction. These maximum and minimum values are about 26% larger and 31% smaller than the fiducial M BH value, respectively.
It is clear that the systematic uncertainties exceed the statistical uncertainty (≈1%) and the uncertainty associated with the distance to the galaxy (≈ 8%; Tonry et al. 2001). Considering that the uncertainties from the dust correction and the host galaxy mass modeling dominate the total systematic uncertainty, we adopt the BH mass of 1.47 × 10 8 M e from our fiducial model and the aforementioned uncertainties as our estimate for M BH . Therefore, the expected range of M BH in NGC 1380 is (1.02-2.04) × 10 8 M e .
We optimized two different dynamical models for NGC 6861, which we refer to as models E and F in Table 2. The difference between them is the input stellar circular velocity profile, based on one of the two NGC 6861 MGE models described in Section 3.3 that model the nuclear region in the Hband image with either a Nuker (E) or Core-Sérsic (F) model. These two dynamical models used a uniform turbulent velocity dispersion across the entire disk and are optimized over the annular region described in Section 4.
2.2. Model E converges on M BH = 1.13 × 10 9 M e and ϒ H = 2.52 with 1.987 2 c = n , while model F returns values of M BH = 2.89 × 10 9 M e and ϒ H = 2.14 with 2.004 2 c = n
. The ϒ H values are higher than the ranges predicted by the SSP models of Vazdekis et al. (2010). Similarly, the range of ϒ I = (5.7-6.3) determined observationally by Rusli et al. (2013) is higher than SSP model predictions. All other free parameters remain consistent between the two models, although the inclination angle i slightly increases from 72°.7 to 73°. 6 from model E to F.
We created moment maps and a major axis PVD, and extracted line profiles for model E to compare with the ALMA data. The residuals between the data and model v LOS maps show that our thin disk model emulates the data's observed v LOS well within our designated fitting region, but discrepancies in excess of ∼60 km s -1 are seen at larger radii. These discrepancies highlight kinematic substructure within the disk at these larger radii that our models are unable to reproduce, although given that the maximum value of σ 0 /v c across the NGC 6861 disk is ∼ 0.02, our treatment of the disk as dynamically cold is justified. The PVD also shows discrepancies along the major axis, as structural differences between the data and model PVDs are prevalent at radii larger than ∼3″, which corresponds to the semimajor axis of our elliptical fitting region. The extracted line profiles in Figure 10 show that our models are able to reproduce the observed line profile shapes well, although slight inconsistencies in the peak amplitude, in terms of both overall height and velocity channel, are evident in some spectra.
We conducted Monte Carlo simulations to estimate fitting uncertainties for both model E and F as we did for model B. We generated 150 realizations of each model cube by adding noise to each pixel as described in Section 5.1. For model E, we found a distribution centered at its best-fit M BH value of 1.13 × 10 9 M e with a standard deviation of 4 × 10 7 M e , or 1.4% of M BH . Model F's Monte Carlo simulation was centered at M BH = 2.89 × 10 9 M e and also had a standard deviation of 4 × 10 7 M e , or 3.5% of M BH . We chose to use the standard deviations of each of the free parameters from model E as representative of statistical uncertainties associated with these values and list them in Table 2.
The Monte Carlo simulations show that the statistical modelfitting uncertainties are significantly smaller than the systematic uncertainty associated with our choice of host galaxy model, which return values of M BH that are different by a factor of ∼ 3. Because of this large difference in M BH between models E and F, we did not perform extensive systematic tests on these models as we did for model B in Section 5.2, as the uncertainty associated with our choice of host galaxy model dominates the total error budget.
To determine a lower limit on M BH , we adjusted the central flux in NGC 6861ʼs H-band image in the same manner as was done for NGC 1380ʼs H-band image in Section 3.2. We corrected our Nuker interpolation MGE model under the assumption that the disk resides in the midplane of the galaxy and that only the starlight originating from behind the dust disk experiences any extinction. Our modified Nuker interpolation model explored the extreme limit where A H = ∞ , signifying that all of the light behind the disk is lost, and raising the innermost value of the major axis surface brightness profile by 0.75 mag arcsec 2 -. Even with this maximally peaked surface brightness model, the value of M BH was nonzero and converged on 9.7 × 10 7 M e , which serves as our measurement's lower limit. We emphasize that even with the assumption that the central region of NGC 6861 is optically thick (which is highly unlikely given the J -H color map) and the absence of dynamical information within the inner 1″, our dynamical models still require a central compact mass to reproduce the observations.
Our molecular gas-dynamical measurements are the first and second attempts to determine the masses of the central BHs in NGC 1380 and NGC 6861, respectively. For each galaxy, the presence of dust limits the measurement precision on M BH . In NGC 1380, we find M BH = 1.47 × 10 8 with an uncertainty of ∼40%, which is dominated largely by the dust corrections. The measurement precision for NGC 6861ʼs BH is even more limited due to the lack of dynamical tracers within the BH's sphere of influence, as the resulting values of M BH differ by a factor of ∼3 depending on the model used for the host galaxy. Below, we discuss the importance of resolving the BH's sphere of influence and accounting for the presence of dust in both galaxies. We also discuss how parameter degeneracies emerge within our models from these factors, and we compare our measured M BH values to predictions from the BH-host galaxy scaling relations given by Kormendy & Ho (2013).
The precision of a BH mass measurement is highly dependent on how well the observations resolve the radius of the BH's dynamical sphere of influence, r g . Within r g , the central BH is the dominant contributor to the total gravitational potential. If r g is unresolved, then the measured value of M BH depends heavily on the accuracy of the assumed host galaxy model, and dynamical models are susceptible to parameter degeneracies. For each galaxy, we estimated r g in two distinct ways. The first was to determine the radius where M BH is equal to the enclosed stellar mass, and the second was to calculate r GM g BH 2 s » using our measured values of M BH and known literature values of σ å for both galaxies.
For NGC 1380, we used the best-fit value of M BH = 1.47 × 10 8 M e from our fiducial model to determine r g . The radius where the enclosed stellar mass (derived from the MGE used in the fiducial model) equaled M BH was r g = 18 pc (0 22), which is nearly identical to our average beam size of 0 21. If we instead calculate r g as r GM g BH 2 s = , using the average value of σ å = 215 km s -1 from Hyperleda (Makarov et al. 2014), we find r g = 14 pc (0 17). We display both the total enclosed mass profile and the separate contributions from each component to our dynamical model in Figure 11. These estimates demonstrate it is likely that the observations only marginally resolve r g .
While there is evidence of a slight central upturn in gas velocity within the innermost ∼0 2 of NGC 1380, given that the observations do not fully resolve r g , it is unsurprising that our measurement of M BH carries a large uncertainty of about 40% (dominated mostly by the uncertainties from the dust correction and host galaxy modeling), and that our dynamical models have a degeneracy between BH and stellar mass. In essence, our dynamical models' ability to distinguish their separate contributions is severely limited when r g is not fully resolved. This limitation was demonstrated by our test removal of the innermost component of the host galaxy MGE model, as the increase in M BH was commensurate with the decrease in stellar mass. While the total enclosed mass is well-constrained at large radii, the mass contributions from the BH and the stars within the central regions are not. This degeneracy is highlighted in Figure 5, which shows stellar mass versus radius for our dust-masked and dust-corrected host galaxy models. With each progressive increase in assumed dust extinction, the stellar mass becomes a larger fraction of the total enclosed mass. While the differences in our stellar mass profiles are minimal at radii greater than ∼100 pc, it is their differences within the innermost ∼30 pc that lead to decreases in M BH and variations in ϒ H . Given that the total enclosed mass is tightly constrained by the well-resolved kinematics at large radii, the cuspier surface brightness models require smaller M BH values.
For NGC 6861, we estimated r g using the best-fit values of M BH from both model E and model F. The BH mass and the enclosed stellar mass were equal at 47 pc (0 36) if we adopt the stellar mass profile from model E, and at 94 pc (0 72) for model F. If we instead use the velocity dispersion of σ å = 389 km s -1 from Rusli et al. (2013), we obtain r g = 32 pc (0 25) when adopting the best-fit M BH from model E and 82 pc (0 63) from model F. Considering that the average ALMA beam size for the NGC 6861 data is 0 28, these estimates suggest that r g would be resolved in NGC 6861.
Despite having observations that could in principle resolve r g if the gas disk extended to the center in NGC 6861, the lack of CO emission within the innermost ∼1″ precludes a highprecision BH mass measurement, as our two dynamical models found BH mass values that differed by a factor of ∼3 and were degenerate with ϒ H . We attribute this large disparity to the central hole, and to the differences between our Nuker and Core-Sérsic interpolation MGE models in the dust-affected regions. These differences, especially in the slope of the surface brightness profile, led to distinct ϒ H values in the dynamical models. Additionally, the dearth of CO emission within the innermost 1″ meant that model fits were optimized over pixels that were more sensitive to differences in the stellar mass distribution. Figure 11 shows the separate and combined enclosed mass profiles of the stars, gas, and BH for each host galaxy model used. Using these mass profiles, we determined the total enclosed mass within the central hole (which is divided between contributions from the BH and the stars, due to the absence of gas) to be 7.46 × 10 9 M e (M BH = 1.13 × 10 9 M e ) when using the Nuker interpolation and 8.74 × 10 9 M e (M BH = 2.87 × 10 9 M e ) for the Core-Sérsic interpolation, which are both consistent with results from Boizelle et al. (2017). Considering that the available dynamical information is restricted to radii extending beyond the hole radius, and that our best-fit values of M BH represent a minor fraction of the total dynamical mass within the hole, it is unsurprising that our two dynamical models find very different but ostensibly precise values of M BH . This precision is seen in the results of our Monte Carlo simulations, and it highlights the importance of accounting for these types of systematic uncertainties when making gas-dynamical BH mass measurements in this regime.
Although there is no prior dynamical BH mass measurement for NGC 1380, Pota et al. (2013) did predict M BH = 2.2 10 0.9
1.8 8 -+ based on velocity dispersion measurements of the globular clusters in the galaxy, which is consistent with our findings. Rusli et al. (2013) measured the BH mass in NGC 6861 through stellar-dynamical modeling and found M BH = (2.0 ± 0.2) × 10 9 M e , which is contained within our range of estimates for M BH .
Using Equations ( 6), ( 7), and (8) from Kormendy & Ho (2013), we derived predictions of M BH from the estimated total K-band luminosity, the estimated bulge mass, and the stellar velocity dispersion for each of the galaxies. For the M BH -L bul,K relation, we converted the total K-band apparent magnitudes from 2MASS into corresponding total K-band luminosities using our assumed luminosity distances. For NGC 1380, we adopted the measured R-band B/T = 0.359 from Gao et al. (2019) to determine L bul,K , while we used B/T = 1 for NGC 6861 adopting its classification as an elliptical galaxy. Schulz et al. (2003) studied the wavelength dependence of B/T in ETGs through evolutionary synthesis modeling, and found that while B/T does change substantially from the U through I bands, the changes diminish at redder wavelengths. In addition, Li et al. (2011) determined that the observed B -I, V -I, and R -I profiles in NGC 1380 remained flat over a large radial range. Based on this information, we expect that the difference between R and K-band B/T values are relatively small, and given that there are no presently available B/T measurements in the K-band, we use the R-band value as an estimate in our calculations of L bul,K . To compare with the BH -M bul relation, we calculated L bul,H for each galaxy by assuming H -K = 0.2 mag based on SSP models and an absolute H-band (K-band) magnitude of 3.37 (3.27) for the Sun (Willmer 2018), and multiplied L bul,H by our best-fit ϒ H values to derive an estimate for M bul . Finally, to compare our BH mass measurements to the M BH -σ å relation, we used the σ å values of 215 km s -1 and 389 km s -1 for NGC 1380 and NGC 6861, respectively. We note that the M BH -L bul,K , M BH -M bul , and M BH -σ å relations of Kormendy & Ho (2013) have intrinsic scatters of 0.28, 0.30, and 0.28 dex.
With a best-fit value of M BH = 1.47 × 10 8 M e and an associated uncertainty of about 40%, our estimate of M BH in NGC 1380 generally agrees with predictions made by the BHhost galaxy scaling relations. For NGC 1380, we derived a total K-band bulge luminosity of 3.8 × 10 10 L e and a bulge mass of M bul = 5.0 × 10 10 M e to use in the M BH -L bul,K and M BH -M bul relations. These relations predict ranges of M BH = (1.4-2.0) × 10 8 M e and M BH = (1.8-2.5) × 10 8 M e , respectively. The range predicted from the stellar velocity dispersion of 215 km s -1 is (3.7-4.8) × 10 8 M e . Thus, our measurement of M BH in NGC 1380 directly overlaps with the predictions made from the M BH -L bul,K and M BH -M bul relations, and lies slightly below and outside the scatter from the M BH -σ å relation.
For NGC 6861, our range of M BH = (1-3) × 10 9 M e is consistent with the previous measurement of M BH = (2.0 ± 0.2) × 10 9 M e by Rusli et al. (2013), and generally agrees with predictions from the BH-host galaxy scaling relations. We derived a total K-band bulge luminosity of 1.2 × 10 11 L e , which estimates M BH = (0.6-0.8) × 10 9 M e . Our estimated bulge mass of M bul = 2.8 × 10 11 M e is higher than the value provided in Kormendy & Ho (2013) of M bul = 1.8 × 10 11 M e . Using our value, the predicted range of M BH is (1.4-2.0) × 10 9 M e , while using the lower Kormendy & Ho (2013) estimate of M bul gives (0.8-1.1) × 10 9 M e . The stellar velocity dispersion of 389 km s -1 in NGC 6861 is one of the highest measured in an ETG. It predicts (5.4-11.4) × 10 9 M e from the M BH -σ å relation, which is higher than our measured value, although our measurement is still contained within the relation's intrinsic scatter. Machacek et al. (2010) suggested that NGC 6861 may have had strong gravitational encounters in its past with neighboring galaxies that has elevated its central velocity dispersion, and that it could be the dominant galaxy in a galaxy subgroup that is merging based on Chandra X-ray observations. While this possibility has yet to be confirmed, our measured range of M BH and the measurement by Rusli et al. (2013) suggest that the M BH -σ å relation slightly overpredicts M BH in NGC 6861. Whether this excess should be attributed to larger intrinsic scatter at the high-σ å end of the relation or physical mechanisms that have affected the growth and evolution of NGC 6861 and its BH remains unclear.
We present gas-dynamical measurements of the BH masses in NGC 1380 and NGC 6861 using ALMA CO(2-1) observations at 0 21 and 0 28 resolution, respectively. We find evidence for gas disks exhibiting regular rotation in the central regions of both galaxies. For NGC 1380, a slight central increase is observed in its maximum LOS velocity, reaching approximately ±300 km s -1 relative to the systemic velocity of the galaxy, as expected of rotation around a central BH. In NGC 6861, the presence of a rotating gas disk with ring-like structure is observed with peak LOS velocities of ∼ 500 km s -1 , but a central ∼1″ hole in the CO distribution precludes a precise measurement of M BH .
For NGC 1380, we determine M BH = 1.47 × 10 8 M e with an uncertainty of about 40% by optimizing thin disk models to the ALMA observations with four different host galaxy models. We find that our measured values of M BH are degenerate with the enclosed stellar mass, and that the uncertainties associated with the dust corrections and host galaxy models dominate the error budget. Given the slight central rise in observed LOS velocity, it is possible that higher-resolution ALMA observations could provide a more confident determination of the BH mass by lifting the stellar and BH mass degeneracy.
In the case of NGC 6861, we optimize dynamical modeling fits to the ALMA CO data using two different host galaxy models, and find that the results for M BH differ by a factor of ∼ 3 due to the lack of dynamical tracers within the innermost 1″, and to the structural differences in the shape of the dustcorrected surface brightness and stellar mass profiles of the host galaxy. Given the large difference between the two results, the value of M BH in NGC 6861 cannot be precisely constrained, although we find that our models suggest a plausible range of M BH = (1-3) × 10 9 M e and a lower limit of ∼1 × 10 8 M e derived by assuming an unlikely amount of central dust extinction. This range encompasses the stellar-dynamical mass measurement of (2.0 ± 0.2) × 10 9 M e determined by Rusli et al. (2013).
When comparing our measurements of M BH to the BH-host galaxy scaling relations determined by Kormendy & Ho (2013), we find that our measured values of M BH for NGC 1380 and NGC 6861 are generally consistent with the M BH -L bul and M BH -M bul relations, and are below the expected values predicted by the M BH -σ å relation. Although for NGC 6861, the measured value of M BH remains within this relation's scatter. More precise BH mass measurements on the high-mass end of these scaling relations are needed to further our understanding of the differences among them and to determine whether there is more intrinsic scatter than previously thought.
Our work highlights a number of factors that limit gasdynamical BH mass measurements with ALMA. Factors such as dust obscuring the stellar light, a lack of high-velocity emission within r g , and the presence of a hole in the CO distribution lead to degeneracies among model parameters and large systematic uncertainties that are important to account for.
A key area of improvement would to incorporate realistic 3D radiative transfer modeling codes (De Geyter et al. 2013;Camps & Baes 2015) to recover the intrinsic stellar surface brightness of the host galaxy from NIR images. This goal is especially important for ALMA observations of dusty ETGs that do not resolve gas deep within the sphere of influence. Nevertheless, ALMA observations in this regime provide highresolution information on circumnuclear disks in ETGs and meaningful constraints on M BH . These BH mass measurements will continue to add valuable information to both local BH demographics and our understanding of galaxy evolution. Note. The NGC 6861 A H = 0.00 mag (Core-Sérsic Model) MGE solution built from the HST H-band image. This MGE solution was used to optimize dynamical model F. The first column is the component number, the second is the central surface brightness corrected Galactic extinction and assuming an absolute solar magnitude of M e,H = 3.37 mag (Willmer 2018), the third is the Gaussian standard deviation along the major axis, and the fourth is the axial ratio, which was constrained to have a minimum value of 0.400 to allow for a broader range in the inclination angle during the deprojection process. Primes indicate projected quantities. Note. Additional NGC 1380 MGE solutions built from the HST H-band image. These MGEs were used to optimize dynamical models A, B, and D. These solutions have assumed central extinction values of A H = 0.00, 0.75, and 1.50 mag. The first column is the component number, the second is the central surface brightness corrected for Galactic extinction and assuming an absolute solar magnitude of M e,H = 3.37 mag (Willmer 2018), the third is the Gaussian's standard deviation along the major axis, and the fourth is the axial ratio, which was constrained to have a minimum value of 0.400 to allow for a broader range in the inclination angle during the deprojection process. Primes indicate projected quantities.
The Astrophysical Journal, 934:162 (21pp), 2022 August 1 Kabasares et al.