# The Galaxy Activity, Torus, and Outflow Survey (GATOS). VI. Black hole mass estimation using machine learning

R. Poitevineau, F. Combes, S. Garcia-Burillo, D. Cornu, A. Alonso Herrero, C.

Ramos Almeida, A. Audibert, E. Bellocchi, P G Boorman, A J Bunker, et al.

► To cite this version:

R. Poitevineau, F. Combes, S. Garcia-Burillo, D. Cornu, A. Alonso Herrero, et al. The Galaxy Activ-
ity, Torus, and Outflow Survey (GATOS). VI. Black hole mass estimation using machine learning. As-
tronomy & Astrophysics - A&A, 2025, 693, pp.A311. 10.1051/0004-6361/202347566. hal-05421787

HAL Id: hal-05421787

https://college-de-france.hal.science/hal-05421787v1

Submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Galaxy Activity, Torus, and Outflow Survey (GATOS)

## VI. Black hole mass estimation using machine learning

R. Poitevineau<sup>1,2,*</sup>, F. Combes<sup>1,3</sup>, S. Garcia-Burillo<sup>4</sup>, D. Cornu<sup>1</sup>, A. Alonso Herrero<sup>5</sup>, C. Ramos Almeida<sup>6,7</sup>, A. Audibert<sup>6,7</sup>, E. Bellocchi<sup>8,9</sup>, P. G. Boorman<sup>10</sup>, A. J. Bunker<sup>11</sup>, R. Davies<sup>12</sup>, T. Díaz-Santos<sup>13,14</sup>,

I. García-Bernete<sup>11</sup>, B. García-Lorenzo<sup>6,7</sup>, O. González-Martín<sup>15</sup>, E. K. S. Hicks<sup>16</sup>, S. F. Hönig<sup>17</sup>, L. K. Hunt<sup>18</sup>, M. Imanishi<sup>19,20</sup>, M. Pereira-Santaella<sup>21</sup>, C. Ricci<sup>22,23</sup>, D. Rigopoulou<sup>11,14</sup>, D. J. Rosario<sup>24</sup>, D. Rouan<sup>25</sup>, M. Villar Martin<sup>4</sup>, and M. Ward<sup>26</sup>

*(Affiliations can be found after the references)*

Received 26 July 2023 / Accepted 23 November 2024

## ABSTRACT

The detailed feeding and feedback mechanisms of active galactic nuclei (AGNs) are not yet well known. For low-luminosity AGNs, obscured AGNs, and late-type galaxies, the masses of their central black holes (BH) are difficult to determine precisely. Our goal with the GATOS sample is to study the circum-nuclear regions and, in the present work, to better determine their BH mass, with more precise and accurate estimations than those obtained from scaling relations. We used the high spatial resolution of ALMA to resolve the CO(3–2) emission within ~100 pc around the supermassive black hole (SMBH) of seven GATOS galaxies and try to estimate their BH mass when enough gas is present in the nuclear regions. We studied the seven bright ($L_{AGN}$(14–150 keV) ≥ 10<sup>42</sup> erg/s) and nearby (<28 Mpc) galaxies from the GATOS core sample. For the sake of comparison, we first searched the literature for previous BH mass estimations. We also made additional estimations using the $M_{BH}-\sigma$ relation and the fundamental plane of BH activity. We developed a new method using supervised machine learning to estimate the BH mass either from position-velocity diagrams or from first-moment maps computed from ALMA CO(3–2) observations. We used numerical simulations with a large range of parameters to create the training, validation, and test sets. Seven galaxies had sufficient gas detected, thus, we were able to make a BH estimation from the ALMA data: NGC 4388, NGC 5506, NGC 5643, NGC 6300, NGC 7314, NGC 7465, and NGC 7582. Our BH masses range from 6.39 to 7.18 log($M_{BH}/M_\odot$) and are consistent with the previous estimations. In addition, our machine learning method has the advantage of providing a robust estimation of errors with confidence intervals. The method has also more growth potential than scaling relations. This work represents the first step toward an automatized method for estimating $M_{BH}$ using machine learning.

**Key words.** galaxies: active – galaxies: ISM – galaxies: kinematics and dynamics – galaxies: nuclei – galaxies: spiral

## 1. Introduction

The centers of massive galaxies host supermassive black holes (SMBHs) and black holes (BHs) with masses ranging from $10^6 M_\odot$ to $10^{10} M_\odot$. Those BHs have grown together with their host galaxies over time by accreting matter. The release of this considerable gravitational energy produces shocks and emits powerful radiation in the central region of the galaxies, a phenomenon known as active galactic nuclei (AGNs). There is an empirical relation between the mass of the SMBH ($M_{BH}$) and the bulge mass of their galaxy host, the latter often measured by its central velocity dispersion, the $M_{BH}-\sigma$ relation (e.g., Kormendy & Ho 2013; Shankar et al. 2016, 2019; Marsden et al. 2020). This suggests either common fueling mechanisms and/or a feedback effect of AGNs on their host. This scaling relation demonstrates a symbiosis and regulation from the BH, which transcends the merger and/or accretion history of the BH.

Overall, AGN feedback is often invoked to moderate or even stop star formation in their massive hosts when supernovae feedback is no longer efficient. Indeed, in the standard $\Lambda$CDM cosmological model, simulations of galaxy evolution in their dark matter halos over-predict the number of the most massive galaxies, with large stellar masses. Through halo abundance match-

ing, Behroozi et al. (2013) observe that around $z \sim 3$ massive galaxies become less efficient at forming stars. The comparison between the observed luminosity function and the simulated one for galaxies without feedback, reveals that star formation has not been efficient both at the low and high mass ends: galaxies have not been able to form as many stars as in the simulations. From the observed halo mass function, one can derive the stellar mass-to-halo mass ratio, showing a peak at 20% of the universal baryon to dark matter ratio, meaning that >80% of the baryonic matter resides outside of galaxies. While at the low mass end, supernovae feedback can eject most of the baryons, only AGN feedback is efficient enough at the high mass end. While at high redshift, dense and cold gas can be accreted efficiently in dark matter halos, with shallow potential wells (Dekel et al. 2009), at low redshift, gas is shock-heated when entering massive galaxies and deep potential wells, and only sparse cold gas is infalling. However, the hot gas could still cool down in the center of the structure and efficiently feed star formation. Then AGN feedback is still necessary to moderate star formation (Fabian 2012).

AGN feedback can be classified into two modes: radiative (quasar) mode and kinetic (radio) mode. In the first case, luminous AGNs are powered by thin accretion disks. If the AGN bolometric luminosity is significant with respect to the Eddington luminosity, which is proportional to the BH mass, the radiation pressure creates a wind from the accretion disk. This is the quasar mode, where the AGN luminosity is higher than 1%

<footnote>* Corresponding author; remi-poitevineau@epfl.ch</footnote>
Eddington (Fabian 2012). In the second case, for low-luminosity AGNs, the pressure and dragging of the relativistic jet and consequent shocks heat the gas and a bubble/cocoon is inflated by the jet. This is the radio mode.

Estimating the SMBH mass is therefore fundamental to improving our understanding of their formation, in symbiosis with the galaxy growth. One of the best ways to measure its mass is to measure its gravitational influence on its surroundings, through the kinematics of the matter inside its sphere of influence. Another way is to rely on the empirically determined $M_{BH}-\sigma$ scaling relation, but with more scatter at the low-mass end (e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000).

For nearby galaxies, it is possible to measure the BH mass ($M_{BH}$) from the stellar kinematics, as in M31 (Bender et al. 2005). The nuclear stellar disk, inside the sphere of influence of the black hole (SoI), is dominated by Keplerian rotation. The kinematics of central ionized gas can also be used for $M_{BH}$ determinations. It offers multiple advantages. First, ionized gas emission lines are much easier to detect than stellar absorption lines or even molecular lines. Secondly, the velocity dispersion of gas is much lower than that of stars and it is easier to model its rotating Keplerian disk (e.g., Kormendy & Ho 2013). However, the gas dynamics has some major drawbacks. Enough gas should be detected to properly sample the nuclear dynamics, and it should be well distributed within the SoI of the BH. Gas is modeled as a collisional fluid, meaning that non-gravitational perturbations must be taken into account. Perturbations, such as outflows, can affect the gas kinematics locally and must be disentangled to obtain the $M_{BH}$. The true velocity dispersion must be distinguished from a beam-smeared rotational velocity gradient to properly correct and obtain the circular velocity for the $M_{BH}$. This is why it is important to compare $M_{BH}$ estimations made with ionized or molecular gas with masses obtained using other methods.

Over the past few years, the domain of machine learning (ML) has vastly expanded in the astrophysical community. It is a very powerful and flexible tool for data analysis that has demonstrated excellent performance on a wide variety of tasks using complex data. One of the main domains of the application of ML is image analysis. Convolutional neural networks (CNN) are efficient and powerful artificial neural network (NN) architectures to handle images. For example, they are used to classifying large numbers of galaxies (Huertas-Company et al. 2020), in view of future missions such as Euclid. They are also used to classify supernovae (Lochner et al. 2016) and find strong lensing features (Lanusse et al. 2018), transient phenomena (Mahabal et al. 2019), photometric redshifts (Pasquet et al. 2019), or star clusters (Castro-Ginard et al. 2019). They have however not yet been used to estimate the mass of SMBHs.

In this paper, we present a new ML method to estimate the mass of SMBH in galaxies. The masses obtained with this new method are compared with SMBH masses found in the literature and with other estimations made using the fundamental plane of BH activity (FPBH) and the $M_{BH}-\sigma$ relation. We describe the galaxy sample and our selection criteria in Section 2. Section 3 focuses on the different methods to estimate the $M_{BH}$. We first introduce the $M_{BH}-\sigma$ and the reverberation mapping methods that were used to make previous estimations and calibrations. We also discuss the FPBH method used to complete the diverse masses from the literature. Section 4 presents how we combine ALMA data cubes, numerical simulations, and a dedicated ML approach to estimate $M_{BH}$ values more precisely from the CO(3–2) gas dynamics. A large number of simulated models are used to evaluate the quality of our estimations. In Section 5,

the results found with our ML method are presented. Finally, Section 6 and 7 summarize and discuss our findings, along with possible future improvements to our methodology.

## 2. Sample

The GATOS international collaboration has for aims to study topics related to the physics taking place in the nuclear region of AGN. These encompass the gas flow cycle, the emission of polar dust, the properties of the torus and obscuring material, and the interplay between star formation activity and AGN phenomena. Here, we study the ten galaxies from the GATOS core sample studied by García-Burillo et al. (2021). The core sample is composed of nearby galaxies (<28 Mpc) and luminous nuclei, $L_{AGN}$ (14–150 keV) $\ge 10^{42}$ erg/s. The upper limit on the distance is set to have a sufficient spatial resolution to study molecular tori as small as ~10 pc in radius while the limit on the luminosity is set to avoid overlap with the ongoing ALMA surveys of nearby Seyferts such as the NUGA (Nuclei of Galaxies survey) (e.g., Combes et al. 2019; Audibert et al. 2019).

The present work uses the CO(3–2) data cubes from ALMA observations. All the details of the GATOS core sample and their observations are given in García-Burillo et al. (2021). We do not consider the three additional GATOS galaxies NGC 1068 (e.g., García-Burillo et al. 2019), NGC 1365 (e.g., Combes et al. 2019), and NGC 3227 (e.g., Alonso-Herrero et al. 2019).

Not all ten data cubes are appropriate for the determination of the central kinematics and $M_{BH}$ estimation. Two of the galaxies, NGC 6814 and NGC 7213 do not present enough CO emission in the very center; their central pixels are empty or with some emission with an signal-to-noise ratio (S/N) below 2. For another galaxy, moment maps are not usable for a $M_{BH}$ estimation with the ML approach; NGC 4941 has too little gas detected, and its distribution is too sparse. There remain therefore 7 galaxies: NGC 4388, NGC 5506, NGC 5643, NGC 6300, NGC 7314, NGC 7465, and NGC 7582.

