Gas-dynamical Mass Measurements of the Supermassive Black Holes in the Early-type Galaxies NGC 4786 and NGC 5193 from ALMA and HST Observations

We present molecular gas-dynamical mass measurements of the central black holes in the giant elliptical galaxies NGC 4786 and NGC 5193, based on CO (2−1) observations from the Atacama Large Millimeter/submillimeter Array (ALMA) and Hubble Space Telescope near-infrared imaging. The central region in each galaxy contains a circumnuclear disk that exhibits orderly rotation with projected line-of-sight velocities of ∼270 km s−1. We build gas-dynamical models for the rotating disk in each galaxy and fit them directly to the ALMA data cubes. At 0.″31 resolution, the ALMA observations do not fully resolve the black hole sphere of influence (SOI), and neither galaxy exhibits a central rise in rotation speed, indicating that emission from deep within the SOI is not detected. As a result, our models do not tightly constrain the central black hole mass in either galaxy, but they prefer the presence of a central massive object in both galaxies. We measure the black hole mass to be (MBH/108M⊙)=5.0±0.2[1σstatistical]−1.3+1.4[systematic] in NGC 4786 and (MBH/108M⊙)=1.4±0.03[1σstatistical]−0.1+1.5[systematic] in NGC 5193. The largest component of each measurement’s error budget is from the systematic uncertainty associated with the extinction correction in the host galaxy models. This underscores the importance of assessing the impact of dust attenuation on the inferred M BH.


INTRODUCTION
Supermassive black holes (BHs) are thought to reside at the centers of most, if not all, massive galaxies.With masses between a million to over a billion times that of the Sun, supermassive BHs gravitationally dominate the orbits of objects within their sphere of influence (SOI).The radius of the SOI is often defined as either r SOI ≈ GM BH /σ 2 ⋆ , where σ ⋆ represents the stellar velocity dispersion of the spheroidal component of the galaxy, or the radius where the enclosed galaxy and BH mass are equal M ⋆ (r SOI ) = M BH .While the SOI of a supermassive BH is limited to the innermost part of a galaxy, research has identified scaling relations between the mass of the central supermassive BH and large-scale properties of the host galaxy such as the spheroidal component's stellar velocity dispersion, luminosity, and mass.These relations suggest that galaxies and their supermassive BHs coevolve through cosmic time and regulate each other's growth (Gültekin et al. 2009;Kormendy & Ho 2013;McConnell & Ma 2013;Saglia et al. 2016).
Increasing the sample of reliably measured BH masses is vital to understanding the scaling relations and BHhost galaxy coevolution as a whole.Over the last three decades, there have been approximately 100 supermassive BH mass measurements obtained primarily through either stellar or ionized gas-dynamical modeling (Kormendy & Ho 2013).A key factor in securing a robust dynamical BH mass measurement is using observations that can probe scales near or within the projected radius of the BH's SOI.For extragalactic sources, H 2 O megamaser emission can be observed at very high angular resolution through Very Long Baseline Interferometry and have been used to secure some of the most precise extragalactic BH masses to date (Miyoshi et al. 1995;Kuo et al. 2011Kuo et al. , 2018)), but given that they are very rare (van den Bosch et al. 2016), other tracers must be used to understand BH demographics.
Since it became operational within the past decade, the Atacama Large Millimeter/submillimeter Array (ALMA) has been used to observe the rotation of molecular gas disks to constrain BH masses in nearby galaxies.In the best cases, ALMA has provided BH mass measurements with uncertainties of ≤10% (Barth et al. 2016a;Boizelle et al. 2019;North et al. 2019).In this paper, we add to the growing number of BH mass measurements in early-type galaxies (ETGs) with ALMA (Barth et al. 2016a;Davis et al. 2017;Onishi et al. 2017;Davis et al. 2018;Smith et al. 2019Smith et al. , 2021;;Boizelle et al. 2021;Cohn et al. 2021;Kabasares et al. 2022;Ruffa et al. 2023) by dynamically measuring the masses of the central BHs in the ETGs NGC 4786 and NGC 5193.
These two galaxies were selected as part of an ALMA program designed to identify rapidly rotating molecular gas on scales comparable to r SOI .The selection was based on the identification of a smooth and morphologically round circumnuclear dust disk from optical Hubble Space Telescope (HST) images.ALMA CO(2−1) observations presented by Boizelle (2018) revealed that their molecular gas kinematics are dominated by disklike rotation, which we use to constrain the masses of their central BHs.The optical dust disks in both galaxies are relatively small.In angular size, the dust disk in NGC 4786 has a radius of about 0. ′′ 6 or 181 pc and the disk in NGC 5193 has a radius of about 1 ′′ or 221 pc given our assumed angular diameter distances (discussed further in this section).The projected line-of-sight (LOS) velocities of the molecular gas in both disks exhibit similar features as they are in excess of ∼270 km s −1 in the outer parts of the disk and remain somewhat flat at radii extending toward the disk center.
NGC 4786 is classified as a cD pec galaxy in the Third Reference Catalog of Bright Galaxies (RC3; de Vaucouleurs et al. 1991).Redshift-independent distances for this galaxy in the NASA/IPAC Extragalactic Database 1 (NED) are between 65 − 75 Mpc when using a ΛCDM cosmology with H 0 = 67.8km s −1 Mpc −1 .We adopt a Hubble-flow distance, using a value of H 0 = 73 km s −1 Mpc −1 based on more recent estimates of H 0 from nearby (< 100 Mpc) galaxies (Blakeslee et al. 2021;Riess et al. 2022;Kenworthy et al. 2022), a recessional velocity of cz = 4623 km s −1 from preliminary gas-dynamical models, Ω M = 0.31, and Ω Λ = 0.69.These assumptions set a luminosity distance of 64.1 Mpc, an angular diameter distance of 62.1 Mpc, and an angular scale where 1 ′′ corresponds to 301 pc.The measured M BH scales linearly with the assumed distance, so any differences in assumed distance to the galaxy will correspond to an equivalent rescaling of the measured M BH .There are no previous studies that have measured the mass of the supermassive BH in this galaxy.
NGC 5193 is classified as an E pec galaxy by RC3.In past literature, this galaxy has been associated with and sometimes identified as the brightest cluster galaxy in the Abell 3560 (Abell et al. 1989) cluster with a recessional velocity of cz ∼ 3800 − 4000 km s −1 (Lauer et al. 1998;Okoń & Harris 2002).However, Vettolani et al. (1990) and Willmer et al. (1999) instead associate NGC 5193 with the Abell 3565 cluster based on findings from Melnick & Moles (1987) that indicate Abell 3560 has a recessional velocity of cz ∼ 14,850 km s −1 .This implies that Abell 3560 is a background galaxy cluster and that NGC 5193 is not a member.Tonry et al. (2001) determined a distance modulus of m − M = 32.66 ± 0.29 mag for NGC 5193 with ground-based I-band surface brightness fluctuation (SBF) data while distinguishing it as separate from Abell 3560.This distance modulus translates to a luminosity distance of 34.0±4.5 Mpc, but they list this measurement as uncertain.Jensen et al. (2003) independently measured a distance modulus of m − M = 33.35± 0.15 mag, using SBF measurements from HST NICMOS F160W data.The corresponding luminosity distance is 46.8 ± 3.2 Mpc, and we adopt this distance for our dynamical modeling purposes.Using a recessional velocity of cz = 3705 km s −1 from initial gasdynamical models, this gives an angular diameter distance of 45.7 Mpc, where 1 ′′ corresponds to 221 pc.As in the case of NGC 4786, there are no previous works that have constrained the central BH mass in this galaxy.
We organize our paper as follows.In Section 2, we describe our ALMA and HST observations as well as the data calibration and reduction process.In Section 3, we model the observed 2D surface brightness distributions of each galaxy and build host galaxy models that account for dust extinction.A description of our dynamical modeling formalism is described in Section 4, and the results are presented in Section 5. We discuss the challenges and limitations of our measurements in Section 6, and conclude in Section 7.  Mullin et al. 2007), and then imaged into a data cube with 20 km s −1 velocity channel spacings (with respect to the rest frequency of the CO(2−1) emission line at 230.538 GHz) using a robust parameter of 0.5.We chose a pixel size of 0. ′′ 05 to adequately sample the synthesized beam's minor axis.The beam's position angle is 67.3 • measured East of North.The major axis full width at half maxmium (FWHM) of the synthesized beam is 0. ′′ 35, whereas the minor axis has a FWHM of 0. ′′ 27, giving it a geometric mean FWHM of 0. ′′ 31.
NGC 5193 was observed on 15 January 2018 for approximately 29 minutes with a maximum baseline of 2386 m.The observation targeted the redshifted CO(2−1) emission line along with a corresponding spectral window for the continuum.The emission line spec-tral window had a channel resolution of 3.904 MHz, and covered the frequency range of 226.84−228.71GHz.The continuum windows had a channel resolution of 31.25 MHz, GHz.The uv-plane visibilities were calibrated in CASA version 5.1.1,and then imaged into a data cube with 10 km s −1 velocity channel spacings, with a pixel size of 0. ′′ 035.The synthesized beam has a position angle of 63.1 • measured East of North, has a major axis FWHM of 0. ′′ 33, and has a minor axis FWHM of 0. ′′ 29, giving it a geometric mean FWHM of 0. ′′ 31.

