© The Authors 2025

1. Introduction

The centers of massive galaxies host supermassive black holes (SMBHs) and black holes (BHs) with masses ranging from 10⁶ M_⊙ to 10¹⁰ M_⊙. 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–σ 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 Λ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 matching, Behroozi et al. (2013) observe that around z ∼ 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% 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–σ 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–σ 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–σ 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)≥10⁴² 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. Since the symmetrization tends to keep the highest gradients on the PVD, we can expect that the M_BH estimations might be biased high.

Fig. 1. Example of NGC 4388 symmetrized PVD.

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 -σ relation

The tight correlation between BH and bulge masses drawn from the M_BH–σ 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–σ relation (e.g., Silk & Rees 1998; Granato et al. 2004). Nevertheless, the observational M_BH–σ relation is not well constrained for low-mass galaxies (M < 10⁹ M_⊙). 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–σ relation:

$$\begin{array}{c}\hfill \log M_{\mathrm{BH}}(M_{\odot })=8.5\pm 0.05+(4.41\pm 0.29)\log \left(\frac{\sigma }{200\mathrm{km}/\mathrm{s}}\right),\end{array}$$(1)

where σ 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:

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=\frac{f(\mathrm{\Delta }{V}^2)R}{G},\end{array}$$(2)

where M_BH is the mass of the central BH, G is the gravitational constant, Δ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 ∼ 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

$$\begin{array}{c}\hfill log\left(\frac{M_{\mathrm{BH}}}{M_{\odot }}\right)=a+blog\left(\frac{L}{{10}^{44}{\mathrm{erg}\mathrm{s}}^{-1}}\right)+clog\left(\frac{\mathrm{FWHM}}{{\mathrm{km}\mathrm{s}}^{-1}}\right),\end{array}$$(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β is the most commonly used line. For redshift between 0.75 and 2, the Mg IIλ2800 line is a common choice (McLure & Jarvis 2002). And for z ≳ 2, the best line to use is C IVλ1549 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 Ṁ_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

$$\begin{array}{cc}\hfill \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}/\mathrm{s}}\right)\hfill \\ \hfill & +(-0.{59}_{-0.15}^{+0.16})\log \left(\frac{L_x}{{10}^{40}\mathrm{erg}/\mathrm{s}}\right),\hfill \end{array}$$(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 = νL_ν = (ν = 5 GHz)×4πD²F₅, F₅ 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–σ 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⁶.

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₀, 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₀ 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,

$$\begin{array}{c}\hfill M_0={\displaystyle \int I_{\nu }}\mathrm{d}\nu .\end{array}$$(5)

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

$$\begin{array}{c}\hfill M_1=\frac{{\displaystyle \int V}{I}_{\nu }\mathrm{d}\nu }{M_0}.\end{array}$$(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.

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 overfitting (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:

$$\begin{array}{c}\hfill \mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n}{(Y_{\mathrm{pred},i}-Y_{\mathrm{real},i})}^2}{n}},\end{array}$$(7)

with n the number of predictions made, Y_pred, i the predicted values, and Y_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)¹. It is a statistical value that represents the proportion of a dependent variable variance that is predictable from the independent variables:

$$\begin{array}{c}\hfill R^2=1-\frac{{\displaystyle {\sum }_{i=1}^n}{(Y_{\mathrm{pred},i}-Y_{\mathrm{real},i})}^2}{{\displaystyle {\sum }_{i=1}^n}{(\overline{Y_{\mathrm{pred}}}-Y_{\mathrm{real},i})}^2},\end{array}$$(8)

where the notations are the same as the RMSE ones and $\overline{Y_{\mathrm{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² 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):

$$\begin{array}{c}\hfill R{E}_i=\frac{Y_{\mathrm{pred},i}-Y_{\mathrm{real},i}}{Y_{\mathrm{real},i}}.\end{array}$$(9)

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

Fig. 2. Histograms of the RE of our models when estimating the MBH of the 5000 simulations from their respective test set. We represent in blue the models trained to estimate the MBH with the PVD, and in red the models trained with the moment maps.

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_⊙) 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_⊙) = 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_⊙)∼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_⊙) 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⁵ 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.

Fig. 3. Probability distributions of the MBH 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.

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.

Fig. 4. Comparison of the probability distributions of the MBH 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.

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