## 3. Classical methods for estimating $M_{BH}$

### 3.1. Dynamical estimation

The estimation of SMBH mass through dynamical models involves studying the movements of celestial objects, such as stars or gas, around the SMBH. These observed motions, often referred to as tracers, serve as constraints for the dynamical models. These models discern the contribution of the BH to the potential from that of the surrounding galaxy, thus facilitating the inference of the SMBH mass by fitting these models to observations (Kormendy & Richstone 1995; Ferrarese & Merritt 2000; Kormendy et al. 2004; Gebhardt et al. 2000; Kormendy & Ho 2013).

A notable constraint of this methodology is its reliance on resolving the kinematics of objects situated near the SMBH. Consequently, it becomes impractical for galaxies located at considerable distances. Therefore, alternative methodologies are required to estimate SMBH masses in such distant galaxies.

Across various dynamical models, an essential consideration is the necessity for a realistic representation of the gravitational potential of the galaxy. Additionally, the data characterizing the tracers motions must be robust enough to detect the presence of a BH (Gebhardt et al. 2003; Ferrarese & Ford 2005; Kormendy & Ho 2013). Among the multitude of intricate processes connected to SMBHs and their immediate surroundings, the study of stellar phenomena proves more amenable to analysis
and modeling. This feasibility arises from stars often being treatable as point masses, predominantly influenced by gravity. The substantial mass disparity between stars and $M_{BH}$ simplifies the problem further, allowing stars to be considered test masses on intermediate scales. Here, they operate within the $M_{BH}$ potential, where gravitational interactions among stars are negligible, yet they are distant enough to avoid significant tidal effects from the $M_{BH}$.

The molecular gas surrounding the SMBH serves as another dynamical tracer for estimating the $M_{BH}$. Gas-dynamical mass measurements present distinct advantages compared to stellar-dynamical modeling. The simplicity inherent in gas-dynamical modeling stems from the gas Keplerian rotation within a dynamically cold disk. Unlike the intricate orbit-based computations necessary for treating stars, analyzing gas dynamics in a rotating disk proves computationally less demanding. The reduced complexity allows neglecting factors such as orbital anisotropy, triaxiality, or the influence of the dark matter halo. The gas-dynamical approach assumes a thin, rotating disk following circular orbits within the principal plane of the galaxy potential (Marconi et al. 2001, 2003). This methodology aims to calculate a model velocity field aligning with observed velocities, velocity dispersions, and the line emission's surface-brightness distribution. The gravitational potential of the galaxy comprises contributions from stars, measured by the projected stellar surface brightness, and assumes a mass-to-light ratio along with the presence of a BH (Melchior & Combes 2011; Combes et al. 2019).

However, gas dynamics introduce several complications (Kormendy & Ho 2013). Properly sampling the BH SoI necessitates distributing gas across radii. Simultaneously, the gas kinematics must display sufficient orderliness for interpretation. Furthermore, gas differs from stars as it behaves like a collisional fluid, responsive to non-gravitational disturbances such as turbulence, shocks, radiation pressure, and magnetic fields. Another challenge in gas dynamics involves dust absorption, potentially rendering the gas distribution opaque. Consequently, assumptions about observing through the gas in projection may not hold under such circumstances. Each galaxy requires meticulous examination to confirm whether the gas has reached an equilibrium configuration, influenced mainly by gravitational effects, considering the opacity caused by dust absorption.

To investigate the gas dynamics of CO(3–2), a useful source of information is the position-velocity diagram (PVD) obtained along the major axis, which is widely employed for this purpose, for instance by the WISDOM project (e.g., Onishi et al. 2017; Davis et al. 2017, 2018).

To estimate the $M_{BH}$ from the PVD, we fit the sum of the contributions from the various galaxy central components: bulge, BH, disk, and possible bar and stellar nucleus.

To model the stellar distribution in each galaxy, we use the Salo et al. (2015) Galfit (Peng et al. 2010) decomposition from the Spitzer Survey of Stellar Structure in Galaxies 3.6 µm infrared images (Muñoz-Mateos et al. 2015). Two of our galaxies, NGC 5643 and NGC 6300, were not in the Salo et al. (2015) sample, so we made the Galfit decomposition ourselves. The physical characteristics of the different galaxy components of our sample are shown in Table 1.

The models employed to fit the galaxies exhibit symmetry, while the observed PVD does not. To address this asymmetry, we adopt a symmetrization approach by generating a symmetrical version of the observed PVD by mirroring it about its central pixel, followed by the summation of the original and mirrored versions. An illustration of the process employed to obtain the symmetrized PVD of NGC 4388 is provided in Fig. 1.

<table><caption>**Table 1.** Stellar mass components.</caption><thead><tr><th>ID</th><th>Comp</th><th>R<sub>*</sub> (kpc)</th><th>M<sub>*</sub>(10<sup>10</sup> M<sub>☉</sub>)</th><th>Sersic n</th><th>B/T</th></tr></thead><tbody><tr><td rowspan="2">NGC 4388</td><td>Disk</td><td>1.68</td><td>2.44</td><td rowspan="2">1.0</td><td rowspan="2">0.84</td></tr><tr><td>Nucleus</td><td></td><td>0.46</td></tr><tr><td rowspan="4">NGC 4941</td><td>Bulge</td><td>0.67</td><td>0.50</td><td rowspan="2">11.5</td><td rowspan="4">0.24</td></tr><tr><td>Disk</td><td>1.13</td><td>0.95</td></tr><tr><td>Bar</td><td>10.92</td><td>0.66</td></tr><tr><td>Nucleus</td><td></td><td>0.01</td></tr><tr><td rowspan="4">NGC 5643</td><td>Bulge</td><td>0.20</td><td>0.89</td><td rowspan="2">2.9</td><td rowspan="4">0.13</td></tr><tr><td>Disk</td><td>9.12</td><td>4.61</td></tr><tr><td>Bar</td><td>3.75</td><td>1.17</td></tr><tr><td>Nucleus</td><td></td><td>0.01</td></tr><tr><td rowspan="4">NGC 6300</td><td>Bulge</td><td>4.07</td><td>0.875</td><td rowspan="2">1.7</td><td rowspan="4">0.19</td></tr><tr><td>Disk</td><td>44.11</td><td>1.375</td></tr><tr><td>Bar</td><td>27.80</td><td>2.240</td></tr><tr><td>Nucleus</td><td></td><td>10<sup>-4</sup></td></tr><tr><td rowspan="2">NGC 7314</td><td>Disk</td><td>1.70</td><td>2.105</td><td rowspan="2">1.0</td><td rowspan="2">0.0</td></tr><tr><td>Nucleus</td><td></td><td>0.012</td></tr><tr><td rowspan="3">NGC 7465</td><td>Bulge</td><td>0.34</td><td>0.486</td><td rowspan="2">1.0</td><td rowspan="3">0.19</td></tr><tr><td>Disk</td><td>2.51</td><td>0.630</td></tr><tr><td>Bar</td><td>1.35</td><td>1.124</td></tr><tr><td rowspan="3">NGC 7582</td><td>Bulge</td><td>0.19</td><td>1.132</td><td rowspan="2">3.7</td><td rowspan="3">0.35</td></tr><tr><td>Disk</td><td>2.24</td><td>1.989</td></tr><tr><td>Bar</td><td>12.54</td><td>0.135</td></tr></tbody></table>

**Notes.** Stellar mass components from the Salo et al. (2015) Galfit decomposition from $S^4G$ 3.6 µm except for NGC 5643 and NGC 6300 where the Galfit decomposition is our own fit from the WCF3 Hubble Legacy Archive image.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Example of NGC 4388 symmetrized PVD.</figcaption></figure>

Since the symmetrization tends to keep the highest gradients on the PVD, we can expect that the $M_{BH}$ estimations might be biased high.

In a Keplerian rotating disk, the gas rotation velocity is inversely proportional to the square root of the radius. When the gas is just entering the SoI of the BH, the main contribution to the velocity is the BH, but there is also some contribution from the bulge stars. Thus, the measured velocity, corrected from velocity dispersion, can put an upper limit on the $M_{BH}$. If no gas is observed in the SoI, an estimation can be done from the closest gas detected instead, assuming an accurate estimation of the bulge mass. Large error bars can be expected if the gas is far from the SoI, and/or if the bulge mass is not precisely known.
## 3.2. $M_{BH}$-$\sigma$ relation

The tight correlation between BH and bulge masses drawn from the $M_{BH}$-$\sigma$ relation (e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000) suggests that there is a mechanical and/or radiative feedback from the AGN on the star formation in galaxies, all along their growth through accretion and mergers. The origin of this coupling is not yet fully understood, but there are analytic and numerical models created to constrain the physical mechanism and reproduce the $M_{BH}$-$\sigma$ relation (e.g., Silk & Rees 1998; Granato et al. 2004). Nevertheless, the observational $M_{BH}$-$\sigma$ relation is not well constrained for low-mass galaxies ($M < 10^9 M_\odot$). The tight correlation is mostly derived from biased galaxy samples with a high $M_{BH}$, which have a larger SoI easier to resolve with common angular resolutions (e.g., Silk & Rees 1998; Costa et al. 2014). In this work, we use the Kormendy & Ho (2013) correction of the $M_{BH}$-$\sigma$ relation:

$$ \log M_{BH}(M_{\odot}) = 8.5 \pm 0.05 + (4.41 \pm 0.29) \log \left( \frac{\sigma}{200 \text{ km/s}} \right), \quad (1) $$

where $\sigma$ is the stellar central velocity dispersion. For our sample, it was taken from the Hyperleda compilation (e.g., Makarov et al. 2014). This relation has an intrinsic error of 0.3 dex.

## 3.3. Reverberation mapping

The radii of BLRs are several light days to weeks across; thus, they are too small to be spatially resolved. However, they can be resolved in time with the reverberation mapping technique (Blandford & Payne 1982; Bahcall et al. 1972; Lyutyi & Cherepashchuk 1972; Cherepashchuk & Lyutyi 1973). The SMBH accretion disk shocks emit UV radiation that heats the BLR gas clouds, ionizing them and exciting emission lines. Thus, when the UV continuum source varies, the emission lines vary as well (i.e., “reverberate”) but with a delay of $R_{BLR}/c$. Therefore, the BLR radius can be estimated by measuring the time delay between the change in emission from the central source and its impact on the BLR lines. Then, assuming that the width of the BLR emission lines is due to the gravitational influence of the central SMBH, we can derive the following relation from the Virial theorem:

$$ M_{BH} = \frac{f (\Delta V^2) R}{G}, \quad (2) $$

where $M_{BH}$ is the mass of the central BH, $G$ is the gravitational constant, $\Delta V$ is the virial velocity of the BLR, $R$ is the BLR radius, and $f$ is a factor depending on the BLR geometry. The fact that the BLR velocity is dominated by the BH mass (virial velocity and gravitational redshift) and not by possible gas outflows, was validated by Gaskell (1988).

Most of the uncertainties of this method come from the factor, $f$. Kollatschny & Zetzl (2011) observed that the BLR structure changes with the emission-line turbulent dispersion and rotational width, implying that the value of $f$ depends on these two factors. The findings of Gaskell (2008, 2011) indicate another inherent constraint, which suggests that the continua of AGN can exhibit localized flares in off-center regions situated relatively close to the inner boundary of the BLR radius.

Even in the absence of the expensive and time-consuming monitoring observations, it is possible to use this method, thanks to the AGN-Luminosity vs BLR-size scaling relation, calibrated from reverberation mapping measurements (Wandel et al.

1999). Refining the method with a large reverberation mapping database, and comparing it to the other methods, Peterson et al. (2004) find that the precision of the method is ~30%, comparable to the gas or stellar kinematical method.