Circumnuclear Disk Properties
A detailed description of the properties of these ALMA datasets was presented by Boizelle (2018).We briefly describe some of their properties below.As seen in Figure 1, the CO emission is cospatial with the optical dust disk.The CO surface brightness extends about 0. ′′ 65 in radius along the disk's major axis in NGC 4786 and about 1. ′′ 2 for NGC 5193.Both disks display orderly rotation about their centers, with the projected LOS velocities reaching ∼270 km s −1 in each disk.Examination of their respective position-velocity diagrams (PVDs) extracted along the major axis reveals that the CO velocities are relatively flat and decrease slightly towards their respective disk center.The PVDs also highlight a lack of CO-bright gas well within the expected BH SOIs.
To incorporate the mass of the gas disk in the total gravitational potential of our dynamical models, we convert the integrated CO(2−1) flux measurements over each disk into M gas profiles.These M gas profiles are dominated by molecular hydrogen and helium and are calculated as M gas = M H2 (1 + f He ), where we set f He = 0.36.
The process of generating M gas profiles starts with the construction of an integrated CO(2−1) flux map.To build this map, we multiply the data cube with a 3D mask generated by the 3DBarolo program (Di Teodoro & Fraternali 2015) and sum the channel maps along the spectral axis to generate a 2D map of integrated flux density.Upon converting the map into units of integrated flux, we average the flux on elliptical annuli centered on the disk centers.If we sum the integrated flux across the entire region of each disk, we find I CO(2−1) = (6.90± 0.14) Jy km s −1 for NGC 4786 and I CO(2−1) = (7.10 ± 0.05) Jy km s −1 for NGC 5193.These statistical uncertainties are calculated through Monte Carlo simulations, but there is an additional 10% systematic uncertainty that stems from the flux scale.For each elliptical annulus, the integrated CO(2−1) flux measurements are converted into CO(1−0) luminosities using: (Carilli & Walter 2013) assuming a CO(2−1)/CO(1−0) line ratio of 0.7 (Lavezzi et al. 1999).Then, a mass of H 2 is obtained by multiplying the CO(1−0) luminosities by α CO = 3.1 M ⊙ pc −2 (K km s −1 ) −1 (Sandstrom et al. 2013) and then multiplying this result by 1.36 as described above to generate an estimate of M gas , though it should be noted that the most appropriate α CO value for ETGs is unknown, and thus the estimated M gas values should be taken as approximations.We find M gas = 6.9 × 10 7 M ⊙ for the disk in NGC 4786 and M gas = 3.9 × 10 7 M ⊙ in NGC 5193.Along with the unknown ideal value of α CO , the uncertainty in D L contributes ≈15% based on the range of redshiftindependent distances in NED for NGC 4786 and ≈7% from the SBF distance measurement to NGC 5193 by Jensen et al. (2003).This distance uncertainty is in addition to the 10% systematic and (1 − 2)% statistical uncertainty associated with the integrated flux measurements.As will be discussed in Section 5.3, the inclusion or exclusion of the gas mass in the total gravitational potential of the system only contributes a small amount to the total error budget on M BH in each galaxy.

HST Observations
All the HST data used in this paper can be found in MAST: 10.17909/gf2w-xv03.For each galaxy, we retrieved F110W (J-band) and F160W (H-band) images taken with the Wide Field Camera 3 (WFC3).For NGC 4786, the H-band image was taken with a fourpoint dither pattern with four separate exposures that lasted 249 seconds each.The H-band images for NGC 5193 employed a similar observation strategy with four separate exposures lasting 299 seconds each.For the Jband images, a two-point dither pattern was used and two separate 249-second exposures were taken for NGC 4786, whereas two 128-second exposures were taken for NGC 5193.Further details regarding the data mosaicing process, and construction of dust-masked MGEs for a larger sample of ETGs will soon be available (Davidson et al., in prep), but we summarize the key steps below.
We processed our data through the calwf3 pipeline and used AstroDrizzle to combine and align the separate exposures.To start, we drizzled and aligned the flat-fielded H-band images to a pixel scale of 0. ′′ 08, and used this image as the reference image when drizzling and aligning the J-band image.We determined offsets between the H and J-band images from the luminosityweighted galaxy center coordinates of each image, and interpolated to align them to within ∼0.2 subpixels of accuracy based on inspection of the J − H maps. Additionally, we constructed TinyTim (Krist & Hook 2004) model point-spread functions (PSFs) that were dithered and drizzled in the same fashion as the H-band image.These PSFs are needed to construct the host galaxy models.

HOST GALAXY MODELING
An accurate host galaxy model that accounts for the mass of the galaxy's stars is crucial when measuring the mass of a central BH.These host galaxy models are constructed by modeling and deprojecting the observed 2D surface brightness profiles from the drizzled H-band images.Given the presence of a circumnuclear dust disk, we chose to model the galaxy's surface brightness in the H-band because dust attenuation is reduced at longer wavelengths.We used the Multi-Gaussian Expansion (MGE) method (Emsellem et al. 1994;Cappellari 2002) to parameterize the surface brightness of both galaxies with a sum of concentric Gaussian functions.

