Black Hole Mass Measurements of Early-type Galaxies NGC 1380 and NGC 6861 through ALMA and HST Observations and Gas-dynamical Modeling^*

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 distortion-corrected image with a pixel scale of 008 pixel⁻¹. 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 008 pixel⁻¹ 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.

2.2.1. Observations and Data Reduction

2.2.2. Circumnuclear Disk Properties

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⁻¹ and ∼500 km s⁻¹ 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 52 (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 ∼01, 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_⊙ pc⁻² (K km s⁻¹)⁻¹ (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⁷ M_⊙ and (25.6 ± 8.9) × 10⁷ M_⊙ 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 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.

Zoom In Zoom Out Reset image size

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 55 in NGC 6861.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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-to-one 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 .

Zoom In Zoom Out Reset image size

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⁻². 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 dust-masked and dust-corrected MGEs to model the host galaxy’s light following methods used for similar galaxies by Boizelle et al. (2019, 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 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, one-half, 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 low-extinction 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 = 25 (≈ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  H-band MGE Parameters

k	log10 IH,k (L⊙ pc−2)	(arcsec)		log10 IH,k (L⊙ pc−2)	(arcsec)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
	NGC 1380			NGC 6861
	AH = 0.31 mag			AH = 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	⋯	⋯	⋯	1.662	33.017	0.641
15	⋯	⋯	⋯	1.684	52.142	0.883
16	⋯	⋯	⋯	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 χ². 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_⊙,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.

Download table as:  ASCIITypeset image

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 55 (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.

Zoom In Zoom Out Reset image size

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 = 031(≈ 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 H-band 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 =35(≈455 pc), r_e = 79(≈1 kpc), γ = 0.61, and α = 2.13. We replaced the pixels in the dust disk region in the original H-band 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 H-band 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.

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. 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) = σ₀, 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 flux-scaling factor F₀, 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 when deriving v_⋆,MGE for each galaxy. At each model iteration, 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⁷ M_⊙ and (25.6 ± 8.9) × 10⁷ M_⊙ 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⁶–10⁷ M_⊙, 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 scikit-image 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².

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 χ² 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 χ² 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. Noise Model

4.2.2. Fit Regions

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.

Table 2. Dynamical Modeling Results

Model	MGE	MBH	ϒH	i	Γ	σ0	xc	yc	vsys	F0
			(M⊙/L⊙)	(°)	(°)	(km s−1)	(arcsec)	(arcsec)	(km s−1)
NGC 1380
	(AH mag)	(108 M⊙)
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)	(109 M⊙)
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.

Download table as:  ASCIITypeset image

Models A–D yield best-fit values of M_BH in the range of (1.02–1.85) × 10⁸ M_⊙, ϒ_H between 1.30 and 1.42, and a range in reduced χ² 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 ∼05. Moment 1 (v_LOS) maps show that data and model velocity fields are also in good agreement, although differences of ∼30 km s⁻¹ are noticeable toward the disk edge, outside the fit region, and at the disk center, where the impact of beam-smearing 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 block-averaged 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⁸ M_⊙ with a standard deviation of 2 × 10⁶ M_⊙, 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⁸ M_⊙. 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⁸ M_⊙, representing a ∼20% decrease from the initial fit. This particular model also shows a decrease in the resultant , as model A with the initial MGE returns , while model B with A_H = 0.31 mag yields . Increasing the extinction further to A_H = 0.75 mag and A_H = 1.50 mag further decreases M_BH, as models C and D converge to values of 1.27 × 10⁸ M_⊙ and 1.02 × 10⁸ M_⊙, i.e., ∼31% and ∼45% decreases from the model A value, respectively. However, both models C and D result in higher 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 χ² 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.

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⁸ M_⊙, 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 σ₀ = 10.5 km s⁻¹, we find a maximum value of σ₀/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) = σ₀, we also tried incorporating a Gaussian turbulent velocity dispersion profile 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 σ₀ = 10.6 km s⁻¹, σ₁ = 0.21 km s⁻¹, r₀ = 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⁻¹, and is nearly identical to the fiducial model’s spatially uniform σ₀ = 10.5 km s⁻¹. The modified model yields M_BH = 1.47 ×10⁸ M_⊙ and , 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 χ². The first adjustment was expanding the fitting ellipse to cover the entirety of the gas disk, with semimajor axis length 41 and semiminor axis length 1″. For this fitting ellipse, the models converged on M_BH = 1.63 × 10⁸ M_⊙, an increase of 10.9% from the fiducial model, and . 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 137 and 033. The resultant value of M_BH was 1.52 × 10⁸ M_⊙, an increase of 3.4% relative to our fiducial model, with . 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⁸ M_⊙, a decrease of 1.4% relative to model B, with a higher of 1.543, as expected for no pixel oversampling. At s = 4, the resulting M_BH was 1.47 × 10⁸ M_⊙, identical to our fiducial model result, with a slightly improved 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⁸ M_⊙, a decrease of 2.7%, and the same 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⁸ M_⊙, and an identical . 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 = 014 = 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⁸ M_⊙, an increase of 38.7% from the fiducial BH mass (1.47 × 10⁸ M_⊙) found when we included the innermost component, but the reduced χ² 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⁸ M_⊙ (assuming ϒ_H = 1.42), which is slightly higher (by 2 × 10⁶ M_⊙) 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⁸ M_⊙ and the maximum (1.85 × 10⁸ M_⊙) and minimum M_BH (1.02 × 10⁸ M_⊙) 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⁸ M_⊙ 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⁸ M_⊙.

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 H-band 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⁹ M_⊙ and ϒ_H = 2.52 with , while model F returns values of M_BH =2.89 × 10⁹ M_⊙ and ϒ_H = 2.14 with . 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 727 to 736 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⁻¹ 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 σ₀/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.

Zoom In Zoom Out Reset image size

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⁹ M_⊙ with a standard deviation of 4 × 10⁷ M_⊙, or 1.4% of M_BH. Model F’s Monte Carlo simulation was centered at M_BH = 2.89 × 10⁹ M_⊙ and also had a standard deviation of 4 × 10⁷ M_⊙, 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 model-fitting 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 . Even with this maximally peaked surface brightness model, the value of M_BH was nonzero and converged on 9.7 × 10⁷ M_⊙, 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.

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 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⁸ M_⊙ 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 (022), which is nearly identical to our average beam size of 021. If we instead calculate r_g as , using the average value of σ_⋆ = 215 km s⁻¹ from Hyperleda (Makarov et al. 2014), we find r_g = 14 pc (017). 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 .

Zoom In Zoom Out Reset image size

While there is evidence of a slight central upturn in gas velocity within the innermost ∼02 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 (036) if we adopt the stellar mass profile from model E, and at 94 pc (072) for model F. If we instead use the velocity dispersion of σ_⋆ = 389 km s⁻¹ from Rusli et al. (2013), we obtain r_g = 32 pc (025) when adopting the best-fit M_BH from model E and 82 pc (063) from model F. Considering that the average ALMA beam size for the NGC 6861 data is 028, 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 high-precision 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⁹ M_⊙ (M_BH = 1.13 ×10⁹ M_⊙) when using the Nuker interpolation and 8.74 ×10⁹ M_⊙ (M_BH = 2.87 × 10⁹ M_⊙) 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 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⁹ M_⊙, 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 M_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⁻¹ and 389 km s⁻¹ 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⁸ M_⊙ and an associated uncertainty of about 40%, our estimate of M_BH in NGC 1380 generally agrees with predictions made by the BH–host galaxy scaling relations. For NGC 1380, we derived a total K-band bulge luminosity of 3.8 × 10¹⁰ L_⊙ and a bulge mass of M_bul = 5.0 × 10¹⁰ M_⊙ 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⁸ M_⊙ and M_BH = (1.8–2.5) × 10⁸ M_⊙, respectively. The range predicted from the stellar velocity dispersion of 215 km s⁻¹ is (3.7–4.8) × 10⁸ M_⊙. 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⁹ M_⊙ is consistent with the previous measurement of M_BH =(2.0 ± 0.2) × 10⁹ M_⊙ 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¹¹ L_⊙, which estimates M_BH = (0.6–0.8) × 10⁹ M_⊙. Our estimated bulge mass of M_bul = 2.8 × 10¹¹ M_⊙ is higher than the value provided in Kormendy & Ho (2013) of M_bul = 1.8 × 10¹¹ M_⊙. Using our value, the predicted range of M_BH is (1.4–2.0) × 10⁹ M_⊙, while using the lower Kormendy & Ho (2013) estimate of M_bul gives (0.8–1.1) × 10⁹ M_⊙. The stellar velocity dispersion of 389 km s⁻¹ in NGC 6861 is one of the highest measured in an ETG. It predicts (5.4–11.4) ×10⁹ M_⊙ 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.