This can be explained by considering that the reverberation mapping technique is notably influenced by the inclination angle of the BLR region, contrary to the gravitational redshift method. As demonstrated in the study by Liu et al. (2017), the substantial disparity between gravitational and virial mass estimates can be reconciled by adjusting the parameter $f$ to a value between 8 and 16 for the virial measurement. This adjustment of $f$ depends on radiation pressure and accretion rate, in addition to the inclination of the BLR. The conventional value of $f \sim 1$ corresponds to an inclination angle of around 30 degrees.

## 3.4. Single epoch method

The single epoch method uses the relation between the BLR size and the AGN optical/UV continuum luminosity empirically found from reverberation mapping (Peterson et al. 2004; Kaspi et al. 2007; Bentz et al. 2009), as well as the tight correlation between the continuum luminosity and that of broad emission lines (e.g., Shen et al. 2011). With these considerations, the $M_{BH}$ can be expressed as

$$ \log \left( \frac{M_{BH}}{M_{\odot}} \right) = a + b \log \left( \frac{L}{10^{44} \text{ erg s}^{-1}} \right) + c \log \left( \frac{\text{FWHM}}{\text{km s}^{-1}} \right), \quad (3) $$

where the coefficients $a$, $b$, and $c$ are empirically calibrated against local AGN with reverberation mapping masses or using different lines. Then, $L$ and FWHM are the line luminosity and width. At low redshift, below 0.75, $H\beta$ is the most commonly used line. For redshift between 0.75 and 2, the Mg II $\lambda2800$ line is a common choice (McLure & Jarvis 2002). And for $z \ge 2$, the best line to use is C IV $\lambda1549$ Vestergaard (2002). Since their first estimation, the relation for these lines had been revised by Vestergaard & Peterson (2006), Wang et al. (2009), Shen et al. (2011), Xiao et al. (2011) and others.

The single epoch method has the advantage of being inexpensive in telescope times since it can estimate the SMBH mass from just a spectrum. However, it is a more indirect method than reverberation mapping and does not contain as many indicators that the measured velocities are virialized ones. This method has an uncertainty of ~0.3 dex.

## 3.5. Fundamental plane of BH activity

The FPBH is a strong correlation between the $M_{BH}$, its 5 GHz radio continuum luminosity, and its 2–10 keV X-ray power-law continuum (Merloni et al. 2003). When accretion disks go through episodes of low-luminosity advection-dominated accretion flow, they are observed to be in a hard state (hard X-ray radiation). In this state, they are often accompanied by relativistic jets (Gallo et al. 2003; Narayan 2005). During this state of inefficient radiation, the accreting material stops flowing in a disk at some distance from the accreting object. This new configuration results in the creation of a hot, optically-thin gas corona around the BH. Photons, from the now truncated disk, are Compton-up scattered to tens to hundreds of keV when they traverse this corona. Furthermore, the synchrotron mechanism in the relativistic jets produces emission primarily in the radio band (Blandford 1984). The entire broadband combined emissions of the accretion disk, corona, and relativistic jets range from the radio to the X-ray bands. This emission can also be
dominated by the jet itself, particularly if accretion rates are higher than ~0.01$M_{Edd}$ (Falcke et al. 2004). A fainter thermal component emanating from the shortened accretion disk may also be identifiable alongside the primary non-thermal corona emission. The jet is linked to the accretion process, and it has been shown that the luminosity of the X-ray and radio emissions are correlated and the disk-jet mechanism is independent of the $M_{BH}$. By considering these arguments, Merloni et al. (2003) and Falcke et al. (2004) probed a large sample of galactic BH and SMBH and found a strong correlation between the X-ray luminosity (between 2 and 10 keV), the radio luminosity (5 GHz) and the $M_{BH}$. This correlation is called the fundamental plane of BH activity. According to Gültekin et al. (2019) the $M_{BH}$

$$ \log\left(\frac{M_{\mathrm{BH}}}{10^8 M_{\odot}}\right) = 0.55 \pm 0.22 + (0.19 \pm 0.10)\log\left(\frac{L_R}{10^{38} \mathrm{erg/s}}\right) + (-0.59^{+0.16}_{-0.15})\log\left(\frac{L_x}{10^{40} \mathrm{erg/s}}\right), \quad (4) $$

where $L_x$ is the 2–10 keV X-ray power-law continuum and $L_R$ the 5 GHz radio luminosity calculated as $L_R = \nu L_\nu = (\nu = 5 \mathrm{GHz}) \times 4\pi D^2 F_5$, $F_5$ being the 5 GHz radio flux density. This relation has an intrinsic error of ~0.40 dex.

This FPBH scaling relation is a useful method to distinguish between X-ray binaries, intermediate mass BH, and SMBH, and to determine the $M_{BH}$ in a Type 2 AGN, or in a host galaxy with disturbed morphology, where the $M_{BH}-\sigma$ method is unusable. The only requirement to estimate the $M_{BH}$ is to have X-ray and radio data of sufficient quality to resolve the central source and also confirm its spectral characteristics (i.e., that it is indeed in the hard state). The data must have a high angular resolution to avoid contamination from other sources and the host. This is especially important to distinguish low luminosity AGN from X-ray binaries.

## 4. ML approach to estimate $M_{BH}$

We have chosen to employ ML techniques to estimate the BH mass using the ALMA data cubes of the GATOS core sample. These 3D data cubes can be effectively transformed into 2D images, enabling us to apply more easily CNNs for analysis. Although everyday data tend to be more complex, CNNs have proven to be highly effective in such cases, leading us to expect their suitability for astrophysical data as well. This study aims to serve as a proof of concept that ML techniques can be used to infer $M_{BH}$ from observational data.

As of now, our sample consists of only seven galaxies, which implies that applying an ML method to such a small dataset may not yield the most optimal results. However, we anticipate the availability of a larger sample in the future. Additionally, the presence of data in the ALMA archives presents further opportunities for analysis and investigation, thus motivating our efforts toward automating the BH mass estimation process. Although scaling relations offer a certain level of automation, their intrinsic errors restrict the full exploitation of the high-resolution data. The ML method presented in this study serves as a proof of concept that ML can be applied to ALMA data to estimate $M_{BH}$ with physical significance and is a first step toward an automated approach that leverages the data resolution to achieve precise $M_{BH}$ estimations.

### 4.1. Supervised learning method

An artificial NN can be seen as a sequence of mathematical operations, organized in layers and composed of adjustable weights, that map an input to an output. In a supervised learning method, a dataset of inputs with known outputs, called the training set, is used to parameterize the weights so the NN infers the mapping function. It is done by feeding the training set to the NN and comparing its outputs to the expected ones. A “loss function” is used to compare the NN output to the expected one. The loss function serves as a measure of how good the NN is at predicting the expected outcome. The more accurate the output, the smaller the loss function. The goal of the training phase is to find a set of weights that minimize the loss function. To do so, the loss function is computed for each input in the training dataset. Then, the gradient of the loss function is used to slightly modify the weights via a backpropagation algorithm. This process is repeated for each input in the training data set. One pass of all the training data through the NN is called an epoch.

Overfitting is a common problem when training an artificial NN with a supervised ML method, and alerts when the NN infers a mapping function that is too specific to the training set, making it unable to perform accurately against unseen data. To avoid falling into such a regime, another dataset, called the validation set, is used. The validation set has the same statistic distribution as the training set but is never used to update the network weights. During the training, the loss function values of the training set and the validation set are regularly monitored to ensure that the network is still learning a generalized mapping function. Finally, after the training process is done, a third dataset, the test set, is used to estimate the performance of the NN on unseen data. This practice aids in preventing the training process from prematurely halting upon an overestimation of model performance on the validation set, which may occur due to limited sample sizes or biased sample selection.

The principle of our approach is to run a large number of numerical simulations, fitted to each studied galaxy. With this method, the gravitational potential of the galaxy disk is assumed to be known from HST images, within some uncertainties. The simulations are subsequently conducted by systematically varying the parameters within the defined uncertainty ranges. This includes the mass and scales of the bulge (plus possible nucleus) and disk (plus possible bar) components, as well as the characteristics of the gas and stellar distributions (notably Sersic index). Inclination and position angles are fixed to the values estimated from the observations (García-Burillo et al. 2021), except in some specified cases, where the nuclear disk is tilted with respect to the large-scale disk. The possible mass of the SMBH is varied in a wide range, compatible with the observations. Numerous numerical simulations with different $M_{BH}$ are used to train an artificial NN model for each galaxy. After the training is done on the simulations, the trained model can then be used to infer the estimated $M_{BH}$ from the real observation.

The sample size of seven galaxies is insufficient for training a NN to generalize across different galaxies due to the high variability and complexity in their physical properties. A larger and more diverse dataset is essential to capture the wide range of conditions present in the universe, which is not feasible with just seven galaxies. However, this limited sample size makes it practical to train a separate NN for each galaxy. By focusing on individual galaxies, we simplify the problem and can determine if our methodology works effectively on a smaller scale. Training separate networks for each galaxy involves creating distinct training, validation, and test sets for each one. For each observed
galaxy, we generated 25 000 simulated data cubes, dividing them into a training set of 16 500 simulations, a validation set of 3500 simulations, and a test set of 5000 simulations. This approach allows us to study the behavior of the networks more precisely, as we can evaluate their performance in a controlled and consistent manner. By doing so, we can identify any potential issues or improvements needed in our model before attempting to generalize the approach to a more extensive and diverse set of galaxies. If successful, this method can serve as a foundation for scaling up to a broader application, ultimately aiming to develop a generalized model capable of accurately inferring $M_{BH}$ across various galaxies.

## 4.2. Numerical simulations

To simulate our data cubes, we used the model described in Melchior & Combes (2011) and Combes et al. (2019). We simulate gas particles in a potential that depends on the physical parameters of the galaxies and their mass components. Since we simulate only the nuclear disk with a radius of typically ~100 pc, the dark matter contribution to the gravitational potential is negligible. Our models are axisymmetric homogeneous gas disks (e.g., Miyamoto & Nagai 1975) with gas particles on nearly circular orbits, with low-velocity dispersion, corresponding to thin cold disks. The sizes of the modelled disks are taken from the observed CO(3–2) disks. The velocity distribution of the gas particles is set such that the Toomre Q parameter of the disk is equal to one and the ratio between tangential and radial velocity dispersion is taken from the epicyclic theory (e.g., Toomre 1964). To ensure sufficient statistics, we fixed the number of particles at $10^6$.

To compare the model with observational data, we generated data cubes by projecting the model onto the sky using the best-fit large-scale inclinations and position angles. We then computed the line-of-sight velocity distribution. The pixel size in the data cubes was chosen to match that of the observed data (ranging from 4 to 7 pc, depending on the galaxy), with velocity channels of 10 km/s. The data were smoothed to match the observed beam, resulting in a cube of dimensions (250, 250, 60) for each simulation. Since the gas distribution exhibits asymmetry and patchiness, which affect the mass-weighted velocity within each observed beam, we applied a normalization process to the model cube. Specifically, we normalized the model cube pixel by pixel in the 2D projection using the zeroth moment map of the CO observations. This normalization serves as a multiplicative filter for our homogeneous gas disks. Consequently, each CO spectrum at every position in the model is normalized to the observed integrated flux at that corresponding position.

The different parameters of the simulations are: the mass, radius, and height of the galaxy bulge ($M_{bulb}$, $R_{bulb}$, $H_{bulb}$); the corresponding properties of the stellar disks ($M_{gal}$, $R_0$, $H_t$); the gas-to-stellar mass ratio within the stellar disk (RAPM); the gas disk radius, height, inclination angle, and position angle ($R_{gas}$, $H_{gas}$, IA, PA) and the SMBH mass ($M_{BH}$).