Dust-Masked and Unmasked MGE Models
We built three unique MGE models for NGC 4786 and NGC 5193 following an approach similar to that used by Cohn et al. (2021).The three models correspond to an unmasked, dust-masked, and dust-corrected MGE model that account for the effects of dust differently, and are used to assess the systematic impact of the chosen MGE model on our BH mass measurement.
As seen in Figure 2, the H-band dust attenuation does not appear severe in either NGC 4786 and NGC 5193, hence we first explored both unmasked and dust-masked MGE models to explore the impact on the measurement of M BH .Prior to fitting MGEs to the drizzled H-band image of each galaxy, we isolated the host galaxy light in each image from extraneous sources such as neighboring galaxies, foreground stars, cosmic rays, and detector artifacts by masking these objects.In addition, we constructed J − H maps for each galaxy in order to better identify and mask out regions where dust obscuration was highest, typically corresponding to areas where J −H > 0.88 mag, or equivalently where the color excess is 0.08 mag higher than in regions just outside the disk.For all reported H and J-band magnitudes in this work, we use the Vega magnitude system.
We first modeled the 2D surface brightness with routines from the MgeFit package in Python (Cappellari 2002).The components of this initial MGE were then used as initial guesses for a second MGE fit using the GALFIT program (Peng et al. 2002).We chose to use GALFIT for our final MGEs because it allows for an asymmetric 2D PSF to be incorporated in the modeling process, in contrast to MgeFit which requires decomposing the PSF into a sum of circular Gaussian functions when used in the MGE construction.For both programs, we accounted for the blurring due to the Hband PSF by incorporating the TinyTim H-band PSF models we built.In addition, our MGEs account for a foreground H-band Galactic reddening of A H = 0.019 mag for NGC 4786 and A H = 0.029 mag for NGC 5193 (Schlafly & Finkbeiner 2011).Our MGE solutions contain fourteen Gaussian components for NGC 4786 and eight Gaussian components for NGC 5193 and are listed in Tables 4 and 5 in the Appendix.
We carefully considered whether to include a PSF component in GALFIT to account for a potential unresolved nuclear source of non-stellar origin.Optical spectra of the nuclear regions showed no evidence of prominent emission lines typically associated with an active galactic nucleus in either galaxy (Jones et al. 2009).While the NGC 5193 H-band surface brightness profile exhibits a cuspy nature, an examination of the galaxy center in multiple wavelength filters revealed that the central light is radially extended and not pointlike.Based on these findings, we decided not to include a central PSF component in our model for either galaxy's central surface brightness distribution.
Our MGE models assume that the galaxy has an oblate and axisymmetric shape, and that each Gaussian component shares the same center and position angle.While individual Gaussian components may not correspond to physically distinct galaxy components, their projected axial ratios q ′ k may converge on values below cos(i) in the MGE optimization process, where i is the inclination of the cirumnuclear disk.A useful proxy for the inclination is i = arccos(b/a) where b/a is the observed axial ratio of the disk as measured from the HST images.This proxy typically agrees with kinematic inclinations derived from dynamical models to within ∼5 • based on previous studies (Boizelle et al. 2019(Boizelle et al. , 2021;;Cohn et al. 2021;Kabasares et al. 2022), and so we set a lower bound on the possible MGE component axial ratios of 0.69 and 0.75 for NGC 4786 and NGC 5193.This enables deprojection of the MGEs down to inclinations as low as 46 • and 51 • , respectively.
Examination of the 2D isophotes in Figure 2 shows that the model isophotes are an excellent match to those seen in the H-band data out to ∼100 ′′ .Within the central dusty regions, the observed H-band isophotes for both NGC 4786 and NGC 5193 remain relatively symmetrical and are modeled well by their dust-masked MGEs.Extracting major axis surface brightness profiles from both the MGE model and the H-band data

Dust-Corrected MGE Models
The final MGE we created for each galaxy involved fitting an MGE model to a dust-corrected H-band image.The process of developing this MGE model follows methods described by Boizelle et al. (2019Boizelle et al. ( , 2021)), Cohn et al. (2021), and Kabasares et al. (2022).We summarize the key steps below.
We fit a 2D Nuker model (Faber et al. 1997) to the innermost 10 ′′ × 10 ′′ of the drizzled H-band image using GALFIT.This fit included the mask we used for the dust-masked MGE, and acts as the starting point of our dust correction.We use a 2D Nuker model to fit the central region of the galaxy as we can tune its parameters to create dust-corrected HST images corresponding to specific values of A H at the center.We prefer using Nuker models to perform these adjustments as op- The surface brightness measurements are made with the Python-implementation of the sectors photometry routine (Cappellari 2002) which performs photometry along evenly spaced sectors from the major axis to the minor axis and averages measurements over the four quadrants of the image.(Bottom) A comparison of the radial ellipticity of the observed H-band photometry and the 2D MGE models for NGC 4786 (left) and NGC 5193 (right).We use the photutils.isophoteroutine (Jedrzejewski 1987;Bradley et al. 2023) which fits elliptical isophotes to a galaxy image to determine the ellipticity of each isophote along the semi-major axis.For each panel, the red squares are the observed values from the H-band image, while blue squares are dust-corrected values described in Sections 3.2.The solid lines in each panel correspond to profiles extracted along the major axis for each of our 2D MGE models.Red lines are for dust-masked MGE models whereas black and blue lines represent dust-unmasked and dust-corrected MGEs, respectively.The dashed lines indicate the dust disk edge and the arrows indicate that the dust extends down to the nucleus.
posed to simply adjusting the observed flux values in the HST image and fitting an MGE to this adjusted image directly because the Nuker models not only provide a way to handle the adjustments at the center, but also provide an estimate of the intrinsic surface brightness in dust-affected regions.This is in contrast to simply adjusting the observed surface brightness values at the galaxy center in the HST image and fitting an MGE to this new image, as the adjustment can create sharp or discontinuous features in the surface brightness profile that are unphysical and can manifest in the resultant MGE model.The Nuker model fits in GALFIT include the H-band PSFs, so the resulting solutions correspond to intrinsic parameters.Nuker models have been shown to accurately model the surface brightness distribution within the innermost few arcseconds of early-type galactic nuclei, and they characterize this distribution with an inner and outer power-law profile (Lauer et al. 2007).Mathematically, the Nuker law has the following form: , with γ and β representing the slopes of inner and outer power laws, respectively.The transition between these two regimes occurs at a given break radius, r b , and the sharpness of this transition is described by the parameter α.For NGC 4786, our GALFIT optimization converged on the following values: an ellipticity of ϵ = 0.20, α = 1.94, β = 1.56, γ = 0.00, and r b = 0. ′′ 46, which are typical values for cored-elliptical galaxies (Lauer et al. 2007).These cores are hypothesized to originate through scouring by massive supermassive BH binaries (Ravindranath et al. 2002;Thomas et al. 2014).The values for our NGC 5193 Nuker model were: an ellipticity of ϵ = 0.22, α = 6.28, β = 1.40, γ = 0.70, and r b = 1.′′ 39, which are values characteristic of power law galaxies (Lauer et al. 2007).
The next step involves estimating how much extinction to the observed H-band stellar light there is at the center of the galaxy.We follow the approach described by Kabasares et al. (2022), which uses the observed J − H color map of each galaxy to determine an estimate of A H , the extinction of the H-band stellar light originating behind the disk.First, we determined a median J − H color of 0.81 mag and 0.82 mag outside the dust disks of NGC 4786 and NGC 5193, respectively, and we determined the color excess, ∆(J − H) = (J − H) − (J − H) median as a function of position along the disk's major axis, averaging over a width of 4 pixels.We note that the J − H map of NGC 5193 is indicative of a central hole in the dust distribution and is supported by the CO(2−1) moment 0 map which displays a deficit of central CO emission as well.
To establish a relationship between extinction and color excess, we used Equations ( 1) and (2) from Boizelle et al. (2019) to generate a curve of ∆(J − H) as a function of V -band extinction, A V (see Figure 4 of Kabasares et al. 2022).This assumes the Viaene et al. ( 2017) embedded-screen model, which effectively models the circumnuclear dust disk as a thin, inclined disk that bisects the galaxy.Along a given LOS, the fraction of light that originates in front of the disk (f ) is unaffected by dust, while the fraction behind it (b) is obscured by screen extinction.The ratio of observed to intrinsic integrated H-band stellar light is represented mathematically as ].This is from Equation (1) in Boizelle et al. (2019), and assumes an intrinsically thin disk, where fractional disk thickness w = 0. Along the major axis of each disk, the fractions of light originating in front of and behind of the dust disk are assumed to be equal (f = b = 0.50).
The next key step in this process is converting the observed ∆(J − H) as a function of position along the disk's major axis into values of A H . Using our curve of ∆(J − H) versus A V , we can associate a unique value of A V (as well as A H ) to an observed ∆(J − H) value.As seen in Figure 4 of Kabasares et al. (2022), this is only valid up to a given turnover point.This is due to the fact that at large (A V > 5 mag) optical depths, variations in color begin to rapidly diminish, and so the same value of ∆(J − H) can correspond to both low and high A V values.Following the procedure outlined by Kabasares et al. (2022), we assumed the lower A V value, as the higher value implies that effectively all of the light originating behind the disk is lost due to extinction.We fit the color excess curve with a third-order polynomial up to the turnover point.To generate predictions of A V as a function of ∆(J − H), we use this polynomial's inverse.Then, we found the lower A V values corresponding to the observed ∆(J − H) along the disk's major axis.Finally, we set A H = 0.175A V based on the standard interstellar extinction law described in Rieke & Lebofsky (1985), which gives us a unique A H value for each observed color excess along the major axis.As stated earlier, this A H value applies only to light originating behind the disk.
Our simple dust correction implies A H = 0.22 mag and A H = 0.18 mag at the centers of the circumnculear disks in NGC 4786 and NGC 5193, which correspond to a reduction of approximately 20% and 15% of the stellar light originating behind the disk.We note that a proper treatment of determining the intrinsic stellar light distribution in the presence of a dusty circumnuclear disk likely requires the usage of radiative transfer models (De Geyter et al. 2013;Camps & Baes 2015) that account for the combined effects of extinction, light scattering, and the geometry of the dust disk itself.Even still, our simple method gives us a relatively straightforward way of producing an estimate of the assumed extinction, and consequently, the impact it has on the measured value of M BH .
With these estimated values of A H for each galaxy, we proceeded to mask the entirety of the dust disk in the H-band images except for the central nine pixels.This was done to anchor the model fit to the observed values at the center.The fluxes of these nine pixels are subsequently boosted by a factor of (0.50 + 0.50 × 10 −A H /2.5 ) −1 .We also tested this procedure with the central four pixels as well and found no significant difference between the two cases.With the entire dust disk masked out except for the pixels that have had their flux values increased, we re-fit the central 10 ′′ ×10 ′′ region again with a new Nuker model, but fixed the values of α, β, and r b to their values from the previous Nuker model to retain the larger scale properties outside of the dusty region.We then adjusted the inner slope parameter γ to find an optimal value where the central pixels of the Nuker model are nearly equal to the scaled flux values of the H-band image.For NGC 4786, this value is γ = 0.11, and for NGC 5193 it is γ = 0.75.
With this new Nuker model, the final steps in our dust correction process are to replace the pixels within the dust disk region in the H-band image with the corresponding pixels in the Nuker model, and to fit this dust-corrected H-band image with a new MGE.As will be discussed in Section 5, these correspond to the MGEs that are used in our fiducial dynamical models.The dust-corrected MGE components are displayed in Table 1.