Fig. 5. 2D and 1D histograms of the RE from linear interpolation and chi squared fitting of MBH on the test set for NGC 7582. On the left, we plotted the RE when estimating the MBH 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(\frac{M_{\mathrm{BH}}}{M_{\odot }})$ and the 1D histograms axis is in number of estimations.

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 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∼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 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_⊙). 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_⊙) 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 iso-velocity 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 simulations 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 methodological consistency in estimating M_BH 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_BH per se but is instead learning to estimate M_BH 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 simulating 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_BH 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 developed 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_BH based on the input. To test the reliability of our models, we used three statistical indicators: the RMSE, the R², and the $\overline{\mathit{RE}}$. These statistical indicators confirmed that our fourteen models can make reliable M_BH estimations on new data. we looked at our models RE distribution when predicting the M_BH of the test sample. we saw that the distributions were symmetrical and centered around 0, indicating no constant over- or underestimation. We also studied the models’ RE of their M_BH 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_BH. This is a strong indication that the models M_BH predictions only depend on the M_BH. We studied the behavior of the M_BH predictions in depth by looking at the predictions RE versus the M_BH.

We managed to make a M_BH estimation with both architectures for seven out of the ten galaxies in the GATOS core sample: NGC 4388, NGC 5506, NGC 5643, NGC 6300, NGC 7314, NGC 7465, and NGC 7582. We made 10⁵ predictions per galaxy using MC dropout. The resulting M_BH estimations are Gaussian-like probability distributions. The scatter of these distributions 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 conditions. 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 represents the first step toward an automatized method for estimating M_BH.

In this study, we demonstrate that even with training sets simulated using simple physics, fixed IA and PA, and basic NN architectures, 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 developing ML techniques for M_BH estimations that can improve scaling with the increasing number of astronomical observations.

Data availability

Appendices B and C, which present the RE of M_BH predictions 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.

Acknowledgments

We would like to thank the referee for very useful comments, that helped us improve and clarify the paper. RP acknowledges financial 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 Hubble 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 funding 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 CONACyT “Frontera de la Ciencia” project CF-2023-G-100. BGL acknowledges support from grants PID2019-107010GB-100 and the Severo Ochoa CEX2019-000920-S funded by MICINN-AEI/10.13039/501100011033. CRA acknowledges 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 NextGenerationEU/PRTR.

References

Appendix A: Moment maps and position-velocity diagrams

Fig. A.1. Left we present the first three moments maps of NGC 4388, obtained from the best-fit model, giving the adopted MBH, 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.

Fig. A.2. Same as A.1 for NGC 5506.

Fig. A.3. Same as A.1 for NGC 5643.

Fig. A.4. Same as A.1 for NGC 6300.

Fig. A.5. Same as A.1 for NGC 7314.

Fig. A.6. Same as A.1 for NGC 7465.

Fig. A.7. Same as A.1 for NGC 7582.

All Tables

All Figures

	Fig. 1. Example of NGC 4388 symmetrized PVD.
In the text

	Fig. 2. Histograms of the RE of our models when estimating the MBH of the 5000 simulations from their respective test set. We represent in blue the models trained to estimate the MBH with the PVD, and in red the models trained with the moment maps.
In the text

	Fig. 3. Probability distributions of the MBH 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.
In the text

	Fig. 4. Comparison of the probability distributions of the MBH 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.
In the text

Fig. 5. 2D and 1D histograms of the RE from linear interpolation and chi squared fitting of MBH on the test set for NGC 7582. On the left, we plotted the RE when estimating the MBH 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(\frac{M_{\mathrm{BH}}}{M_{\odot }})$ and the 1D histograms axis is in number of estimations.

In the text

Fig. A.1. Left we present the first three moments maps of NGC 4388, obtained from the best-fit model, giving the adopted MBH, 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.

In the text

	Fig. A.2. Same as A.1 for NGC 5506.
In the text

	Fig. A.3. Same as A.1 for NGC 5643.
In the text

	Fig. A.4. Same as A.1 for NGC 6300.
In the text

	Fig. A.5. Same as A.1 for NGC 7314.
In the text

	Fig. A.6. Same as A.1 for NGC 7465.
In the text

	Fig. A.7. Same as A.1 for NGC 7582.
In the text