In the case of NGC 5643, it is observed that the PA changes at approximately 30 pc. To model this behavior, we introduced a second PA that replaces the initial one after the 30 pc mark.

For each simulation, all the parameters are taken randomly from a uniform distribution between the limits with physical meaning (see Table 1), except the IA and the PA. We can deduce from the moment maps a range for the radial scale of the gas $R_{gas}$ and since the molecular disk is thin, we choose to keep the ratio $H_{gas}/R_{gas} < 0.1$. For the stellar disk $R_0$ and $M_{gal}$ are large with respect to the scale probed by our observations/simulations,

and their possible variations are small. The IA and PA are set constant in all the simulations for a galaxy. They were taken as estimated by García-Burillo et al. (2021) based on the gas kinematics.

## 4.3. Preparing the inputs

The GATOS CO(3–2) molecular gas observations are in the form of spectroscopic data cubes which are 3D arrays that combine spatial and spectral information. Two of the dimensions represent the spatial coordinates and define the position in the sky. The third dimension represents the intensity or flux of light at each spatial position and each wavelength or frequency.

Gathering information directly from spectroscopic data cubes can be a challenging task, especially when these data cubes are utilized as inputs in ML processes (Hartley et al. 2023). The introduction of a third dimension into the input data significantly amplifies data complexity and can lead to considerable increases in computational requirements. To mitigate the resulting computational load and efficiently extract valuable insights from the data cube, it becomes essential to optimize and fine-tune the ML architectures employed.

However, we can take an easier approach and pre-process the data cubes to put the meaningful information more straightforwardly in the shape of 2D images. This will make it easier for the NN to learn how to estimate the $M_{BH}$ from the simulated data cubes. Two commonly used representations to study the CO(3–2) gas dynamics are the PVD and the moment maps (see Fig. A.7).

The PVD allows the gas velocity and concentration as a function of the radius in a chosen direction to be seen in detail. Since the PVD is used to better see the effect of the SMBH on the gas dynamics, we select the PVD along the axis with the higher velocity gradient at the smallest radius possible. Nevertheless, this representation has certain limitations. The PVD solely provides insight into a particular direction of the gas dynamics, thereby leading to the partial depiction of information. Consequently, the selection of a specific direction can influence the bias in the information presented in the PVD. Furthermore, studying the gas in only one direction makes the PVD more sensitive to a lack of gas at a specific radius, such as the center of NGC 4388 (Fig. A.1).

In addition to studying the gas in one specific diagram, we consider a representation that encapsulates more information about the entire gas dynamic distribution by calculating the first two moment maps. The zeroth-moment map is the integrated intensity over the spectral axis, proportional to the gas surface density,

$$M_0 = \int I_{\nu} d\nu. \qquad (5)$$

The 1st moment represents the velocity, V, weighted by the intensity,

$$M_1 = \frac{\int V I_{\nu} d\nu}{M_0}. \qquad (6)$$

The PVD and the first moment maps are calculated for each simulated cube. Finally, moment by moment, and for the PVD, we normalize the value of each pixel by dividing it by the highest one. The PVD and first-moment maps of the galaxies in our sample are presented in Appendix A. Initially, the estimation of the SMBH mass is conducted with the PVD as input and then the first-moment map is used as the input.
<table><caption>**Table 2.** NN architecture for $M_{BH}$ estimations with PVD.</caption><thead><tr><th>Layer</th><th>Filters</th><th>Size</th><th>Stride</th></tr></thead><tbody><tr><td>Conv1</td><td>16</td><td>7×7</td><td>2</td></tr><tr><td>Conv2</td><td>20</td><td>7×7</td><td>1</td></tr><tr><td>Conv3</td><td>20</td><td>7×7</td><td>2</td></tr><tr><td>Conv4</td><td>24</td><td>7×7</td><td>1</td></tr><tr><td>Conv5</td><td>64</td><td>5×5</td><td>2</td></tr><tr><td>Conv6</td><td>80</td><td>5×5</td><td>1</td></tr><tr><td>Conv7</td><td>80</td><td>5×5</td><td>2</td></tr><tr><td>Conv8</td><td>128</td><td>5×5</td><td>1</td></tr><tr><td>Max pooling</td><td></td><td>3×3</td><td>2</td></tr><tr><td>Dense1</td><td>2048</td><td></td><td></td></tr><tr><td>Dense2</td><td>1024</td><td></td><td></td></tr><tr><td>Dense3</td><td>512</td><td></td><td></td></tr><tr><td>Dense4</td><td>1</td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 3.** NN architecture for $M_{BH}$ estimations with moment maps.</caption><thead><tr><th>Layer</th><th>Filters</th><th>Size</th><th>Stride</th></tr></thead><tbody><tr><td>Conv1</td><td>32</td><td>11×11</td><td>4</td></tr><tr><td>Max Pooling 1</td><td></td><td>3×3</td><td>2</td></tr><tr><td>Conv2</td><td>64</td><td>5×5</td><td>1</td></tr><tr><td>Max Pooling 2</td><td></td><td>3×3</td><td>2</td></tr><tr><td>Conv3</td><td>64</td><td>5×5</td><td>1</td></tr><tr><td>Max Pooling 3</td><td></td><td>3×3</td><td>2</td></tr><tr><td>Conv4</td><td>64</td><td>5×5</td><td>1</td></tr><tr><td>Max Pooling 4</td><td></td><td>3×3</td><td>2</td></tr><tr><td>Conv5</td><td>96</td><td>3×3</td><td>1</td></tr><tr><td>Conv6</td><td>96</td><td>3×3</td><td>1</td></tr><tr><td>Conv7</td><td>96</td><td>3×3</td><td>1</td></tr><tr><td>Max Pooling 5</td><td></td><td>3×3</td><td>2</td></tr><tr><td>Dense1</td><td>2048</td><td></td><td></td></tr><tr><td>Dense2</td><td>1024</td><td></td><td></td></tr><tr><td>Dense3</td><td>512</td><td></td><td></td></tr><tr><td>Dense4</td><td>1</td><td></td><td></td></tr></tbody></table>

### 4.4. ML architecture

To make our $M_{BH}$ estimations via ML, we used a CNN-supervised training method (e.g., LeCun et al. 2015). Our NNs were developed using Python's library Keras (Chollet 2015). The PVD and the moment maps represent the information differently, so to get the best out of their distinct characteristics, we need different CNN architectures. Thus, we will use two different architectures to make the $M_{BH}$ estimations: one taking the PVD as input and one taking the first-moment maps.

The NN taking the PVD as inputs is composed of a succession of eight convolution layers, a max pooling layer, and at the end, four dense layers (i.e., Table 2). The NN that takes the first moment maps as inputs is composed of a succession of convolution and max pooling layers followed by four dense layers (i.e., Table 3). In both of these architectures, the fourth dense layers return the $M_{BH}$.

Dropout was added to the second dense layer in both architectures with a drop rate of 20%. It means that in this layer, 80% of the neurons are randomly selected to remain active at each path through the network. Usually, the dropout technique is used as a regularization method during training to help prevent over-fitting (Srivastava et al. 2014). It can also improve the predictive

performance of NNs in several tasks. Another way to use dropout is to keep it activated at inference time. At each inference, the dropout will randomly select a different combination of neurons, resulting in a different estimation. It is equivalent to an approximation of the probabilistic deep Gaussian process. Thus, making multiple inferences for the same input will not give one value but a probability distribution on the predicted value, similar to what is produced by a Markov chain Monte Carlo (MCMC) process. This distribution can be studied to characterize the model uncertainties for that input.

We used the root mean square error function (RMSE) as loss function. It is a common loss function used in regression problems:

$$RMSE = \sqrt{\frac{\sum_{i=1}^{n}(Y_{\text{pred},i} - Y_{\text{real},i})^2}{n}}, \quad (7)$$

with *n* the number of predictions made, $Y_{\text{pred},i}$ the predicted values, and $Y_{\text{pred},i}$ the actual values. To help stabilize the weights in our models during the training, we normalized the output targets (i.e., the $M_{BH}$) by dividing them by the highest value among the output targets.

We trained our models during 2000 epochs and took the weights of the epoch with the lowest RMSE on the validation set. We used an initial learning rate of 0.0001 and the Adam optimizer (Kingma & Ba 2014). The training was made with a mini-batch size of 128 images.

### 4.5. Estimating the reliability of models

Multiple statistical indicators can be used to evaluate the quality of a regression model. These indicators can be used as metrics to help compare the training results of the different models. A first indication of the quality can be obtained by looking at the value of the loss function on the test set. However, studying the loss function alone is not enough to evaluate a model because it can be highly degenerate. Furthermore, the RMSE function goes from 0 (best fit possible) to +∞ (worst fit possible). The lack of an upper bound can make the RMSE interpretation difficult but in our case, the error directly measures the average difference between the estimated $M_{BH}$ and the actual ones.

A more reliable indicator would be the coefficient of determination (Wright 1921)<fnref n="1" />. It is a statistical value that represents the proportion of a dependent variable variance that is predictable from the independent variables:

$$R^2 = 1 - \frac{\sum_{i=1}^{n}(Y_{\text{pred},i} - Y_{\text{real},i})^2}{\sum_{i=1}^{n}(\bar{Y}_{\text{pred}} - Y_{\text{real},i})^2}, \quad (8)$$

where the notations are the same as the RMSE ones and $\bar{Y}_{\text{pred}}$ is the mean of the predictions. In a regression problem, the coefficient of determination measures how well the model fits the real data. It has the advantage of being confined by construction between 0 and 1, 0 being the worst fit possible and 1 meaning that the model perfectly fits the data. By studying the RMSE and the $R^2$ we can have a global evaluation of the capacity of a model to make the correct prediction.

To gain more information about the in-depth behavior of a model, we also use the relative errors (RE):

$$RE_i = \frac{Y_{\text{pred},i} - Y_{\text{real},i}}{Y_{\text{real},i}}. \quad (9)$$

<footnote n="1">¹ https://search.nal.usda.gov/permalink/01NAL_INST/178fopj/alma9916341892807426</footnote>
Studying the distribution of the RE of a model when predicting the test set can reveal any systematic bias during the predictions. Furthermore, since our datasets are made of simulations, we can check if any of the simulation parameters primarily influence the results of the predictions by looking at the relations between the parameter values and the RE.

# 5. Results on $M_{BH}$ estimation

## 5.1. Reliability of ML models

As described above (Sects. 2, 7, of the 10), the galaxies in our sample have sufficiently robust CO(3–2) data in the central regions to enable $M_{BH}$ estimations with our ML method: NGC 4388, NGC 5506, NGC 5643, NGC 6300, NGC 7314, NGC 7465, and NGC 7582. Before comparing our $M_{BH}$ estimations with those from the literature, we need to verify their reliability with the tools described in Sect. 4.5.

We report in Table 4 the measurement of the RMSE, the $R^2$, and the mean RE of our models on their test set for both the ones estimating the $M_{BH}$ via the PVD and the ones using the first-moment maps. The analysis reveals that the diverse evaluations of the fitting quality exhibit a high degree of reliability and demonstrate consistency across the various indicators for both architectures. The RMSE values are at worst 0.0367 and the $R^2$ are all above 0.946, confirming that all the models are capable of making a reliable estimation of the $M_{BH}$ on simulated data that were not present during the training process.

In Fig. 2, we show the models’ RE distributions estimations on the 5000 simulations that composed their respective test sets. As shown in Table 4, the mean RE is ~0 for the fourteen models, two for each galaxy, one with the PVD as input and one with the moment maps as input. The histograms are centered at the value of zero and prominently feature a pronounced peak at this center point. The locations of the centers around 0 show that there is no systematic bias during the predictions and the relatively symmetrical distributions indicate that there is no constant over- and underestimation of the $M_{BH}$.