DYNAMICAL MODELING
The BH masses are measured by modeling the observed gas kinematics in the ALMA data cubes with a thin, rotating disk model.Our models use a minimum of nine free parameters which include BH mass M BH , H-band mass-to-light ratio Υ H , disk inclination angle i, disk position angle Γ, disk dynamical center (x c , y c ), CO flux normalization factor F 0 , and turbulent velocity dispersion σ 0 .The details of our modeling process are described by Kabasares et al. (2022), but we summarize the key aspects here.
We create synthetic data cubes that model the observed CO line profiles and fit them directly to the ALMA data cubes.To start, we build a model circu-lar velocity field on a grid that is oversampled by a factor of s = 3 relative to the size of the spatial pixels in the data cube.The circular velocity at each grid point is a function of radius and is determined by the quadrature sum of the velocity contributions from the BH (treated as a point source with mass M BH ), the host galaxy's stars, and from the gas disk itself.We neglect the contribution of dark matter as the enclosed mass profiles on scales comparable to the circumnuclear disks are thought to be dominated by the stars.To derive the circular velocity profile due to the stellar mass, we deprojected the MGEs of each galaxy described in Section 3 using routines from the JamPy package (Cappellari 2008) under the assumptions that NGC 4786 and NGC 5193 are oblate, axisymmetric, and have inclination angles of 70 • and 60 • , respectively, based on initial gas-dynamical modeling fits.Each of these deprojections represents a possible three-dimensional luminous density for the galaxy.With each model iteration, the stellar mass circular velocity profiles are scaled by √ Υ H , the square root of the H-band mass-to-light ratio.
The circular velocity contribution of the gas disk is obtained by numerically integrating projected gas mass surface densities that are calculated on the same annular regions described in the calculation of M gas in Section 2.1.1.We assumed the gas was distributed in a thin disk and determined the midplane circular velocity of a test particle orbiting in the disk's gravitational potential using Equation 2.157 in Binney & Tremaine (2008).
For a given disk inclination angle i, position angle Γ, and assumed (fixed) distance to the galaxy, D, we calculate the LOS projection of the model velocities as seen on the plane of the sky.At each point in the disk, we model the emergent CO line profile as a Gaussian and use the projected LOS velocity and an assumed spatially uniform turbulent velocity dispersion term, σ 0 , to build the model line profile.Given that the spectral axis of the ALMA data cubes is in frequency units, we transform the projected LOS velocity and velocity dispersion into a frequency line centroid and line width using the redshift, z obs (related to the free parameter of the systemic velocity of the galaxy through v sys /c = z obs ).
Once the model data cubes are constructed, the Gaussian line profiles must be weighted by a CO flux map and each image slice needs to be convolved with the ALMA synthesized beam.The flux map is created by using a 3D mask generated by the 3DBarolo program (Di Teodoro & Fraternali 2015), multiplying this mask with the ALMA data cube, and summing this product along the spectral axis to form a 2D image.To deconvolve this image, we use an elliptical Gaussian PSF with major and minor FWHMs that match the speci-   The final step in the modeling process is to minimize χ 2 between our model and the ALMA data.We compute χ 2 on elliptical spatial regions centered on the disk centers shown in Figures 5 and 6.Given the small number of resolution elements across each disk, we opted to fit our models to nearly the full extent of each disk, though we explore the systematic impact of this choice in Section 5.3.Our elliptical fitting regions have an axial ratio of q = 0.50 and semimajor axis lengths of a = 0. ′′ 55 and a = 1.′′ 0 for NGC 4786 and NGC 5193, respectively.Additionally, to mitigate the impact of neighbor-ing pixel-to-pixel correlation, we further block-averaged the cubes by a factor of 4. On these block-averaged scales, we minimized 2 , where d i , m i , and σ i represent the flux density of the 3D data, model, and noise cubes at a given pixel.The σ i values are determined from a 3D noise model that was generated for each galaxy described in detail in Section 4.2.1 of Kabasares et al. (2022), but we briefly describe its construction here.We first calculated the noise in each frequency channel as the standard deviation of emissionfree pixels on the block-averaged scale.This is done in a version of the data cube prior to primary beam correction as the noise is spatially uniform at this stage.Then, we took the primary beam cube generated during data calibration and reduction and block-averaged it by the same factor.The 3D noise model is then generated by dividing the uniform rms noise value by the primary beam for each frequency channel.We fit our dynamical models to 31 frequency channels in the NGC 4786 data cube that correspond to velocities of |v LOS − v sys | ≤ ∼300 km s −1 .This velocity range slightly extends past the channels with visible CO emission.After block-averaging our model and data cubes, the elliptical spatial region contains 12 pix-els per channel, or a total of 372 data points in the entire model fit.For the NGC 5193 data cube, we fit our models to 55 spectral channels corresponding to |v LOS − v sys | ≤ ∼280 km s −1 .On the block-averaged scale, the elliptical spatial region contains 58 data points per channel which gives an overall total of 3190 data points used in the χ 2 minimization.

NGC 4786 Dynamical Modeling Results
We present the results for three dynamical models (A,B,C) for NGC 4786 in Table 2.The difference among them is the input host galaxy MGE model, which we described in Section 3. In summary, dynamical models A-C yield a range in BH mass that spans (3.9−5.8)×10 8 M ⊙ with reduced χ 2 (χ 2 ν ) values between 1.421 and 1.488 over 363 degrees of freedom.Using the dust-corrected MGE as an input, dynamical model C is the statistically best-fitting model, and we adopt it as our fiducial model to use in our systematic tests of the error budget.With χ 2 ν = 1.421, model C is not a formally acceptable fit when considering the degrees of freedom, which should achieve χ 2 ν ≤ 1.125 when using a significance level of 0.05.The χ 2 values may be high in part due to inherent limitations of deconvolving the CO flux map.However, while changes to the intrinsic flux map may improve χ 2 , they are unlikely to significantly affect M BH (e.g., Marconi et al. 2006;Walsh et al. 2013;Boizelle et al. 2021).An in-depth analysis of the systematic uncertainties of the measurement is described in Section 5.3.We present the major axis PVD, moment maps, and CO line profiles extracted from the data and fiducial model cubes in Figures 4, 5, and 7.The comparisons highlight good overall agreement between the data and model.The model PVD is able to emulate features such as the broad distribution in rotational velocity that is observed within r ≤ 0. ′′ 5 on both the approaching and receding sides of the disk, as well as the decrease in velocity towards the center.An analysis of the moment maps shows that the model velocity field captures most of the behavior seen in the outer portions of the disk, where discrepancies in LOS velocity are < ∼20 km s −1 , although noticeable disagreement is seen at the center, particularly along the minor axis.As described by Barth et al. (2016b), along an inclined disk's minor axis, the projected distance between the nucleus of a galaxy and a point along the minor axis is compressed by a factor of cos(i), and poor spatial resolution across the minor axis can lead to pronounced beam smearing.Given the large ALMA beam (FWHM = 0. ′′ 31 = 93 pc) relative to this small disk, it is unsurprising that beam-smearing effects are most evident in this region.The moment 0 map also highlights discrepancies between data and model in the central region, with the moment 0 map indicating a higher model CO surface brightness than what is observed.The CO line profiles can display complex characteristics, but even though the fine details within the broad and asymmetric data line profiles may be missed, our model CO line profiles generally capture their overall shape.
As for other free parameters in the model, our dynamically determined Υ H values, especially for model A, are higher than typical Υ H values (Υ H ≤ 1.30 M ⊙ /L ⊙ ) seen in single stellar population models that assume either a Kroupa (2001) or a Chabrier (2003) initial-mass function (IMF) for an old (10-14 Gyr) stellar population with solar metallicity, but our measurement is consistent within systematic uncertainties with models assuming a Salpeter (1955) IMF (Υ H ∼ 1.51 − 1.84 M ⊙ /L ⊙ ).This is expected, given that galaxies with large σ ⋆ tend to follow Salpeter-like or even heavier IMFs (e.g., Cappellari et al. 2013;Smith 2020).This is also consistent with Mehrgan et al. (2024), who find that the centers of massive ETGs typically have Salpeter-like or heavier mass-to-light ratios.Mehrgan et al. ( 2024) also find Υ gradients in the majority of massive ETGs, which can affect the measured M BH .However, the Υ gradients are measured at radii of ∼ 0.2−1.0kpc, beyond the physical scale of the dust disks of the galaxies in this study.
The σ 0 parameter remains fairly low between 9.3 − 10.8 km s −1 , which is less than the data cube's channel spacing and is consistent with other gas-dynamical modeling of ALMA data (Barth et al. 2016a;Boizelle et al. 2019Boizelle et al. , 2021;;Cohn et al. 2021;Kabasares et al. 2022).The flux normalization factor F 0 is found to be between 1.54 − 1.56 in our models and is higher than values seen in previous works, where it is typically closer to unity.Upon inspecting Figure 5, a comparison of the data and best-fit model's moment 0 map reveals a noticeable difference in surface brightness, particularly close to the disk's center, where the data appears to have faint CO emission.We explore the systematic effect of the input flux map on the mass measurement's error budget in Section 5.3.
Based on previous dynamical modeling work, we expect the statistical uncertainty on a given dynamical model fit to be much smaller than the uncertainty associated with the extinction corrections in our host galaxy models.To determine the statistical uncertainty for the NGC 4786 BH mass measurement, we used a Monte Carlo resampling procedure.We generated 100 noiseadded realizations of our fiducial model by adding noise to each pixel of the model cube.At each pixel, we drew a randomly sampled value from a Gaussian distribution with a mean of zero and a standard deviation equal to the value of our 3D noise cube at the same pixel.We re-optimized a dynamical model to each of our 100 altered realizations, using the values in Table 2 as our initial guesses.We list the standard deviation of each free parameter as the 1σ uncertainty under the results for model C.For the BH mass, the distribution has a mean of M BH = 5.0 × 10 8 M ⊙ and a standard deviation of 0.2 × 10 8 M ⊙ or approximately 4% of the mean.