The mock datasets used for training, validation, and testing were produced via a numerical simulation with ten free parameters. We can use the RE to make sure that the confidence of our models $M_{BH}$ predictions are not influenced by any other parameters than the $M_{BH}$. For each parameter, we plotted the RE of the $M_{BH}$ predictions as a function of the parameter values. We observed the same behavior for each galaxy so we present in the Appendix only the case of NGC 6300 (i.e., see Appendix B). For every parameter, except for the $M_{BH}$, the distribution is seen to be uniform and centered around zero. This is strong evidence for considering that the models prediction does not depend on any other parameter other than $M_{BH}$.

To further investigate the behavior of our models, we represent the 2D histograms of the RE as a function of $M_{BH}$ (i.e., see Appendix C). All the distributions have the same shape. From $log(M_{BH}/M_\odot)$ of 4.0 to ~4.5 dex, the distribution means increase linearly from <0 to ~0 and their scatter stays constant. At ~4.5 dex, the scatters get bigger and we see outliers that are underestimated. Then, the distribution scatters get thinner and the outliers number decreases as the $M_{BH}$ predicted increases until ~6.5 dex or ~7 dex and stays constant afterward. We note the exception of the NGC 5643 model using the moment maps. Unlike the other models, its mean is already at 0 for $log(M_{BH}/M_\odot) = 4$ dex and stays constant with $M_{BH}$. Its scatter is also constant and uniform. As the $M_{BH}$ gets lower, fewer and fewer pixels contain information about it. Thus, it becomes

harder for the NN to extract meaningful information out of the data. This may explain why we observe the high scatter below $log(M_{BH}/M_\odot) \sim 6$ dex.

## 5.2. $M_{BH}$ estimation results

The statistical indicators showed that the models succeeded in predicting precisely the $M_{BH}$ on unseen simulated data when the $log(M_{BH}/M_\odot)$ are above ~6 dex. Thus, if the numerical simulation can accurately emulate the observations, our models should be able to have the same behavior when estimating the $M_{BH}$ of their corresponding observations. In Fig. 3, we show our model estimations on their respective observations. We made $10^5$ predictions per galaxy. The observed scatter on the $M_{BH}$ predictions comes from the MC dropout. These scatters represent the models confidence in their predictions and do not have a predefined shape. The observed Gaussian shapes come from the fact that our models are confident in their estimations. Moreover, it is important to note that the models are trained on numerical simulations that do not encompass all the underlying physics observed in reality. Thus, the incertitude shown here cannot be associated with physical ones but only with the models confidence in their estimations.

We can observe that the distributions are all Gaussian, except for the NGC 5643 model using the moment maps. This particular model seems to be the sum of two Gaussians with different dispersions, but centered around the same value (6.99 dex). The width of the Gaussians given by the PVD models is broader compared to the moment maps models. Furthermore, the models scatter is constant for the moment map models, whereas it is galaxy-dependent for the PVD models. This can be explained if the PVD shows only a part of the gas dynamic. The effectiveness of the NN will rely on the extent to which the presented information is valuable for estimating the $M_{BH}$. On the contrary, the moment maps show the entirety of the gas dynamics, so helpful information is always represented. On one hand, the moment maps contain information primarily concentrated in the values of the central pixels, providing a localized representation of the $M_{BH}$-related information. On the other hand, the PVD captures information through the velocity distribution as a function of position. Particularly, in the shape of this distribution in the vicinity of the center. This results in a more sparse distribution of information which makes it more difficult for the network to relate this information to the $M_{BH}$. Another explanation is that we use different architectures for the PVD and the moment maps.

We can observe (in Fig. 4) that for four out of the seven galaxies, the estimations of the PVD and the moment maps do not overlap. This can be attributed to the same hypotheses as to why the PVD distributions are broader than the moment maps.

We compare our $M_{BH}$ estimations with the ones previously made using classical methods in Table 5 and Fig. 4. For the seven galaxies in our sample, the $M_{BH}$ predicted by our models is coherent with previous estimations. There is some overlap with at least one previous estimation, except for NGC 7314 and NGC 7582. For these, our estimations fall between the value from Wang & Zhang (2007) and the three reported in Bian & Gu (2007), van den Bosch (2016), She et al. (2017).

## 5.3. Comparison with classical method

In this study, we used 25 000 moment maps of NGC 7582 previously simulated datasets, originally generated for ML methods, to explore traditional non-ML approaches for estimating $M_{BH}$, specifically linear interpolation and chi-squared fitting. These non-ML methods provided a useful comparison to ML
<table><caption>**Table 4.** ML models statistics.</caption><thead><tr><th rowspan="2">ID</th><th colspan="2">RMSE</th><th colspan="2">R<sup>2</sup></th><th colspan="2">RE</th></tr><tr><th>PVD</th><th>Mom</th><th>PVD</th><th>Mom</th><th>PVD</th><th>Mom</th></tr></thead><tbody><tr><td>NGC 4388</td><td>0.0341</td><td>0.0341</td><td>0.986</td><td>0.946</td><td>0.0020</td><td>0.0007</td></tr><tr><td>NGC 5506</td><td>0.0335</td><td>0.0367</td><td>0.986</td><td>0.947</td><td>0.0060</td><td>-0.0032</td></tr><tr><td>NGC 5643</td><td>0.0194</td><td>0.0178</td><td>0.995</td><td>0.996</td><td>-0.0019</td><td>0.0001</td></tr><tr><td>NGC 6300</td><td>0.0243</td><td>0.0180</td><td>0.995</td><td>0.996</td><td>-0.0002</td><td>0.0005</td></tr><tr><td>NGC 7314</td><td>0.0312</td><td>0.0277</td><td>0.984</td><td>0.990</td><td>-0.0029</td><td>-0.0041</td></tr><tr><td>NGC 7465</td><td>0.0246</td><td>0.0255</td><td>0.993</td><td>0.992</td><td>0.0010</td><td>0.0014</td></tr><tr><td>NGC 7582</td><td>0.0263</td><td>0.0220</td><td>0.992</td><td>0.994</td><td>-0.0002</td><td>-0.0031</td></tr></tbody></table>

**Notes.** RMSE, R² and mean RE of our models on the test sets. The RMSE (>0), R² (between 0 and 1) and the RE are defined in Section 4.5.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Histograms of the RE of our models when estimating the $M_{BH}$ of the 5000 simulations from their respective test set. We represent in blue the models trained to estimate the $M_{BH}$ with the PVD, and in red the models trained with the moment maps.</figcaption></figure>

models, as they rely on conventional estimation techniques rather than data-driven learning. Linear interpolation estimates unknown values by assuming a linear relationship between nearby data points, while chi-squared fitting minimizes the discrepancies between observed and predicted data. Both methods were applied to these simulations to evaluate their performance in estimating $M_{BH}$.

Linear interpolation works by assuming that the changes between adjacent data points are linear, allowing unknown values within the dataset to be estimated. In our simulations, the data was split into 20 000 training points and 5000 test points. However, due to the high dimensionality of the data, memory limitations arose when attempting to interpolate across the entire dataset. To overcome this, we adopted a localized interpolation strategy. Instead of interpolating over the full dataset, we identified the K nearest neighbors in multidimensional space (those points closest to the test data) and performed linear interpolation on this subset. This localized approach reduced both memory usage and computation time. Despite its efficiency, the linear

interpolation method was less accurate than the ML models. As shown in Fig. 5, the histogram of RE from linear interpolation indicates a peak near zero, but the errors were generally higher than those observed with ML estimations. Moreover, the error distribution was asymmetric, implying that the linear assumption in high-dimensional space may not fully capture the data's complexity. The 2D histograms of RE provide additional insight into the limitations of linear interpolation (Fig. C.7). Unlike the ML model, which consistently estimated $M_{BH}$ starting from around 5.5 dex, the linear interpolation method only showed reliable performance at the upper end of the $M_{BH}$ distribution, approximately from 8 dex onwards. This suggests that linear interpolation struggles to produce accurate estimates across the entire $M_{BH}$ range, particularly in the lower regions where relationships between variables are likely more complex and non-linear. As a result, while linear interpolation may offer computational efficiency, its ability to generalize across different regions of the $M_{BH}$ distribution is limited, making it less reliable than ML approaches for comprehensive mass estimation.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Probability distributions of the $M_{BH}$ predicted by our models for the seven galaxies in our sample. We have plotted the estimations made with the PVD as input in blue, and the estimations made with the first moment maps as inputs in red.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of the probability distributions of the $M_{BH}$ predicted by our models and the ones found in the literature. We represent in blue the results of the estimation made from the PVD and in red the estimation made from the moment maps. The other Gaussians represent the estimations found in the literature. The y-axis upper limit is fixed at 2 arbitrarily fixed at 2 for clarity purposes.</figcaption></figure>

Chi-squared fitting in 2D works by comparing observed data points to a theoretical model and minimizing the differences between them. For each pair of observed and expected values, the squared difference is calculated and normalized by the variance of the observed data. These individual terms are then

summed to produce the chi-squared statistic, which quantifies the overall goodness of fit. The method iterates through possible model parameters, adjusting them to reduce the chi-squared value. Figure 5 presents 2D histograms illustrating the distribution of chi-squared results, which exhibit a similar shape to those
<table><caption>**Table 5.** Summary of $M_{BH}$ estimations.</caption><thead><tr><th>ID</th><th>logM<sub>BH</sub> PVD</th><th>logM<sub>BH</sub> Mom</th><th>M-σ</th><th>Reverberation</th><th>FPBH</th><th>Ref</th></tr></thead><tbody><tr><td>NGC 4388</td><td>7.18 ± 0.16</td><td>6.61 ± 0.02</td><td>7.17<sup>8.16</sup><sub>6.85</sub></td><td>—</td><td>7.62</td><td>3,4,5,6,7,8</td></tr><tr><td>NGC 4941</td><td>—</td><td>—</td><td>6.83<sup>7.14</sup><sub>6.83</sub></td><td>6.34</td><td>—</td><td>2,3,6,7</td></tr><tr><td>NGC 5506</td><td>6.40 ± 0.14</td><td>6.63 ± 0.02</td><td>7.61<sup>8.24</sup><sub>6.65</sub></td><td>6.30</td><td>8.15</td><td>1,3,4,6,7,8,9</td></tr><tr><td>NGC 5643</td><td>6.13 ± 0.10</td><td>6.99 ± 0.02</td><td>6.73<sup>7.05</sup><sub>6.44</sub></td><td>6.38<sup>6.45</sup><sub>6.30</sub></td><td>—</td><td>2,3,5,6,7</td></tr><tr><td>NGC 6300</td><td>7.15 ± 0.12</td><td>6.39 ± 0.01</td><td>6.63<sup>7.14</sup><sub>5.45</sub></td><td>6.29</td><td>—</td><td>3,4,5,6,7</td></tr><tr><td>NGC 6814</td><td>—</td><td>—</td><td>6.96<sup>7.32</sup><sub>6.42</sub></td><td>7.28</td><td>7.45</td><td>1,4,5,7,9</td></tr><tr><td>NGC 7213</td><td>—</td><td>—</td><td>8.49<sup>8.82</sup><sub>8.06</sub></td><td>7.35<sup>7.99</sup><sub>6.88</sub></td><td>—</td><td>1,3,7</td></tr><tr><td>NGC 7314</td><td>6.84 ± 0.05</td><td>6.41 ± 0.02</td><td>6.48<sup>7.83</sup><sub>5.59</sub></td><td>6.70</td><td>7.18</td><td>1,4,6,7,9</td></tr><tr><td>NGC 7465</td><td>6.75 ± 0.05</td><td>6.62 ± 0.02</td><td>6.79<sup>7.07</sup><sub>6.51</sub></td><td>—</td><td>—</td><td>6,7</td></tr><tr><td>NGC 7582</td><td>6.62 ± 0.03</td><td>6.40 ± 0.01</td><td>7.56<sup>7.76</sup><sub>7.24</sub></td><td>5.59</td><td>—</td><td>2,3,5,6,7</td></tr></tbody></table>

**Notes.** Table summarizing the different $M_{BH}$ estimations found in the literature and the one we made with ML. The $M_{BH}$ is in $M_{\odot}$, and only the log of values are indicated. When multiple $M_{BH}$ estimations using the same method were available, we report here their mean, the highest values in superscript, and the lowest in the subscript. (1) Galaxies ID; (2) Mean of the $M_{BH}$ probability distribution obtained with ML and the PVD as input (the errors are at $1\sigma$); (3) Same as the previous column but with the moment maps as input; (4) Mean of the $M_{BH}$ estimations with the $M-\sigma$ (intrinsic errors of 0.3 dex); (5) $M_{BH}$ calculated with reverberation mapping (intrinsic errors of 0.4 dex); (6) $M_{BH}$ estimations references 1: Peterson et al. (2004), 2: Bian & Gu (2007), 3: Wang et al. (2007), 4: Beckmann et al. (2009), 5: van den Bosch (2016), 6: She et al. (2017), 7: $M-\sigma$ estimation made with the Hyperleda data, 8: FPBH mass estimated from the XMM-Newton and the FIRST data, 9: FPBH mass estimated from the XMM-Newton and the NVSS data).

obtained using linear interpolation. Both methods display a comparable constant scatter, not centered around zero, indicating a shared unreliability in estimating $M_{BH}$ below approximately 8 dex. Notably, the scatter in the lower mass range is more pronounced in the chi-squared 2D histograms. However, the 1D histograms reveal that the RE distribution for the chi-squared method is more symmetric, suggesting improved consistency in error distribution when compared to linear interpolation.

In conclusion, both the linear interpolation and chi-squared fitting methods exhibit significant challenges in reliably estimating $M_{BH}$ within the 4 to 8 dex range. In contrast, our ML approach demonstrates improved performance, effectively mitigating this limitation and confining the estimation inaccuracies to a narrower range of 4 to 6 dex. This highlights the enhanced robustness and precision of the ML method, particularly for mass estimates in the lower dex intervals.

# 6. Individual galaxies

The three moments of the CO(3–2) cubes from ALMA observations (top panel) are compared with the best fit of our simulations (bottom panel) on the left sides of Figs. A.1 to A.7. On the right, the observed PVD is overlaid with the contours of the second model, optimized for this. The models are fitted with the same methodology as in Combes et al. (2019).

## 6.1. NGC 4388

The galaxy is almost edge-on (IA of 79°), and its nuclear molecular disk reveals a ring of ~100 pc radius, which corresponds to the inner Linblad resonance of the bar. The ring is not totally empty and shows clumps that enter the SOI, allowing a $M_{BH}$ estimation. Some molecular gas is expelled along the minor axis and might be dragged by the ionized wind and radio jet (Falcke et al. 1998).

Figure A.1 shows the model best fit found for NGC 4388. It is clear in this figure that the velocity field of the galaxy is highly asymmetric, and that our symmetric gravitational potential model cannot reproduce the observed first-moment map.

While the gas with negative velocity distribution is well reproduced, the positive part is not. On the PVD panel, the model fits better than the symmetrized observation. The observed second moment shows a sudden drop of velocity dispersion at $P \sim 200$ pc that is not present in the model. Moreover, the model has a large velocity gradient at the center contrary to the observations. This comes from the lack of gas in the SOI, so it is difficult to precisely measure the SMBH mass in this galaxy. It can be noted that the lowest red contour in the PVD is less reliable, related to a low level of emission.

## 6.2. NGC 5506

NGC 5506 is also almost edge-on (IA of 80°), and possesses a nuclear molecular ring of radius ~50 pc, corresponding to the ILR of its strong bar. There is however dense molecular gas near the AGN, and HCO$^+$ has been detected in the center (García-Burillo et al. 2021). The estimation of $M_{BH}$ is therefore possible.

Figure A.2, reveals in the center of the observed first-moment map that the gas velocity distribution has an S-like shape. This is probably due to the presence of a bar. Except for these details, the model can well reproduce the observed gas velocity distribution. However, it cannot do as well for the velocity dispersion distribution. The model gives a more or less constant velocity dispersion whereas two clumps with high-velocity dispersion are observed near the center and a low-velocity dispersion elsewhere. This effect is also seen in the PVD and could be also due to the bar. The low gas density in the center makes the SMBH mass estimation uncertain. This is reflected in the low-emission red contour of the model, which proposes in general a higher $M_{BH}$ than the moment-based model.

## 6.3. NGC 5643

This is an almost face-on barred spiral, with an IA of 29°. The strong stellar bar is almost E-W, with $PA = 85°$. At ~100 pc
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** 2D and 1D histograms of the RE from linear interpolation and chi squared fitting of $M_{BH}$ on the test set for NGC 7582. On the left, we plotted the RE when estimating the $M_{BH}$ of the test set as 2D histograms. On the right, we plotted the same RE but as a 1D histogram. On top, we plotted the RE from the linear interpolation and at the bottom, we plotted the RE from the chi squared fitting. The y-axis represents the RE values. The x-axis of the 2D histograms is in $\log(M_{BH}/M_{\odot})$ and the 1D histogram axis is in number of estimations.</figcaption></figure>

scale, NGC 5643 CO(3–2) moment zero map shows the presence of a nuclear bar, prolonged by a spiral structure, but the CO bar is almost N-S, with $PA = 5°$, cf. Fig. A.3. The very center reveals a tilted molecular torus, of diameter 26 pc, oriented almost N-S, perpendicular to the ionized gas outflow and radio jet (Alonso-Herrero et al. 2018; García-Bernete et al. 2021). Dense gas traced by HCO$^+$ has been detected in the very center (García-Burillo et al. 2021).

The NGC 5643 model reproduces well the shift in the velocity orientation towards the nucleus (Fig. A.3). However, the second-moment map reveals a high-velocity dispersion in the nucleus, which is horizontal in the model while almost vertical in the observation. In the PVD, the model does not distribute the gas clumps at the same positions as in the observed diagram; however, the velocity gradient in the nucleus is reproduced, and the SMBH mass estimation appears robust. Note that the moment-based model leads to a higher $M_{BH}$ in this galaxy, which is unusual, maybe due to the strong misalignment.

## 6.4. NGC 6300

This is a strongly barred galaxy, moderately inclined (57°). At ~100 pc scale, the CO(3–2) emission shows a detached nuclear disk, with a central concentration, allowing a determination of $M_{BH}$, cf. Fig. A.4. This strong concentration corresponds to a molecular torus of 25 pc radius, also detected in HCO$^+$. The torus oriented at $PA = 85°$, is perpendicular to the ionized gas outflow, detected by Scharwächter et al. (2016).

For the NGC 6300 fit (Fig.A.4), the first-moment map is well reproduced by the model, except for the wiggle at the center of the disk, likely due to the bar. The modeled second-moment map underestimates the velocity dispersion at the center but overestimates it near the edges. Furthermore, the central region with high-velocity dispersion is horizontal in the observation while it is almost vertical in the model, a likely consequence of the bar. The PVD panel reveals a reasonable fit, with a sufficient gas sampling of the SoI.

## 6.5. NGC 7314

This barred spiral galaxy looks relatively inclined (70°) on the sky, but is kinematically less (55°) since its morphological shape is biased by the bar. At ~100 pc scale, the molecular disk reveals a central depletion, surrounded by two peaks, suggesting the existence of a ring of radius 20 pc, that could be an ILR resonant ring. Its PVD shows very little gas in the SoI but it is enough to estimate an upper limit of 6.2 dex on the $\log(M_{BH}/M_{\odot})$. In this galaxy, the molecular nuclear disk does not appear strongly misaligned with the large-scale one. It is also oriented along the radio source (Thean et al. 2000), which might be dominated by star formation.

The NGC 7314 fit (Fig. A.5) overestimates the velocity gradient in the outer parts of the disk. The model also overestimates the central velocity dispersion along the minor axis, while the dispersion is more spatially distributed in the observation. The PVD reveals a high central velocity dispersion, leading to a high estimation of the SMBH mass. This is however provided by the lowest emission red contour.

## 6.6. NGC 7465

This barred early-type spiral is moderately inclined (60°) on the sky and its morphological major axis of $PA = 160°$ at large-scale is roughly perpendicular to the molecular disk at ~100 pc scale, with $PA = 45°$. The circum-nuclear disk in the regular rotation is dense, of 70 pc radius, and prolonged by a sparse spiral structure. There is enough dense gas towards the center to estimate $M_{BH}$ through the ML method. The inferred $\log(M_{BH}/M_{\odot})$ is 5.9 ± 0.8 dex.

The model for NGC 7465 (Fig.A.6) reproduces very well the observed first-moment map. The velocity dispersion is however gathered along the minor axis, while it is more spatially distributed in the observed second moment. The low-velocity gradient in the SoI region leads to a low mass estimation for the SMBH.
## 6.7. NGC 7582

NGC 7582 is a barred early-type spiral, highly inclined (68°), with a ~100 pc scale molecular disk aligned with the large-scale galaxy disk, with $PA = 156°$. The molecular gas is condensed in a 200 pc radius ring, which corresponds to an ILR ring. However, there is enough CO emission right into the very center, making it possible to estimate $M_{BH}$. There is even a concentration of CO emission in a small disk of radius 30 pc, corresponding to a molecular torus. An outflow of ionized gas is observed in the perpendicular direction (Davies et al. 2016).

The NGC 7582 model fits the observed velocity gradient in the disk relatively well, although its morphology is too regular (Fig. A.7). The observed first-moment map shows that the isovelocity morphology gets thinner towards the nucleus, with a large central velocity gradient, which might be due to the bar. The velocity dispersion reflects also this morphology and is not reproduced in the model. The PVD reveals a low-velocity gradient in the SoI and a low SMBH mass estimation.

## 7. Discussion

To provide more robust estimations of $M_{BH}$, there are three key areas of improvement: input data, numerical simulations, and the ML method itself. First of all, the model precision is limited by the observation's spatial resolution. The 2D histograms displayed in Appendix C demonstrate that NN face challenges in providing accurate estimations when the $log(M_{BH}/M_⊙)$ is below approximately 6 dex. This issue arises because, as $M_{BH}$ decreases, fewer pixels in the data carry information about it. Consequently, the NN struggles to extract meaningful insights, leading to increased variability in estimations for smaller $M_{BH}$ values. Therefore, observations with higher resolutions provide more pixels with relevant $M_{BH}$ information, enhancing the NN precision and reducing the sensitivity to small $M_{BH}$ variations. The role of the SoI in the context of NN remains somewhat unclear. While the goal is to increase the number of informative pixels about the BH, these pixels do not need to be situated where the BH dominantly influences gas dynamics. Rather, these pixels need only exhibit a noticeable effect on the NN, even if it is relatively small.

Another limiting factor comes from the physics included in the numerical simulations used to make the datasets. We used numerical simulations to mimic real observations, but we simulated an axisymmetric and homogeneous thin gas disk, with no clumpy structure, nor any spiral or bar structures. Furthermore, in our models, we adopted the same IA and PA for the gas and stellar disks, while some misalignments of the nuclear molecular disks have been observed with ALMA (Combes et al. 2019). Adding these misalignments and more complex physics, such as the presence of a bar, in the simulations would add free parameters to influence directly the $M_{BH}$ estimations. This discrepancy in realism can result in the network exhibiting excessive confidence in its predictions, while there is a strong likelihood that the uncertainty in the actual data is being underestimated.

Additionally, as previously elaborated in Section 5.2, the probability distributions illustrated in Fig. 3 lack physical significance. This discrepancy arises from the inherent differences in the underlying physics between the numerical simulations and the actual observational data. However, by applying a reverse perspective, as the physical fidelity of the simulations approaches that of the observations, the probability distributions of $M_{BH}$ become more representative of the actual physical uncertainties. Therefore, enhancing the realism of the numerical sim-