NGC 5193 Dynamical Modeling Results
We present three dynamical models (D,E,F) for the NGC 5193 data cube in Table 2.As in the case of NGC 4786, these dynamical models used three unique host galaxy MGEs that treat the effects of the dust on the stellar light differently.Dynamical models D-F yield a range of M BH = (1.4− 2.9) × 10 8 M ⊙ , Υ H = 1.46 − 1.69 M ⊙ /L ⊙ , σ 0 = 3.1 − 6.7 km s −1 , and F 0 = 1.14 − 1.16 with χ 2 ν = 2.096 − 2.541 over 3181 degrees of freedom.While our models are generally successful at reproducing the observed kinematics over a majority of the disk, a formally acceptable fit would achieve χ 2 ν ≤ 1.042 assuming a level of significance of 0.05.The dynamically determined Υ H values are consistent with single stellar population models assuming a Salpeter (1955) IMF.As with NGC 4786, the σ 0 values are small and are less than the channel spacing in the data cube.The range in F 0 are closer to unity than what was found for NGC 4786, though a lack of central CO emission is noticeable in the data visualizations.The major axis PVD, moment maps, and CO line profiles that compare the data and the best-fit model (model F), are shown in Figures 4, 6, and 7.The PVDs and moment maps show that model F emulates the observed PVD structure and disk kinematics over nearly the full extent of the disk.Similarly to NGC 4786, the most noticeable differences are in the central parts of the disk, as the data appears to be more CO-faint than our model.Additionally, the observed line profiles seen in Figure 7 extracted near the center of the disk show complex and asymmetric structure that our model cannot fully describe, as channel-to-channel variations in the amplitudes of the data line profiles are not entirely reproduced, although model F does generally capture their overall shapes well.
We performed a Monte Carlo simulation to determine the statistical uncertainty of the measurement.As we had done for NGC 4786, we created 100 realizations of model F by adding Gaussian noise.Then, we optimized each realization to produce a new estimate of M BH and create a distribution.The distribution of M BH was centered around a mean of 1.4×10 8 M ⊙ , and had a standard deviation of 0.03 × 10 8 M ⊙ , or 2% of the mean.We list the standard deviation of all other free parameters determined by this simulation under the results of model F in Table 2.