ulations serves a dual purpose, it not only improves the accuracy of error estimates but also confers them with a deeper and more grounded physical significance.

The ML method itself can also be improved. Here, we chose to work with one galaxy at a time which means creating a different training sample and training a different model for each galaxy. This approach has the advantage of simplifying the training set creation and the learning process but it is also very inefficient time-wise and this approach is effective solely due to the limited number of galaxies available in our dataset. It takes ~8 hours to create a dataset and ~4 hours to train one NN.

To enhance the efficiency and applicability of our method, it is essential to generalize it to work across multiple galaxies rather than being constrained to individual ones. By creating a training sample that includes moment maps from multiple different galaxies, we can develop a NN capable of generalizing and working for any galaxy. Generalization requires a robust training dataset that includes moment maps and PV diagrams from a wide variety of galaxies. This diversity is crucial for the NN to learn the different galactic structures, dynamics, and environments it might encounter.

Moreover, when working with multiple galaxies simultaneously, we expect the network to be able to mutualize or share learned features across different galaxies. This allows the NN to extract common patterns and relationships between galactic properties, enhancing its ability to generalize. Consequently, while the architecture needed for such a task might be more complex due to the need to handle diverse galactic environments, it does not necessarily have to be much larger in terms of parameters. The mutualization of features could lead to a more efficient learning process, allowing the network to converge without significantly increasing its size. Therefore, the inference time should not see a substantial increase, ensuring that the method remains computationally feasible even as it scales to a broader range of galaxies.

However, our current limited sample size allows us to center the moment maps and PV diagrams around the BH, a process that becomes challenging with a larger number of galaxies. To address this, incorporating a BH localization step within the network can further improve performance. This step involves identifying the position of the BH in the moment maps and PV diagrams before estimating its mass, enabling the model to handle cases where the BH is not centrally located.

To effectively learn and predict additional physical parameters and/or accurately detect the position of BH, it is necessary to make modifications to the NN architecture to address these challenges. Moreover, a larger and more diverse training sample will be required. When aiming to predict multiple parameters, it is essential to ensure that the training sample is sufficiently extensive to comprehensively map the entire feature space. Furthermore, to successfully learn how to localize the BH in the input data, the positional information of the BH must be known for each example. However, the specific size and diversity requirements of the sample are empirical and contingent upon the particularities of the problem at hand. At the very least, several hundred distinct galaxy molecular gas observations are anticipated to be necessary.

Then, we can choose to train our NN on the observations or simulations made from the observations. Training the NN on the observations has the advantage that we do not have problems with the realism of the physics used in the numerical simulations. However, it requires us to have $M_{BH}$ estimation for all the observed galaxies. In addition, it is important to note that a single galaxy may yield multiple $M_{BH}$ estimations derived from various methods, often resulting in discrepancies between these
estimates. This situation highlights the necessity for method-
ological consistency in estimating M<sub>BH</sub> across all galaxies in the
sample to avoid potential confusion for the NN. However, this
uniformity in estimation methods also implies that the NN is not
learning to estimate M<sub>BH</sub> per se but is instead learning to esti-
mate M<sub>BH</sub> based on the predetermined method. Consequently, if
the selected method contains inherent flaws or inaccuracies, the
NN is likely to incorporate and replicate these shortcomings in
its estimations. Working with numerical simulations does have
its disadvantages, particularly concerning the intricacies of simu-
lating the physics, as elaborated earlier. However, it also presents
the advantage of reducing the volume of observations necessary
for the method generalization, as it allows for the creation of
multiple simulations per galaxy.

Nevertheless, obtaining the exact position of the BH may not
be a prerequisite; instead, it may suffice to ensure that the BH is
contained within the image. This approach offers the advantage
of simplifying the acquisition of the required sample. However,
it is important to note that training the model using this method
may be harder.

8. Conclusion

The high sensitivity and spatial resolution of ALMA allow us
to resolve the CO(3-2) molecular disks at the center of nearby
galaxies. In the case of the GATOS core sample used in this
work, the distances and ALMA angular resolutions allow for the
resolution on typical physical scales of 10 pc, which probe their
SoI. Since the nuclear disks entered the sphere of influence of
the BH, we were able to make M<sub>BH</sub> estimations. In this article,
we provide a proof of concept for a novel approach to estimate
the mass of SMBH using ALMA observations of the circum-
nuclear disk. Our method employs supervised ML techniques on
data obtained from numerical simulations.

Gas dynamics can be studied by looking at the PVD and/or
the moments map calculated from a spectral data cube. we devel-
oped two artificial NN architectures. The first one takes the PVD
as input, the other one takes the first moment maps, and both
of them give an estimation of the M<sub>BH</sub> based on the input. To
test the reliability of our models, we used three statistical indica-
tors: the RMSE, the R<sup>2</sup>, and the <span style="text-decoration: overline">RÉ</span>. These statistical indicators
confirmed that our fourteen models can make reliable M<sub>BH</sub> esti-
mations on new data. we looked at our models RE distribution
when predicting the M<sub>BH</sub> of the test sample. we saw that the dis-
tributions were symmetrical and centered around 0, indicating no
constant over- or underestimation. We also studied the models'
RE of their M<sub>BH</sub> estimations made on the test set as a function of
the simulation parameters values. The resulting distributions are
all uniform and center around zero except for the M<sub>BH</sub>. This is a
strong indication that the models M<sub>BH</sub> predictions only depend
on the M<sub>BH</sub>. We studied the behavior of the M<sub>BH</sub> predictions in
depth by looking at the predictions RE versus the M<sub>BH</sub>.

We managed to make a M<sub>BH</sub> estimation with both architec-
tures for seven out of the ten galaxies in the GATOS core sam-
ple: NGC 4388, NGC 5506, NGC 5643, NGC 6300, NGC 7314,
NGC 7465, and NGC 7582. We made 10<sup>5</sup> predictions per
galaxy using MC dropout. The resulting M<sub>BH</sub> estimations are
Gaussian-like probability distributions. The scatter of these dis-
tributions represents the models confidence in their predictions.
We observe that the PVD models give larger scatters from
0.03 dex to 0.16 dex, depending on the galaxy. On the other hand,
the moment maps models have a constant scatter of ~0.02 dex
for all galaxies. The narrow error bars observed in our results
are likely indicative of the network's overconfidence, which

can be attributed to its training on numerical simulations that
encompass simplified physics compared to real-world condi-
tions. Nonetheless, all the estimations made are consistent with
the previous studies except for NGC 7314 and NGC 7582 where
our estimate is situated between the previous ones.

In this study, we show that even with training sets simulated
with simple physics, along with a fixed IA and PA and simple
NN architectures, our approach can produce results that are fully
coherent with results reported in the literature. This work repre-
sents the first step toward an automatized method for estimating
MBH.

In this study, we demonstrate that even with training sets sim-
ulated using simple physics, fixed IA and PA, and basic NN archi-
tectures, our approach can produce results consistent with those
in the literature. However, a key limitation of this work is that we
trained and applied the NNs on one galaxy at a time. As a result, we
had to generate numerical simulations and train a separate model
for each individual galaxy. This approach limits scalability and
the ability to generalize across different galactic environments.
While our method is not necessarily superior to classical methods,
it serves as a proof of concept and an initial step towards develop-
ing ML techniques for M<sub>BH</sub> estimations that can improve scaling
with the increasing number of astronomical observations.

Data availability

Appendices B and C, which present the RE of M<sub>BH</sub> predic-
tions as a function of parameters for NGC 6300 and the 2D
histograms of the models RE, respectively, are available at
https://zenodo.org/records/14244314.

Acknowledgements. We would like to thank the referee for very useful com-
ments, that helped us improve and clarify the paper. RP acknowledges finan-
cial support from the SNSF under the Weave/Lead Agency RadioClusters grant
(214815). AAH and MVM acknowledge support from PID2021-124665NB-I00
by the Spanish Ministry of Science and Innovation/State Agency of Research
MCIN/AEI/ 10.13039/501100011033 and by 'ERDF A way of making Europe.
SGB acknowledges support from the research project PID2019-106027GA-C44
of the Spanish Ministerio de Ciencia e Innovación. Training computation and
inferences have been performed on the MINERVA group from the Observatoire
de Paris GPUs server. AJB has received funding from the European Research
Council (ERC) under the European Union's Horizon 2020 Advanced Grant
789056 "First Galaxies". Based on observations made with the NASA/ESA Hub-
ble Space Telescope, and obtained from the Hubble Legacy Archive, which is
a collaboration between the Space Telescope Science Institute (STScI/NASA)
the Space Telescope European Coordinating Facility (ST-ECF/ESA), and the
Canadian Astronomy Data Centre (CADC/NRC/CSA). MPS acknowledges fund-
ing support from the Ramón y Cajal program of the Spanish Ministerio de
Ciencia e Innovación (RYC2021-033094-I). C.R. acknowledges support from
the Fondecyt Regular grant 1230345 and ANID BASAL project FB210003.
OGM acknowledges support from UNAM PAPIIT project IN109123 and CONA-
CyT "Frontera de la Ciencia" project CF-2023-G-100. BGL acknowledges sup-
port from grants PID2019-107010GB-100 and the Severo Ochoa CEX2019-
000920-S funded by MICINN-AEI/10.13039/501100011033. CRA acknowl-
edges financial support from the European Union's Horizon 2020 research
and innovation program under Marie Skłodowska-Curie grant agreement No
860744 (BiD4BESt), from the State Research Agency (AEI-MCINN) and from
the Spanish MCINN under grants "Feeding and feedback in active galaxies",
with reference PID2019-106027GB-C42, the project "Quantifying the impact of
quasar feedback on galaxy evolution", with reference EUR2020-112266, funded
by MICIN/AEI/10.13039/501100011033 and the European Union NextGenera-
tionEU/PRTR.

References

Alonso-Herrero, A., Pereira-Santaella, M., García-Burillo, S., et al. 2018, ApJ,
859, 144

Alonso-Herrero, A., García-Burillo, S., Pereira-Santaella, M., et al. 2019, A&A,
628, A65

Audibert, A., Combes, F., García-Burillo, S., et al. 2019, A&A, 632, A33
Bahcall, J. N., Kozlovsky, B.-Z., & Salpeter, E. E. 1972, ApJ, 171, 467

Beckmann, V., Soldi, S., Ricci, C., et al. 2009, A&A, 505, 417

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Bender, R., Kormendy, J., Bower, G., et al. 2005, ApJ, 631, 280

Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2009, ApJ, 705, 199

Bian, W., & Gu, Q. 2007, ApJ, 657, 159

Blandford, R. D. 1984, Ann. New York Acad. Sci., 422, 303

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Castro-Ginard, A., Jordi, C., Luri, X., Cantat-Gaudin, T., & Balaguer-Núñez, L. 2019, A&A, 627, A35

Cherepashchuk, A. M., & Lyutyi, V. M. 1973, ApLett, 13, 165

Chollet, F., et al. 2015, Keras, https://keras.io

Combes, F., García-Burillo, S., Audibert, A., et al. 2019, A&A, 623, A79

Costa, T., Sijacki, D., & Haehnelt, M. G. 2014, MNRAS, 444, 2355

Davies, R. L., Dopita, M. A., Kewley, L., et al. 2016, ApJ, 824, 50

Davis, T. A., Bureau, M., Onishi, K., et al. 2017, MNRAS, 468, 4675

Davis, T. A., Bureau, M., Onishi, K., et al. 2018, MNRAS, 473, 3818

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451

Fabian, A. C. 2012, ARA&A, 50, 455

Falcke, H., Wilson, A. S., & Simpson, C. 1998, ApJ, 502, 199

Falcke, H., Koerding, E., & Markoff, S. 2004, A&A, 414, 895

Ferrarese, L., & Ford, H. 2005, Space Sci. Rev., 116, 523

Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9

Gallo, E., Fender, R. P., & Pooley, G. G. 2003, MNRAS, 344, 60

García-Bernete, I., Alonso-Herrero, A., García-Burillo, S., et al. 2021, A&A, 645, A21

García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2019, A&A, 632, A61

García-Burillo, S., Alonso-Herrero, A., Ramos Almeida, C., et al. 2021, A&A, 652, A98

Gaskell, C. M. 1988, ApJ, 325, 114

Gaskell, C. M. 2008, Rev. Mex. Astron. Astrofis. Conf. Ser., 32, 1

Gaskell, C. M. 2011, Balt. Astron., 20, 392

Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13

Gebhardt, K., Richstone, D., Tremaine, S., et al. 2003, ApJ, 583, 92

Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580

Gültekin, K., King, A. L., Cackett, E. M., et al. 2019, ApJ, 871, 80

Hartley, P., Bonaldi, A., Braun, R., et al. 2023, MNRAS, 523, 1967

Huertas-Company, M., Guo, Y., Ginzburg, O., et al. 2020, MNRAS, 499, 814

Kaspi, S., Brandt, W. N., Maoz, D., et al. 2007, ApJ, 659, 997

Kingma, D. P., & Ba, J. 2014, ArXiv e-prints [arXiv:1412.6980]

Kollatschny, W., & Zetzl, M. 2011, Nature, 470, 366

Kormendy, J. 2004, Coevolution of Black Holes and Galaxies, 1

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581

Lanusse, F., Ma, Q., Li, N., et al. 2018, MNRAS, 473, 3895

LeCun, Y., Bengio, Y., & Hinton, G. 2015, Nature, 521, 436

Liu, H. T., Feng, H. C., & Bai, J. M. 2017, MNRAS, 466, 3323

Lochner, M., McEwen, J. D., Peiris, H. V., Lahav, O., & Winter, M. K. 2016, ApJS, 225, 31

Lyutyi, V. M., & Cherepashchuk, A. M. 1972, Astron. Tsirkulyar, 688, 1

Mahabal, A., Rebbapragada, U., Walters, R., et al. 2019, PASP, 131, 038002

Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13

Marconi, A., Capetti, A., Axon, D. J., et al. 2001, ApJ, 549, 915

Marconi, A., Axon, D. J., Capetti, A., et al. 2003, ApJ, 586, 868

Marsden, C., Shankar, F., Ginolfi, M., & Zubovas, K. 2020, Front. Phys., 8, 61

McLure, R. J., & Jarvis, M. J. 2002, MNRAS, 337, 109

Melchior, A. L., & Combes, F. 2011, A&A, 536, A52

Merloni, A., Heinz, S., & di Matteo, T. 2003, MNRAS, 345, 1057

Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533

Muñoz-Mateos, J. C., Sheth, K., Regan, M., et al. 2015, ApJS, 219, 3

Narayan, R. 2005, Astrophys. Space Sci., 300, 177

Onishi, K., Iguchi, S., Davis, T. A., et al. 2017, MNRAS, 468, 4663

Pasquet, J., Bertin, E., Treyer, M., Arnouts, S., & Fouchez, D. 2019, A&A, 621, A26

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682

Salo, H., Laurikainen, E., Laine, J., et al. 2015, ApJSS, 219, 4

Scharwächter, J., Dopita, M. A., Shastri, P., et al. 2016, Astrophys. Space Sci. Proc., 42, 263

Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, MNRAS, 460, 3119

Shankar, F., Bernardi, M., Richardson, K., et al. 2019, MNARS, 485, 1278

She, R., Ho, L. C., & Feng, H. 2017, ApJ, 835, 223

Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45

Silk, J., & Rees, M. J. 1998, A&A, 331, L1

Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., & Salakhutdinov, R. 2014, J. Mach. Learn. Res., 15, 1929

Thean, A., Pedlar, A., Kukula, M. J., Baum, S. A., & O'Dea, C. P. 2000, MNRAS, 314, 573

Toomre, A. 1964, ApJ, 139, 1217

van den Bosch, R. C. E. 2016, ApJ, 831, 134

Vestergaard, M. 2002, ApJ, 571, 733

Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689

Wandel, A., Peterson, B. M., & Malkan, M. A. 1999, ApJ, 526, 579

Wang, J.-M., & Zhang, E.-P. 2007, ApJ, 660, 1072

Wang, R., Carilli, C. L., Beelen, A., et al. 2007, AJ, 134, 617

Wang, J.-G., Dong, X.-B., Wang, T.-G., et al. 2009, ApJ, 707, 1334

Wright, S. 1921, J. Agricul. Res., 20, 557

Xiao, T., Barth, A. J., Greene, J. E., et al. 2011, ApJ, 739, 28

<footnote n="1">¹ Observatoire de Paris, LERMA, PSL University, Sorbonne Université, CNRS, F-75014 Paris, France</footnote>

<footnote n="2">² Institute of Physics, Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne (EPFL), 1290 Sauverny, Switzerland</footnote>

<footnote n="3">³ Collège de France, 11 Place Marcelin Berthelot, 75231 Paris, France</footnote>

<footnote n="4">⁴ Observatorio Astronómico Nacional (OAN-IGN)-Observatorio de Madrid, Alfonso XII, 3, 28014 Madrid, Spain</footnote>

<footnote n="5">⁵ Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="6">⁶ Instituto de Astrofísica de Canarias (IAC), Calle Vía Láctea, s/n, 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="7">⁷ Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="8">⁸ Departamento de Física de la Tierra y Astrofísica, Fac. de CC Físicas, Universidad Complutense de Madrid, 28040 Madrid, Spain</footnote>

<footnote n="9">⁹ Instituto de Física de Partículas y del Cosmos IPARCOS, Fac. CC Físicas, Universidad Complutense de Madrid, 28040 Madrid, Spain</footnote>

<footnote n="10">¹⁰ Astronomical Institute, Academy of Sciences, Bôcní II 1401, 14131 Prague, Czech Republic</footnote>

<footnote n="11">¹¹ Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="12">¹² Max-Planck-Institut für Extraterrestrische Physik, Postfach 1312, 85741 Garching, Germany</footnote>

<footnote n="13">¹³ Institute of Astrophysics, Foundation for Research and Technology-Hellas, 71110 Heraklion, Greece</footnote>

<footnote n="14">¹⁴ School of Sciences, European University Cyprus, Diogenes Street, Engomi 1516, Nicosia, Cyprus</footnote>

<footnote n="15">¹⁵ Instituto de Radioastronomía y Astrofísica (IRyA), Universidad Nacional Autónoma de México, Antigua Carretera a Pátzcuaro #8701, Ex-Hda. San José de la Huerta, C.P. 58089 Morelia, Michoacán, Mexico</footnote>

<footnote n="16">¹⁶ Department of Physics & Astronomy, University of Alaska Anchorage, Anchorage, AK 99508-4664, USA</footnote>

<footnote n="17">¹⁷ School of Physics & Astronomy, University of Southampton, Southampton SO17 1BJ, Hampshire, UK</footnote>

<footnote n="18">¹⁸ INAF - Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="19">¹⁹ National Astronomical Observatory of Japan, National Institutes of Natural Sciences (NINS), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="20">²⁰ Department of Astronomy, School of Science, Graduate University for Advanced Studies (SOKENDAI), Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="21">²¹ Instituto de Física Fundamental, CSIC, Calle Serrano 123, 28006, Madrid, Spain</footnote>

<footnote n="22">²² Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile</footnote>

<footnote n="23">²³ Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China</footnote>

<footnote n="24">²⁴ School of Mathematics, Statistics, and Physics, Newcastle University, Newcastle upon Tyne NE1 7RU, UK</footnote>

<footnote n="25">²⁵ LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Paris-Cité University, 5 Place Jules Janssen 92190, Meudon, France</footnote>

<footnote n="26">²⁶ Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>
Appendix A: Moment maps and position-velocity diagrams

<figure><img src="image_7.png" /><figcaption>Fig. A.1. Left we present the first three moments maps of NGC 4388, obtained from the best-fit model, giving the adopted $M_{BH}$, compared with the NGC 4388 CO(3-2) molecular gas deduced from the García-Burillo et al. (2021) ALMA observations. From left to right, the integrated intensity, the velocity field, and the velocity dispersion; the top is the observation, the bottom is the best-fit model. On the right, we present the observed PVD, with the best-fit model overlaid in red contours. The parameters used for the fit are presented in Table A.1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. A.2. Same as A.1 for NGC 5506.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. A.3. Same as A.1 for NGC 5643.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. A.4. Same as A.1 for NGC 6300.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. A.5. Same as A.1 for NGC 7314.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. A.6. Same as A.1 for NGC 7465.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. A.7. Same as A.1 for NGC 7582.</figcaption></figure>
<table><caption>**Table A.1.** Parameters of the moment maps fits</caption><thead><tr><th>PARAMETER</th><th>NGC4388</th><th>NGC5506</th><th>NGC5643</th><th>NGC6300</th><th>NGC7314</th><th>NGC7465</th><th>NGC7582</th></tr></thead><tbody><tr><td><i>M</i><sub>bulb</sub></td><td>0.40</td><td>0.20</td><td>0.10</td><td>0.25</td><td>0.10</td><td>0.60</td><td>0.45</td></tr><tr><td><i>R</i><sub>bulb</sub></td><td>0.10</td><td>0.20</td><td>0.10</td><td>0.10</td><td>0.10</td><td>0.15</td><td>0.10</td></tr><tr><td><i>H</i><sub>bulb</sub></td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td><i>M</i><sub>gal</sub></td><td>8.0</td><td>4.0</td><td>8.0</td><td>9.0</td><td>3.0</td><td>1.0</td><td>5.0</td></tr><tr><td><i>R</i><sub>0</sub></td><td>3.0</td><td>3.0</td><td>1.0</td><td>3.0</td><td>3.0</td><td>10.0</td><td>8.0</td></tr><tr><td><i>H</i><sub>t</sub></td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.003</td><td>0.001</td><td>0.002</td><td>0.001</td></tr><tr><td>log <i>M</i><sub>BH</sub></td><td>6.61</td><td>6.63</td><td>6.99</td><td>6.39</td><td>6.41</td><td>6.62</td><td>6.40</td></tr><tr><td><i>RAPM</i></td><td>0.020</td><td>0.021</td><td>0.020</td><td>0.021</td><td>0.020</td><td>0.020</td><td>0.020</td></tr><tr><td><i>R</i><sub>gas</sub></td><td>0.48</td><td>0.50</td><td>0.55</td><td>0.20</td><td>0.20</td><td>0.20</td><td>0.60</td></tr><tr><td><i>H</i><sub>gas</sub></td><td>0.12</td><td>0.012</td><td>0.011</td><td>0.03</td><td>0.02</td><td>0.12</td><td>0.08</td></tr><tr><td><i>PA</i></td><td>10</td><td>0</td><td>[50,150]</td><td>160</td><td>80</td><td>30</td><td>110</td></tr><tr><td><i>IA</i></td><td>79</td><td>80</td><td>29</td><td>57</td><td>55</td><td>55</td><td>59</td></tr></tbody></table>

**Notes.** Table of the parameters used to obtain the fits of the observed moment maps shown in the figures A.1 to A.7. The 12 parameters are described in the fourth paragraph of the 4.2. For each galaxy, we choose the $M_{BH}$ given by the estimations made with moment maps model. In the case of NGC 5643, we have two different PA in order to model the change of PA seen in the observation. The RAPM parameter has no units, $M_{bulb}$ and $M_{gal}$ are in units of $2.25 \times 10^9 M_{\odot}$, $M_{BH}$ in $log(M/M_{\odot})$, PA and IA are in degrees and $R_{bulb}$, $H_{bulb}$, $R_0$, $H_t$, $R_{gas}$, and $H_{gas}$ are in kpc.