Error Budget of the Mass Measurements
While the statistical uncertainties on M BH derived from our Monte Carlo procedure are small, there are other sources of uncertainty that stem from different aspects of our model construction.It has been shown by previous dynamical-modeling studies (Boizelle et al. 2019(Boizelle et al. , 2021;;Cohn et al. 2021;Kabasares et al. 2022) that the statistical model-fitting uncertainties for a given dynamical model vastly underestimate the total error budget of a given measurement when considering the uncertainties from model systematics.To assess the impact of the systematic uncertainties on the error budget in each galaxy, we took our statistically best-fit dynamical models and performed a number of systematic tests that involved changing aspects of the model construction.We briefly describe and list these changes below: • Dust extinction: We explored the impact of dust extinction by creating our three MGE host galaxy models (unmasked, dust-masked, and dust-corrected) described in Sections 3.1 and 3.2.Dust is clearly present at the centers of both systems, and previous gas-dynamical studies (Boizelle et al. 2019(Boizelle et al. , 2021;;Cohn et al. 2021;Kabasares et al. 2022) have also shown that even in the best cases where the BH SOI is well-resolved, differences in the assumed host galaxy profiles can lead to large discrepancies in M BH .Therefore, while the intrinsic host galaxy mass profile and the uncertainty in its inner slope may be difficult to ascertain, by building a set of MGE models that account for the presence of dust in different ways, we can effectively produce a set of models that bracket the likely range of profiles.
We adopt model C for NGC 4786 and model F for NGC 5193 as the fiducial model for each galaxy and perform the remaining systematic tests on them for a number of reasons.First, these two models are the statistically best-fitting dynamical models for NGC 4786 and NGC 5193.In addition, previous ALMA BH mass measurements (Boizelle et al. 2019(Boizelle et al. , 2021;;Cohn et al. 2021;Kabasares et al. 2022) have shown that the statistical model fitting uncertainties for a given dynamical model vastly underestimate the total error budget of a given measurement when considering uncertainties associated with the host galaxy extinction corrections, and the dust-corrected MGE models used in models C and F are our best estimate of the underlying host galaxy profile.
• Radial motion: Following the approach described by Kabasares et al. (2022), we allow for radial motion in our dynamical models, by incorporating a simple toy model which is controlled by a param-eter α.This parameter controls the balance between purely rotational (α = 1) and purely radial inflow (α = 0) motion.
• Turbulent velocity dispersion: The gas velocity's dispersion is changed from a spatially uniform term of σ(r) = σ 0 to a Gaussian that is a function of radius: σ(r) = σ 0 +σ 1 exp[−(r−r 0 ) 2 /2µ 2 ].This model adds three additional free parameters, with σ 1 , r 0 , and µ representing the Gaussian's velocity amplitude, radial offset from r = 0, and standard deviation, respectively.While this model is not Figure 6.Moment maps for NGC 5193 constructed from the ALMA CO(2−1) data cube (left) and its fiducial model (center, model F; see Section 5.2).Shown are maps of moments 0, 1, and 2, corresponding to surface brightness, LOS velocity vLOS, and LOS velocity dispersion σLOS.The units for the surface brightness map are mJy km s −1 pixel −1 , and the units for the vLOS and σLOS maps are km s −1 .The systemic velocity of 3705 km s −1 estimated from our dynamical models has been removed from vLOS.Maps of (data-model) residuals are shown in the rightmost column.The coordinate system is oriented such that +x corresponds to East and +y corresponds to North.While the line profile fits have been determined at each pixel of the full disk, the elliptical fitting region used in calculating χ 2 is denoted in the top left panel with a yellow ellipse.The synthesized beam is represented by an open ellipse in the bottom left corner of the same image.
motivated by any physical mechanism, it allows for more variation in the form of the velocity dispersion, and previous modeling has showed some preference for a moderate increase in σ toward the center (e.g., Barth et al. 2016b).
• Fit region: We adjust the elliptical spatial region used to optimize our dynamical models described in Section 4. The new fitting regions for NGC 4786 and NGC 5193 are ellipses with semimajor axes of approximately 1.25 ALMA resolution elements centered on the disk's dynamical center.With this setup, the models are fit to parts of the disk that are more sensitive to the BH's gravitational potential as opposed to the host galaxy's, but there is now a larger fraction of pixels that are in the central beam-smeared region of the disk.
• Gas mass: The gas disk's mass contribution to the gravitational potential is removed, and model fits incorporate only the gravitational potential of the BH and the host galaxy.This is done by setting the gas disk's contribution to the total circular velocity to 0 km s −1 .• Block-averaging factor : The process of blockaveraging our dynamical models leads to coarser angular resolution, and could potentially limit the models' ability to constrain M BH .To test for this possibility, we optimized a dynamical model where no block-averaging was performed.
• Oversampling factor : We originally built our models on a grid that was oversampled by a factor of 3 relative to the ALMA data.We set this factor to 1 to test the effect of building our models on a grid that is equal in size to the original ALMA spatial scale.
• Input flux map: We built a different input flux map to weight the CO line profiles.This flux map was constructed by fitting an azimuthally normalized 3D tilted-ring model (Rogstad et al. 1974;Begeman 1989) in the 3DBarolo program to the ALMA data cubes (Di Teodoro & Fraternali 2015).
The program returns a 3D model data cube from which we can produce a flux map in the same manner as a regular ALMA data cube described in Section 4 and use as an input in our dynamical models.
The largest shifts in M BH are seen in the systematic tests that involve changing the host galaxy model to account for extinction.We describe the changes observed from these tests below and list the results of the other systematic tests that do not involve changes to the host galaxy MGE model in Table 3.
The value of M BH in both galaxies is highly sensitive to how we account for dust extinction, which is highlighted by the differences in the unmasked, dust-masked, and dust-corrected MGE models shown in Table 2.For NGC 4786, the value of M BH decreases by about 22% to Note-Results of the systematic tests performed on the fiducial dynamical model for each galaxy, not including tests that involved changes to the host galaxy MGEs.We list the new value of MBH that the optimization converged to, as well as the percent change from the BH masses determined for the fiducial models C and F.
M BH = 3.9 × 10 8 M ⊙ when using the unmasked MGE and increases by 16% to M BH = 5.8×10 8 M ⊙ when using the dust-masked MGE.In the case of NGC 5193, the dynamical model using the unmasked MGE converged to a BH mass of M BH = 1.5 × 10 8 M ⊙ , only 7% higher than when using the fiducial model's dust-corrected MGE, but a large factor of 2 increase to M BH = 2.9 × 10 8 M ⊙ is observed when using the dust-masked MGE.As stated earlier, the dust-corrected MGE models are our best estimate of the intrinsic host galaxy profile.As a final systematic test, we fixed M BH = 0 M ⊙ in our dynamical models and reoptimized with the dustcorrected MGE models used in models C and F to see if the models could still reliably emulate the observed gas kinematics without the need for a BH.Again, this resulted in poorer fits to the data, with χ 2 ν = 1.974 in NGC 4786 as well as a large increase in Υ H to 2.60 M ⊙ /L ⊙ to compensate for the lack of a supermassive BH.As mentioned earlier, this is a higher Υ H value than what is seen in single stellar population models.We see a similar result in the dynamical model for NGC 5193, with the model fit yielding a much higher χ 2 ν = 2.978 and a higher Υ H = 1.58 M ⊙ /L ⊙ .The results show that M BH = 0 M ⊙ is highly disfavored for both galaxies.
As shown in Table 3, most of the systematic tests that did not involve adjustments to the MGE models led to relatively insignificant changes to M BH , but there were a few exceptions.For NGC 4786, the most profound (> 10%) changes were a 15% decrease due to the change in oversampling factor and a 21% increase due to the choice of flux map.The models are constructed on an oversampled grid in order to account for potentially steep velocity gradients in the disk, as pixels near the disk center can contain molecular gas spanning a large velocity range.Without any oversampling, these velocity gradients can be be missed, and can subsequently lead to models converging on a less massive BH.As for using the new flux map, our dynamical models not only converged on a higher BH mass, but also substantially different values of F 0 = 0.93 and i = 58.6 • with an improved χ 2 ν = 1.354.The change of over 10 • in inclination angle, as well as a 40% decrease in the flux normalization factor indicate significant differences when modeling the disk's structure.As stated earlier, this flux map was created from a 3D model data cube generated from an automated 3D tilted-ring fit in 3DBarolo.The tiltedring model allows Γ and i to be different for each ring, and converged on ring inclinations of approximately 60 • , a significant difference from our flat disk models.The differences in the empirically measured and tilted-ring model flux maps can be attributed to the assumptions in the 3DBarolo model fit as well as variations in the CO surface brightness across the disk, especially near the disk's center, where the disk is CO-faint.These features are not encapsulated in the azimuthally normalized tilted-ring model, which has a relatively constant surface brightness across the disk.This leads to large differences in the overall flux normalization factor and the inferred disk structure.We attribute the large (13%) change in M BH in NGC 5193 when changing its dynamical model's flux map to these reasons as well.
These tests highlight a clear degeneracy between the BH and stellar mass components in our dynamical models.The systematic uncertainties from the host galaxy modeling dominate over the statistical (≈2%) and the distance uncertainty of ≈15% in NGC 4786 and ≈7% in NGC 5193, as well as the other systematic uncertainties (≈5 − 20%) not associated with the host galaxy models.While the mass range for M BH found in Table 2 is at the ≈20% level in NGC 4786 and the factor of 2 level in NGC 5193, our dynamical models prefer the presence of a central supermassive BH to reproduce the observed gas kinematics, as opposed to models where no BH is present.
When considering only the fiducial host galaxy MGE model, and including the uncertainties from the systematic tests in Table 3, the distance to the galaxy, and the statistical fluctuations in the data, the range in M BH is (M BH /10 8 M ⊙ ) = 5.0 ± 0.
for NGC 5193 if we add the negative and positive systematic uncertainties from Table 3 in quadrature.However, these ranges in M BH neglect contributions from the uncertainties in the extinction correction of the host galaxy models.
We will incorporate these systematic uncertainties for a number of reasons.First, while our dust-corrected MGE model is our best estimate of the underlying host galaxy model, we must account for the potentially large variation in the inner slope of the stellar mass profile in each galaxy.While dust attenuation on the host galaxy light may not appear obvious when comparing the data and model isophotes or surface brightness profiles in Figures 2 and 3, the resulting deprojected host galaxy models produce a broad range in M BH .
In other works that feature M BH measurements in ETGs (Davis et al. 2017(Davis et al. , 2018;;Smith et al. 2021;Ruffa et al. 2023), the surface brightness of the host galaxy has typically been parameterized with only dust-masked MGEs.While masking out the most dust-obscured features in the image prior to fitting an MGE may yield better models than without any masking, it does not fully address the problem of extinction, and limits the model fit to the remaining pixels that are not completely unaffected by dust.As shown in this work and in other previous ALMA M BH measurements (Boizelle et al. 2019(Boizelle et al. , 2021;;Cohn et al. 2021;Kabasares et al. 2022), uncertainties from the extinction correction far exceed the formal statistical uncertainties, so accounting for a range in the inner stellar mass profile slope and its effect on the inferred M BH should be explored.Therefore, we strongly advocate that future M BH studies with ALMA explore this often overlooked component of an M BH measurement's error budget.
After adding the systematic uncertainties associated with the extinction correction in quadrature with the uncertainties associated with the systematic tests in Table 3, the final measured range in . We compare our measured ranges of M BH with predicted values of M BH from BH-host galaxy scaling relations in the following section.

DISCUSSION
Our ALMA gas-dynamical mass measurements are the first attempts to measure the central BH masses in NGC 4786 and NGC 5193.In the following subsections, we determine a range for the projected radius of the BH SOI in each galaxy, and describe how using ALMA observations that do not fully resolve this scale leads to large systematic uncertainties in the mass measurement.In addition, we compare our measured ranges of M BH in each galaxy to predicted ranges from the scaling relations of Kormendy & Ho (2013).
6.1.The Impact of Resolving the BH Sphere of Influence in ALMA Gas-dynamical Measurements Our results indicate that the BH's projected radius of influence, r SOI , is not well resolved in either galaxy, which leads to a degeneracy between stellar and BH mass in our models.We calculated r SOI in NGC 4786 and NGC 5193 by determining the radius where M ⋆ = M BH for each of the three input MGE models and highlight the range of r SOI in each galaxy in Figure 8.For NGC 4786, r SOI is between 73 pc (0. ′′ 24) and 90 pc (0. ′′ 30), whereas the ALMA beam FWHM is 93 pc (0. ′′ 31), indicating that the SOI is nearly resolved along the major axis of the disk.As pointed out by Barth et al. (2016b), however, the threshold for a successful BH mass measurement rises considerably at higher inclination angles, as projected distances along the minor axis of an inclined disk become compressed by a factor of cos(i).At an inclination angle of 70 • , r SOI in NGC 4786 is unresolved along the disk's projected minor axis by a factor of nearly 3.As for NGC 5193, we find that r SOI is between 11 pc (0. ′′ 05) and 26 pc (0. ′′ 12), and is well below the resolution limit of ∼69 pc (0. ′′ 31).
Following the work of Rusli et al. (2013), we compare our measured r SOI values to the average ALMA beam size through ξ = 2r SOI /θ FWHM .For NGC 4786, we find that ξ = (1.6 − 1.9) and for NGC 5193, we find ξ = (0.3 − 0.8).Davis (2014) suggests that observations that satisfy ξ ∼ 2 are adequate for the purpose of conducting molecular gas-dynamical BH mass measurements.This figure of merit was designed to aid in the planning of observational campaigns focused on estimating central BH masses, and is a less stringent threshold than traditional considerations.However, as our results show, measurements based on observations in this regime and lower (ξ ∼ 1) will have increasingly larger uncertainties.
While the ALMA observations for NGC 4786 approximately satisfy the aforementioned figure of merit, the sensitivity of M BH to the choice of host galaxy model in NGC 4786 is still readily apparent.Detecting the en-  2. Mgas(r) is calculated via: rv 2 gas /G.The red dotted line indicates the angular resolution of the ALMA observations, whereas the black dotted line represents the outer edge of the dust disk as measured along the major axis.The gray shaded region indicates the range of MBH determined by our dynamical models using the different input MGE models, while the yellow shaded region enclosed within the black dotted lines indicate the range of the radius of the black hole sphere of influence, rSOI.We have defined rSOI to be the radius where the stellar and BH masses of a given dynamical model are equal.hancement of a tracer particle's circular velocity due to the presence of a supermassive BH is more difficult when the slope of the stellar mass profile in the central region of a galaxy is not well constrained.As we see in Figure 8, differences in the overall shape of the different host galaxy mass profiles are minimal at the disk edge and beyond, as expected.However, at radii less than the projected ALMA resolution limit, there are noticeable differences in the overall shape of the stellar mass profiles.The process of converting the observed brightness distributions into mass profiles is complicated by the presence of dust, which limits a dynamical model's ability to reliably separate the BH and galaxy contributions to the overall gravitational potential.Figure 8 reveals an inner bump in the M ⋆ (r) profile for the unmasked MGE.The bump is unphysical and is a limitation of MGEs that is exacerbated by both the annular extinction of the dust and the constraint that the enclosed mass of each MGE should match at the disk edge.This aspect of MGEs highlights the need for dust-masked and dust-corrected MGEs when modeling the host galaxy's light in cases where dust is readily apparent.
In the case of NGC 5193, we see nearly a factor of 2 discrepancy in the measured BH mass from our dynamical models.Quantifying the uncertainty in the inner slope of the NGC 5193 host galaxy mass profile is challenging with the central dust disk, and because the projected BH SOI is unresolved, the range in our measured M BH is broad.With the BH SOI unresolved by a factor of 2 to 3, the BH mass is a smaller fraction of the total enclosed mass on scales comparable to the observation's resolution limit, and thus, our measurements are even more sensitive to the differences among the stellar mass models.Similarly to NGC 4786, our stellar mass models for NGC 5193 shown in Figure 8 are well-matched at the edge of the disk, but are noticeably different within the central ALMA resolution element.Moreover, as in the case of NGC 4786, the unmasked MGE displays an unphysical bump at small radii in the M ⋆ (r) profile, which we attribute to the aforementioned reasons provided for NGC 4786.
An additional factor in the measurement uncertainty is the apparent deficit of CO emission in the central region of both galaxies.This deficit is apparent in the moment 0 maps of both NGC 4786 and NGC 5193 shown in Figures 5 and 6, as residuals between the data and best-fit model highlight prominent discrepancies at the center of each disk.Furthermore, the J −H map of NGC 5193 shows additional evidence of a central hole in the circumnuclear disk.In previous work, such as in the case of NGC 6861, these holes can limit the measurement of M BH to using kinematic information beyond r SOI and prevent tight constraints on the M BH (Kabasares et al. 2022).

Comparison to Predictions from BH Scaling Relations
The ALMA observations were designed to probe scales comparable to r SOI using BH mass estimates from the M BH − σ ⋆ relation of Kormendy & Ho (2013) using the definition: r SOI = GM BH /σ 2 ⋆ .For NGC 4786, the relation predicts a BH mass of M BH = 1.5×10 9 M ⊙ , whereas for NGC 5193, it predicts 3.4 × 10 8 M ⊙ when using the σ ⋆ values of 286 km s −1 and 205 km s −1 from Hyperleda (Makarov et al. 2014).The measured M BH values for NGC 4786 and NGC 5193 fall below the predicted value, but given the intrinsic scatter of 0.28 dex in the M BH −σ ⋆ relation, the higher end of our M BH range for NGC 5193 falls within the predicted scatter.Thomas et al. (2016) suggests that for core galaxies, the core radius r b is a more robust proxy for M BH than σ ⋆ .Given that NGC 4786 exhibits properties of a cored-elliptical galaxy, we input our GALFIT value of r b = 0. ′′ 46 = 138.5 pc into the M BH − r b (Nuker) relation of Thomas et al. (2016).We find that the expected M BH is 9.7 × 10 8 M ⊙ .Our mass range of NGC 4786 falls below this relation's predicted value, but within the scatter of the relation.
To compare our mass ranges with the M BH − L bulge,K relation of Kormendy & Ho (2013), we used our dustcorrected MGE models in Table 1 and assumed a color of H − K = 0.2 mag based on stellar population models of an old stellar population with solar metallicity (Vazdekis et al. 2010), as well as absolute solar magnitudes of M ⊙,H = 3.37 mag and M ⊙,K = 3.27 mag (Willmer 2018).By adding the total luminosity of each component in our H-band MGEs and converting them to the K-band, we found that L bulge,K = 3.5 × 10 11 L ⊙ for NGC 4786 and L bulge,K = 1.3 × 10 11 L ⊙ for NGC 5193.These yield predicted values of M BH of 2.5 × 10 9 M ⊙ and 7.4 × 10 8 M ⊙ for NGC 4786 and NGC 5193, respectively.Given the scatter of 0.30 dex, our measured BH mass range for both galaxies falls below their respective predicted value and the scatter of the relation.
Finally, we compared our measured ranges with the M BH − M bulge relation by converting the total H-band luminosities of our dust-corrected MGEs and multiplying them by our Υ H values for these models listed in Table 2.This was done under the assumption that there are no gradients in Υ H .We found M bulge = 5.7×10 11 M ⊙ for NGC 4786 and M bulge = 1.7×10 11 M ⊙ for NGC 5193.As discussed in Zhu et al. (2021), a nonneglible fraction of the mass in elliptical galaxies may not belong to the "bulge" component of the galaxy, and thus, a more nuanced mass decomposition of the host galaxy is needed to reliably calibrate empirical correlations between M BH and bulge properties.Therefore, our calculated masses should be taken as approximations that likely overestimate M bulge .Indeed, our estimated M bulge values correspond to predicted values of M BH = 3.8 × 10 9 M ⊙ and 9.0 × 10 8 M ⊙ in NGC 4786 and NGC 5193.Our measured values are below these predictions and the intrinsic 0.28 dex scatter in the relation.

CONCLUSION
We present the first dynamical measurements of the central BH masses in ETGs NGC 4786 and NGC 5193 using ALMA CO(2−1) observations that each have 0. ′′ 31 resolution.In both galaxies, a circumnculear disk in orderly rotation with LOS velocities of ∼270 km s −1 relative to the systemic velocity is observed.
For NGC 4786, our dynamical models constrain the central BH mass to be (M BH /10 8 M ⊙ ) = 5.0 ± 0.2 [1σ statistical] +1.4 −1.3 [systematic] ± 0.75 [distance] by fitting the ALMA data with three different host galaxy MGE models.Upon conducting numerous systematic tests, we found that the systematic uncertainties associated with the extinction correction of the host galaxy MGE models were the dominant contributor to the overall error budget.
In the case of NGC 5193, we fit dynamical models to the ALMA CO(2−1) data with three host galaxy MGE models as well.With ALMA observations that resolve on scales of a few BH SOI radii, we found that the uncertainty in the BH mass is around the factor of 2 level.Our measured range of (M BH /10 8 M ⊙ ) = 1.4 ± 0.03 [1σ statistical] +1.5 −0.1 [systematic] ± 0.1 [distance] for NGC 5193 highlights the importance of accounting for differences in the host galaxy models due to the presence of dust, as the BH comprises a smaller fraction of the total enclosed mass on scales comparable to the observation's resolution.
While the mass range for M BH is broad and our dynamical models do not tightly constrain M BH in either galaxy, models that contain a central supermassive BH fit the observed ALMA data much better than models without one.Additionally, our models underscore the importance of incorporating a range of plausible inner slopes in the host galaxy mass models when calculating the error budget of an M BH measurement.This incorporation is necessary when conducting M BH measurements in dusty systems with ALMA observations that do not fully resolve the projected BH SOI, as stellar mass and BH mass in a dynamical model can become degenerate.Future higher resolution and higher sensitivity observations with ALMA could potentially break the observed BH and stellar mass degeneracy in our dynamical models for NGC 4786 and NGC 5193 and lead to a tighter range on M BH .The higher sensitivity observations would aid in the detection of possible faint CO emission within r SOI .For now, our measurements place important mass constraints on the central BHs in two ETGs that did not have a prior dynamical M BH measurement, and highlight important limiting factors and considerations for future gas-dynamical M BH measurements with ALMA.& Hook 2004), JamPy (Cappellari 2008), scikit-image (Van der Walt et al. 2014)

APPENDIX
In Tables 4 and 5, we list the components of the unmasked and dust-masked MGEs for NGC 4786 and NGC 5193 described in Section 3.1.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.Primes indicate projected quantities.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.Primes indicate projected quantities.
ALMA Observations We obtained ALMA imaging of NGC 4786 and NGC 5193 from ALMA Programs 2015.1.00878.S and 2017.1.00301.S, respectively.NGC 4786 was observed on 23 July 2016 for approximately 21 minutes with a maximum baseline of 1110 m.The observation consisted of three spectral windows targeting continuum emission and one spectral window targeting the redshifted CO(2−1) emission line.The continuum windows had a channel resolution of 15.625 MHz and covered the following frequency ranges: 227.14 − 231.14 GHz, 239.47 − 243.47 GHz, and 241.78 − 245.78 GHz.The emission line spectral window had a channel resolution of 3.906 MHz and spanned the frequencies between 225.14 − 228.89GHz.The uv-plane visibilities were reduced and calibrated in version 4.5.3 of the Common Astronomy Software Applications package (CASA; Mc-

Figure 1 .
Figure 1.HST F110W (J-band), F160W (H-band), F110W − F160W (J − H), and ALMA CO(2−1) images of NGC 4786 (top) and NGC 5193 (bottom) showing the co-spatial alignment of the gas and dust.The ALMA integrated flux density maps were created by summing channels in the data cubes that displayed visible CO emission.Pixels with emission were identified with an automatically generated mask by the 3DBarolo program (Di Teodoro & Fraternali 2015).In the J − H maps, light regions correspond to redder colors and dark regions are bluer than the surrounding starlight.North is up and East is to the left in each image.Colorbars are shown in magnitude units.

Figure 2 .
Figure 2. 2D isophote maps comparing the observed HST WFC3 F160W isophotes to those of our three MGEs for NGC 4786 (top) and NGC 5193 (bottom).Black contours represent isophotes from the F160W images, while red contours are for the MGE models.For each image, the central ≈100 ′′ × 100 ′′ region is displayed with an inset of the innermost 3. ′′ 5 × 3. ′′ 5 region in the top right corner.The gray ellipse shown within each inset indicates the size and orientation of the circumnuclear dust disk.Arrows in the middle panels indicate the orientation of North and East for each galaxy.in Figure 3 also shows good agreement at intermediate radii, though discrepancies within the dusty region are noticeable.Despite slight mismatches in the model and data surface brightness profiles, the dust in NGC 4786 and NGC 5193 appears to have a less noticeable impact on the observed H-band surface brightness distribution in each galaxy in comparison to what has been seen in previous work, such as in the ETGs NGC 1380 and NGC 6861, where the H-band isophotes become non-elliptical and asymmetric within the dust disk (Kabasares et al. 2022).Additionally, we show the ellipticity as a function of semi-major axis for the H-band data and the various MGE models in the bottom panels of Figure 3.The ellipticity profiles of the unmasked and dustcorrected MGE models are in good agreement with their respective images.The unmasked H-band photometry is poorly behaved within the dust lane.The ellipticity

Figure 3 .
Figure 3. (Top) A comparison of the observed and modeled H-band surface brightness profiles of NGC 4786 (left) and NGC 5193 (right).The surface brightness measurements are made with the Python-implementation of the sectors photometry routine(Cappellari 2002) which performs photometry along evenly spaced sectors from the major axis to the minor axis and averages measurements over the four quadrants of the image.(Bottom) A comparison of the radial ellipticity of the observed H-band photometry and the 2D MGE models for NGC 4786 (left) and NGC 5193 (right).We use the photutils.isophoteroutine(Jedrzejewski 1987;Bradley et al. 2023) which fits elliptical isophotes to a galaxy image to determine the ellipticity of each isophote along the semi-major axis.For each panel, the red squares are the observed values from the H-band image, while blue squares are dust-corrected values described in Sections 3.2.The solid lines in each panel correspond to profiles extracted along the major axis for each of our 2D MGE models.Red lines are for dust-masked MGE models whereas black and blue lines represent dust-unmasked and dust-corrected MGEs, respectively.The dashed lines indicate the dust disk edge and the arrows indicate that the dust extends down to the nucleus.
4786 and NGC 5193 dust-corrected MGE solutions created from the combination of HST H-band images and best-fitting GALFIT Nuker models.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⊙,H = 3.37 mag(Willmer 2018), the third is the Gaussian standard deviation along the major axis, and the fourth is the axial ratio.Primes indicate projected quantities.The GALFIT MGE position angle found for NGC 4786 was −17.0 • East of North and for NGC 5193 was 71.2 • East of North.fications of the ALMA synthesized beam, and applied five iterations of the Richardson-Lucy algorithm implemented in the scikit-image Python package (Van der Walt et al. 2014).To weight the line profiles on the oversampled grid scale, we divide the CO flux map so that a pixel in the map is divided into s 2 = 9 subpixels and is normalized so that the subpixels corresponding to the same original pixel have equal fluxes.The model is then block-averaged down to the original ALMA scale and has each of its slices convolved with the model synthesized beam using the convolution implementation in the astropy framework(Astropy Collaboration et al.  2013;Price-Whelan et al. 2018).

Figure 4 .
Figure 4. PVDs along the major axes of both NGC 4786 (above) and NGC 5193 (below) and their respective best-fit models.Columns show ALMA CO(2−1) data (left), models (center), and fractional residuals (right).The PVDs were extracted with a spatial width equivalent to a resolution element for each cube.

Figure 5 .
Figure5.Moment maps for NGC 4786 constructed from the ALMA CO(2−1) data cube (left) and its fiducial model (center, model C).Shown are maps of moments 0, 1, and 2, corresponding to surface brightness, LOS velocity vLOS, and LOS velocity dispersion σLOS.The units for the surface brightness map are mJy km s −1 pixel −1 , and the units for the vLOS and σLOS maps are km s −1 .The systemic velocity of 4621 km s −1 estimated from our dynamical models has been removed from vLOS.Maps of (data-model) residuals are shown in the rightmost column.The coordinate system is oriented such that +x corresponds to East and +y corresponds to North.While the line profile fits have been determined at each pixel of the full disk, the elliptical fitting region used in calculating χ 2 is denoted in the top left panel with a yellow ellipse.The synthesized beam is represented by an open ellipse in the bottom left corner of the same image.

Figure 7 .
Figure 7. CO line profiles extracted from six spatial locations within the block-averaged NGC 4786 and NGC 5193 data cubes, along with their respective fiducial models (models C and F).The x and y positions are given relative to the disk dynamical center, with +y indicating North and +x indicating East.The gray shaded area represents values in the range of data ±1σ, where the 1σ value is from our 3D noise model used in the χ 2 optimization.

Figure 8 .
Figure 8. Plot of log 10 M⋆(r) and log 10 Mgas(r) vs. log 10 r in NGC 4786 (top) and NGC 5193 (bottom) for the three host MGE models and gas mass profiles used.M⋆(r) is calculated via M⋆(r) = rv 2 ⋆,MGE /G, where the v⋆,MGE values have been scaled by their respective √ ΥH values listed in Table2.Mgas(r) is calculated via: rv 2 gas /G.The red dotted line indicates the angular resolution of the ALMA observations, whereas the black dotted line represents the outer edge of the dust disk as measured along the major axis.The gray shaded region indicates the range of MBH determined by our dynamical models using the different input MGE models, while the yellow shaded region enclosed within the black dotted lines indicate the range of the radius of the black hole sphere of influence, rSOI.We have defined rSOI to be the radius where the stellar and BH masses of a given dynamical model are equal.

Table 2 .
Dynamical Modeling ResultsNote-Best-fit parameter values obtained by fitting thin disk dynamical models to the NGC 4786 and NGC 5193 CO(2−1) data cubes.We derive 1σ statistical uncertainties for the parameters of fiducial models C and F, based on a Monte Carlo resampling procedure and list them under the results for models C and F. These models have 363 and 3181 degrees of freedom, respectively.The major axis PA, Γ, is measured east of north for the receding side of the disk.We found the dynamical center of the fiducial model to be at RA=12 h 54 m 32.4115 s , Dec=−06 • 51 ′ 33.′′ 920 for NGC 4786 and RA=13 h 31 m 53.5289 s , Dec=−33 • 14 ′ 03.′′ 546 for NGC 5193.These are within 0. ′′ 001 of the dynamical centers of the other models.The observed redshift, z obs , is used in our dynamical models as a proxy for the systemic velocity of the disk, vsys, in the barycentric frame via the relation: vsys = cz obs and is used to translate the model velocities to observed frequency units.

Table 4 .
NGC 4786 Unmasked and Dust-Masked MGE Components k log 10 I H,k (L⊙ pc −2 ) σ ′ Note-NGC 4786 unmasked and dust-masked MGE solutions created from fitting MGE models in GALFIT to the H-band image.