Comparing indirect methods for black hole masses in AGN: the good, the bad, and the ugly

The black hole mass MBH is crucial in constraining the growth of supermassive BHs within their host galaxies. Since direct measurements of MBH with dynamical methods are restricted to a limited number of nearly quiescent nearby galaxies and a small minority of active galactic nuclei (AGN), we must rely on indirect methods. In this work, we utilize an unbiased, volume-limited, hard X-ray selected sample of AGN to compare the reliability of some commonly used indirect methods, emphasising those that can be applied to obscured AGN. Based on a subsample of AGN with MBH determined via dynamical methods, our study suggests that X-ray based techniques, such as the scaling method and the one based on the variability measured through the excess variance, are in good agreement with the dynamical methods. On the other hand, the M-sigma correlation based on inactive galaxies tends to systematically overestimate MBH, regardless of the level of obscuration. We provide a correcting factor that produces an acceptable agreement with dynamical values and can be used to quickly correct the MBH computed with this method. We also derive an alternative M-sigma correlation based on this unbiased sample of AGN with a slope considerably shallower than the ones obtained using inactive galaxies, suggesting that the latter correlation may not be appropriate to compute the MBH in AGN. Finally, we find that no quick fix can be applied to correct the MBH obtained from the fundamental plane of black hole activity, casting doubts on the reliability of this method.


INTRODUCTION
It is now widely accepted that Active Galactic Nuclei (AGN), the most luminous persistent sources in the universe, are powered by accretion onto a supermassive black hole at the center of their host galaxy.Depending on the black hole mass,  BH , and the accretion rate, , their bolometric luminosity  bol may range from a tiny fraction of their host galaxy luminosity, as in the case of our Galaxy, to orders of magnitude beyond the host galaxy, as in the case of the most powerful quasars.Since the accretion rate is generally measured in terms of Eddington ratio,  Edd =  bol / Edd (where the Eddington luminosity is defined as  Edd = 1.3 × 10 38 ( BH /M ⊙ ) erg s −1 ), the  BH plays a crucial role in characterizing the intrinsic properties of AGN.
As for all astrophysical sources, the most reliable way to determine the mass of black holes is through dynamical methods: by measuring the orbital parameters of objects or gas gravitationally bound to the BH, one can constrain  BH via the virial theorem (e.g., Ferrarese & Ford 2005).For a few low-power AGN, close enough to have their gravitational sphere of influence spatially resolved by current instrumentation, the  BH may be determined from stellar or gas dynamics measurements (e.g., Gebhardt et al. 2003;Macchetto et al. ★ E-mail: mgliozzi@gmu.edu1997).For brighter unobscured AGN (the so-called type 1 AGN),  BH can be reliably constrained using the reverberation mapping (RM) technique, where the delay between the optical continuum, associated with the accretion disk, and the Balmer lines, produced in the broad line region (BLR), yields the size of the BLR, and the width of the Balmer lines constrains the BLR velocity (e.g., Blandford & McKee 1982;Peterson et al. 2004).Unfortunately, the RM method is heavily resource-and time-consuming and can only be employed for broad-lined AGN that vary on a reasonably short amount of time, which severely limits its applicability.Luckily, the discovery of a tight correlation between the optical luminosity and the size of the BLR (e.g., Kaspi et al. 2005;Bentz et al. 2006) has made it possible to extend the  BH determination to a much larger number of broadlined AGN, using the so-called single-epoch (SE) technique, where the measurement of the line-continuum delay is no longer necessary (e.g., Laor 1998;Wandel et al. 1999).Although this indirect method has proven very successful in extending the estimate of the  BH to thousands of broad-lined AGN, some caution should be considered when this technique is applied well beyond the range probed by the RM sample, such as at high redshift (e.g., Denney 2012;Denney et al. 2016) or in sources where there is substantial obscuration (see, e.g., Mejía-Restrepo et al. 2022 and references therein).
Indeed, there is a major problem in accurately constraining  BH in obscured AGN (type 2 AGN).In this case, the only direct dynamical method available is based on the rare detection of water masers in Keplerian motion within the accretion disk surrounding the BH.This technique has been successfully applied only to a fairly limited number of obscured AGN (e.g., Greenhill et al. 2003;Kuo et al. 2011).In addition, for a few partially obscured AGN with broad lines observable in the near-infrared band, a variation of the SE method, based on the correlation between BLR size and hard X-ray luminosity, was successfully utilized by Ricci et al. (2017a,b) and Onori et al. (2017).However, for the vast majority of obscured AGN, which represent the bulk of the AGN population, no simple direct method to constrain  BH exists, and generally one has to rely on statistical correlations between  BH and some properties associated with the bulge of the host galaxy, such as mass, luminosity, and velocity dispersion  ★ (e.g., Kormendy & Richstone 1995;Magorrian et al. 1998;Gebhardt et al. 2000;Ferrarese & Merritt 2000).Since these correlations have been established using samples of mostly quiescent galaxies and unobscured AGN in the local universe, it is important to test their validity to a larger population of AGN with particular focus on the obscured ones.
In our recent work, using NGC 4151, one of the few AGN with  BH determined with multiple direct dynamical methods, we carried out a comparison between indirect methods that can be used to constrain  BH in obscured AGN (Williams et al. 2023).Here, we expand our previous investigation by performing a systematic comparison between indirect methods using an unbiased sample that can be considered as representative of the AGN population in the local universe.More specifically, we will compare the - ★ correlation method (Kormendy & Ho 2013), the fundamental plane of BH activity (Gültekin et al. 2019), the X-ray scaling method (Shaposhnikov & Titarchuk 2009;Gliozzi et al. 2011), and the method based on X-ray variability as measured by the normalized excess variance (e.g., Papadakis 2004;Ponti et al. 2012;Akylas et al. 2022).For completeness, we will also verify whether the optically based single epoch method is fully consistent with dynamical measurements for broad-lined AGN contained in our sample.To this aim, we utilize a volume-limited sample of AGN, whose selection criteria -very hard X-ray detection, short distance  < 40 Mpc, and X-ray bright with  X > 10 42 erg s −1 -ensure that only bona fide AGN that are unbiased in terms of obscuration are included, and that a sizable fraction of AGN has  BH estimated with direct dynamical methods.Importantly, the near totality of the sample has been observed with NuSTAR which, thanks to its high sensitivity in the hard energy band (3-79 keV) and low background, is ideal to robustly constrain the spectral contributions of absorption and reflection and hence the primary emission, and to measure the intrinsic variability of the AGN.
The structure of the paper is the following.In Sections 2 and 3, we describe respectively the sample selection and the data reduction and analysis of NuSTAR data, as well as that of radio data.In Section 4, we report the  BH values obtained with the different indirect methods, and in Section 5 we carry out a statistical analysis and systematic comparison of these methods.In Section 6, we discuss the main findings and draw our conclusions in Section 7.
In brief, we found that the X-ray scaling method and the X-ray variability one yield  BH values that are fully consistent with the dynamical ones.On the other hand the  BH values obtained using the - ★ correlation of Kormendy & Ho (2013) appear to be systematically overestimated.We provide a correction factor and derive an alternative - ★ correlation with a shallower slope that yields values consistent with the dynamical estimates of this volume-limited sample of AGN.We also confirm that the  BH estimates derived with the FP are unreliable and there is not a simple correction that can make them consistent with the dynamical ones.Finally, for the broad-lined AGN of our sample, we verified that the optically based SE method provides values broadly consistent with the dynamical ones with a tendency to underestimate them in some (but not all) of the AGN that are substantially obscured.

SAMPLE SELECTION
For our work we utilized AGN with hard X-ray luminosity  14−195 keV > 10 42 erg s −1 , taken from Oh et al. (2018), and redshift-independent distance  < 40 Mpc from Koss et al. (2022a).Our sample was selected from the Swift 70-month Burst Alert Telescope (BAT) catalog (Baumgartner et al. 2013), where the detection of AGN at very hard X-ray energies (14-195 keV) ensures that the sample is almost independent of obscuration.Our choice of a hard X-ray luminosity threshold guarantees that all the selected objects are genuine AGN, and the distance criterion makes sure that the sample is not biased towards the brightest AGN.When dealing with nearby AGN, one must consider the potential error associated with the distance, since their actual values can be substantially different from the distances obtained from their radial velocity assuming a uniform expansion and Hubble's law.Luckily, the second data release of the Swift BAT AGN Spectroscopic Survey (BASS DR2, Koss et al. 2022a) provides high-quality redshift-independent distances for all nearby AGN, determined following the approach of Leroy et al. (2019), which we adopted as bona fide distances.
Following these criteria our sample is made up of 15 broad-lined AGN with five of them classified as Seyfert 1.9 galaxies (hereafter, Sy1.9), whose optical spectra are characterized by a narrow H  and a broad H  lines, and 20 narrow-lined AGN (hereafter, Sy2).After excluding the sources without archival NuSTAR data (NGC 4235 and LEDA 89913) and one source that is likely to be jet-dominated (Cen A), our resulting sample is made up of 32 objects.
Since our goal is to test the reliability of  BH estimates obtained with various indirect methods, we searched the literature for direct dynamical measurements and found 11 objects whose  BH was determined either via maser measurements (five objects), gas dynamics (two objects), or reverberation mapping with direct modeling of the BLR via velocity-delay maps (four objects).These 11 objects define our restricted dynamical sample, which comprises the most reliable  BH estimates to be compared with the values obtained with different indirect methods.Ideally, we would compare all the indirect methods tested in this work using this restricted sample.However, some of these methods (namely, the X-ray variability method and the single epoch method based on broad H  line measurements) have  BH estimates for less than ten objects in the restricted dynamical sample.For this reason, we also defined an extended dynamical sample that comprises six additional AGN (for a total of 17 objects with dynamical mass), of which four have the  BH determined with the single epoch method in the near-IR (Ricci et al. 2017a) and two objects with the traditional reverberation mapping, where a virial factor ⟨  ⟩ = 4.3 was adopted.In the remainder of the paper, we will explicitly state when a result refers to the restricted dynamical sample or to the extended dynamical sample.
The basic characteristics (name, BAT id, classification, distance, dynamically determined  BH , and relative reference) of our AGN sample are reported in Table 1.A single asterisk after the name denotes our restricted dynamical sample of 11 objects.A double asterisk denotes the sources added to the restricted sample to form our extended sample of 17 objects.The source for the distances is Koss et al. (2022a).All  BH values obtained with methods that depend on the distance were renormalized using the distances quoted in this table.

DATA REDUCTION
We processed the entire sample with the NuSTAR data analysis script nupipeline, which calibrates and screens the data.Our version of the calibration database CALDB is 20200429.We defined each source region as a circle centered on the source and each background region as a circle nearby in the same image in an area that did not appear to contain source photons.Depending on the brightness of the source, our circles had radii from 30 to 100 arcsec.We used both focal plane modules FPMA and FPMB on the satellite to produce spectra and light curves with the nuproducts script.We list all the NuSTAR observations used in this paper in Table A1 in the Appendix, where we report the effective exposures (which by definition are shorter than the observation durations).In the spectral analysis for sources with multiple NuSTAR exposures, we chose the observations with longer exposure, higher mean count rate, and relatively low variability.The first two criteria favor higher signal-to-noise spectra, whereas the latter reduces the possibility of spectral variability affecting the timeaverage spectra (see Williams et al. 2023 for more details).For the temporal analysis, following Akylas et al. (2022) all the observations with length longer than 80 ks were utilized.

X-ray spectral analysis
One of the  BH estimation methods we wish to test is the X-ray scaling method, which is based on the spectral fitting of the primary X-ray emission of the objects in our sample with the bulk Comptonization model BMC (Titarchuk et al. 1997).To accurately constrain the primary X-ray emission produced in the corona in heavily obscured sources (such as the majority of our sample), it is important to properly parameterize the contributions of absorption and reflection associated with the torus.This goal is achieved by using the broadband X-ray spectrum afforded by NuSTAR and employing physically motivated models of the torus.As in our previous works focused on sources that are part of this sample (Gliozzi et al. 2021;Shuvo et al. 2022), we used as baseline model phabs * (atable(Borus) + MYTZ*BMC) where the first absorption model phabs parameterizes the Galaxy contribution, the borus table model describes the continuum scattering and fluorescent emission line components associated with the torus (Baloković et al. 2018), and the zeroth-order component of MYTorus, MYTZ, models the absorption and Compton scattering that act on the transmitted primary emission (Murphy & Yaqoob 2009).
As explained in more detail in Gliozzi et al. (2021), we use MYTZ be-cause it models the Compton scattering with the Klein-Nishina cross section, which plays an important role at energies above 10 keV in heavily absorbed AGN.The Comptonization model BMC has four parameters: the normalization  BMC , the spectral index  (linked to the photon index by Γ =  + 1), the temperature of the seed photons , and log , which is related to the fraction of scattered seed photons  by the relationship  =  /(1 −  ) (note that log  has very little impact on the determination of  BH ).
To maintain the self-consistency of the torus model, the spectral index of the bmc model was linked to the photon index value of the torus model and  BMC was forced to be equal to the borus powerlaw normalization divided by a factor of 30, based on the empirical relationship obtained in Gliozzi et al. (2021).
For the few objects with relatively low absorption, the MYTZ component is substituted by a zphabs model left free to vary.For a few sources with more complex spectra, additional components were added to the baseline model; in particular, we included a model describing the fraction of primary emission that is directly scattered towards the observer by an optically thin ionized medium (const*bmc).
The main spectral parameters obtained by fitting this baseline model and the 2-10 keV luminosity are reported in Table A2 in the Appendix, along with the details on the sources that required more complex models, which were not already published in our previous works.We performed the X-ray spectral analysis using the xspec v.12.9.0 software package (Arnaud 1996), and the errors quoted on the spectral parameters represent the 1 confidence level.

X-ray temporal analysis
Another indirect method that we want to test is based on X-ray variability.For this reason, we extracted NuSTAR light curves in the 10-20 keV band and estimated the normalised excess variance based on the prescription presented in Akylas et al. (2022).In brief, the light curves obtained for the FMPA and FMPB modules were first background subtracted and then summed using the LCMATH tool within FTOOLS.To measure the variability power of the sources in our sample, we compute the normalised excess variance (e.g., Nandra et al. 1997) using: where  is the number of bins in the light curve, X i and  i are the count rates and their uncertainties, respectively, and  is the unweighted mean count rate.
The 2   has been measured over light curve segments with a duration of Δt = 80 ks.When more than one valid segment is available for a source, we compute the mean of the individual excess variances.

Radio analysis
For the majority of our sources, we either utilized the 5 GHz VLA peak intensities (F peak ) provided by Fischer et al. (2021) or used the C-band (∼ 4.89 GHz) VLA archival1 images to measure the peak intensities, whenever available in A-Configuration 2 .Similarly to Williams et al. (2023), we used VLA (kpc scale) measurements in our analysis because AGN do not follow the fundamental plane of black hole activity when radio intensities are measured on milliarcsecond (subparsec) scales from VLBA observations (Fischer et al. 2021;Shuvo et al. 2022).The declination coverage of VLA is restricted to  > −48 • ; therefore, for some of the sources below this threshold (NGC 1566, ESO 5-4, NGC 6221, NGC 6300) we have used images at 4.8 GHz from the 64-meter Parkes radio telescope, and for a few others (NGC 4945, NGC 7213, Circinus) we have instead used images at 5 GHz from the Australia Telescope Compact Array (ATCA), which are available on the NASA Extragalactic Database (NED)3 .We could not find any archival images at 5 GHz for NGC 678, NGC 5290, or ESO 137−34.
We used the NRAO Common Astronomical Software Applications, also known as casa (McMullin et al. 2007) 4 , to analyze the images obtained from the VLA archive or NED.In casa, we utilized the two-dimensional fitting application through the "viewer" command to determine the peak intensities directly from the image.The radio luminosities (erg s −1 ) shown in column 7 of Table 2 for all sources were calculated using the redshift-independent distances provided in Table 1.

BLACK HOLE MASS MEASUREMENTS
In the vast majority of AGN, dynamical methods cannot be applied either because they are too distant and hence their BH sphere of influence is too small to be spatially resolved or because the central engine and more specifically the BLR is not directly observable.In those cases, indirect methods have a crucial importance.Here, we briefly describe a few of these indirect methods whose accuracy and reliability will be tested against dynamical mass estimates in the following section.We will mainly focus on methods that can be applied to obscured AGN, that is, objects for which the BLR is not detected or hardly detected.However, for comparison purposes, we also test the optically based single epoch method.The black hole mass values obtained with the different indirect methods are reported in Table 2.

Single epoch spectra
The single epoch (SE) method is the most utilized method to estimate  BH of broad-lined AGN at all redshifts.The  BH is obtained from the same formula used in the reverberation mapping method: where Δ is the velocity dispersion of a specific broad emission line (generally measured from full width at half maximum or from the standard deviation  of the Gaussian model used to fit the line),  is the size of the BLR,  the universal gravitational constant, and  the virial factor, which encompasses the unknown geometry and inclination of the BLR.
The key difference between SE and RM is that in the latter case  is obtained from the time delay between variations of the optical/UV ionizing continuum and those observed in the BLR ( = Δ), whereas in the SE method the BLR size is directly obtained from the tight correlation between  and a monochromatic continuum luminosity in the optical band  opt , which was empirically derived from reverberation mapping studies (e.g., Kaspi et al. 2005;Bentz et al. 2006).
This method has been successfully applied to thousands of broadlined AGN.However, some caution should be used since the average virial factor ⟨  ⟩ introduces a systematic uncertainty of ∼ 0.4 dex (e.g., Pancoast et al. 2014).Additionally, although the correlation - opt was established for a nearby sample of RM AGN, the same correlation is extended to very distant AGN using high-ionization lines like C iv, which are often blue-shifted because of winds and hence not reliable tracers of the BLR dynamics (e.g., Trakhtenbrot & Netzer 2012).Finally, it has been recently suggested that the optically based SE method tends to underestimate the  BH in AGN with substantial obscuration (e.g., Caglar et al. 2020;Mejía-Restrepo et al. 2022).For these reasons, we have included the SE method in our statistical comparison with the dynamical methods.
To be more specific, in this work we compute the  BH via the SE method using the broad H α values of FWHM(H α) and L(H α) provided by Mejía-Restrepo et al. (2022) in their table 5, as well as their equation 2:

M-𝜎 ★ relation
Over the last few decades, numerous studies focusing on nearby galaxies have revealed the existence of several relatively tight correlations between properties of the galaxy bulge (for example, luminosity, mass, and velocity dispersion  ★ ) and the mass of the supermassive BH at the center of the galaxy, suggesting a co-evolution between these two elements operating at very different scales (see, e.g., Kormendy & Ho 2013, for a detailed review).Since - ★ is the tightest among these correlations, it has been frequently used to estimate  BH in the absence of more direct dynamical methods.
In our recent work, where we compared different indirect methods using NGC 4151 as a laboratory, we tested different versions of this correlation (Williams et al. 2023).More specifically, we compared the dynamical value of  BH with those derived with the - ★ correlation of McConnell & Ma (2013), which makes a distinction between early and late type galaxies, with the equation of Woo et al. (2013), which was calibrated using an RM sample of AGN, with the correlation of Kormendy & Ho (2013), which was derived from a sample of galaxies with classical bulges and with  BH determined with direct dynamical methods, and finally with the correlation proposed by Shankar et al. (2016), who, assuming that previous correlations were biased towards higher values of  ★ , proposed a new equation that they named the intrinsic correlation.Of all the - ★ correlations tested, only the one proposed by Kormendy & Ho (2013) yielded a value in agreement with the dynamical estimate of NGC 4151.For this reason, and because the same correlation is used by Koss et al. (2022b) to estimate the  BH for the vast majority of the AGN in BASS DR2, throughout this work we will use the Kormendy & Ho (2013)  (4)

Fundamental plane of BH activity
The fundamental plane of black hole activity (hereafter, FP) is a correlation involving X-ray luminosity, radio luminosity, and  BH , which was independently proposed by Merloni et al. (2003) and Falcke et al. (2004) to unify the accretion and ejection properties of black hole systems at all scales.When other methods to constrain  BH are unavailable, the FP has been frequently used in AGN studies, despite its inherent large scatter (nearly 1 dex), which yields order of magnitude estimates rather than tight constraints on  BH .
As for the - ★ correlation, different versions of the FP were developed over the years and three of them were tested in Williams et al. (2023), where the FP proposed by Dong et al. (2014) for relatively highly accreting BH systems yielded the closest estimate to the dynamical mass of NGC 4151.However, when we tried to apply this FP to our volume-limited sample of AGN, we obtained unreasonable values with  BH ranging from a few solar masses to 10 11 M ⊙ .We therefore used the FP recently proposed by Gültekin et al. (2019), which was derived using only dynamically constrained  BH and is described by the equation below.
Although high-resolution VLBA data exist for most of the AGN in our sample, we utilized lower-resolution VLA data for our estimate of  BH based on the FP, because it was recently shown that using higher-resolution radio data the FP breaks down (Fischer et al. 2021;Shuvo et al. 2022).

X-ray scaling method
The X-ray scaling method, introduced by Shaposhnikov & Titarchuk (2009) for X-ray binary systems (XRBs) with an accreting BH, was first successfully applied to unobscured AGN with  BH determined via RM (Gliozzi et al. 2011) and then extended to heavily obscured AGN obtaining  BH consistent with those determined via maser measurements (Gliozzi et al. 2021), making it one of the few truly universal methods, since it can be equally utilized for stellar mass BHs and supermassive ones, as well as for obscured and unobscured BH systems.
In brief, this method determines  BH in AGN by scaling up the dynamically measured mass of a stellar-mass BH (the reference source) under the assumption that the primary X-ray emission from the corona is produced by the same Comptonization process in AGN and XRBs and that Γ can be considered a faithful proxy of the accretion rate of all BH systems.A detailed description of the method is provided in Gliozzi et al. (2021) and references therein; here we just report the equation used to determine the  BH .A derivation of this equation is provided in Williams et al. (2023), where the tabulated values of  BMC for the different reference sources are also provided.
where  BMC is the normalization of the Comptonization model bmc model (measured in units of 10 39 erg s −1 divided by the distance squared in units of 10 kpc) and  the distance in kpc.The application of this method to NGC 4151 in Williams et al. ( 2023) confirms the conclusions derived in previous works that there are two reference sources, GRO J1655-40 during the 2005 decay and GX 339-4 during the 2003 decay, which yield the most reliable estimates of  BH .For this reason, throughout this work we computed the  BH with both reference sources (when available) and then calculated the average, which is the value reported in the table.Only one source (Circinus) has Γ too steep for these reference sources.In that case, we utilized XTE J1550-564 during the 1998 rise as reference source, and multiplied the mass estimate by a factor of 4, which is the correction needed for this reference source to be consistent with dynamically measured  BH (Williams et al. 2023).

X-ray variability method
High persistent variability is one of the defining characteristics of AGN at all wavelengths.However, on relatively short timescales, it is particularly prominent in the X-ray band, which has the additional advantages with respect to the optical/UV bands of being less affected by absorption and more easily separated from star/galaxy contributions.The normalised excess variance  2 NXV , estimator of the intrinsic 'band-variance' of the source, (e.g., Nandra et al. 1997) has been extensively used as a predictor of the  BH of AGN.Past studies have demonstrated the existence of a tight anti-correlation between the normalized excess variance  2 NXV and  BH in both soft (0.5-2 keV) and hard (2-10 keV) energy ranges (e.g., Papadakis 2004;Ponti et al. 2012;Lanzuisi et al. 2014).Very recently, Akylas et al. (2022) extended this technique to the harder X-ray band afforded by NuSTAR, making easier the application of this technique to obscured AGN.To this end, Akylas et al. (2022) utilized a combined sample of AGN with  BH determined via RM and the - ★ correlation method and explored in detail the minimum requirements for the light curves that will be used to compute  2 NXV , so that the resulting BH mass estimates will be as unbiased as possible and of known variance.
Since in this work we want to assess the reliability of several indirect methods (including the - ★ correlation), and since the anti-correlation between  BH and  2 NXV derived by Akylas et al. ( 2022) is based on black hole mass values derived from the - ★ correlation and from the reverberation mapping technique (RM), we will only use the anti-correlation obtained with the RM sample.However, since the RM sample itself is calibrated on the - ★ correlation, we also derive a new anti-correlation using the  BH values of our extended dynamical sample.The use of the extended sample is necessary because the three Compton-thick sources (NGC 1068, NGC 4945, Circinus) and the changing-look AGN NGC 2992 with an anomalously low value of  2 NXV are excluded from the correlation analysis, further reducing the number of data points.Our results are shown in Fig. 1, where the green continuous line represents our best fit, which is similar but slightly offset with respect to the best fit obtained by Akylas et al. (2022) using reverberation mapping data only (short-dashed black line).For completeness, in this figure we also show the best fit obtained using the restricted sample of six objects (long-dashed red line) and that from the extended dynamical sample including the CL AGN (blue long-dashed line).
In the remainder of the paper, when we discuss the  BH estimated with the X-ray variability technique, we will quote both the ones obtained using the RM sample in Akylas et al. (2022) (whose anticorrelation is described by the short-dashed black line in Fig. 1) and those derived from the anti-correlation obtained with the extended dynamical sample in this work (whose anti-correlation is represented by the continuous green line in the same figure).

Comparison with the dynamical sample
A simple way to test the reliability of indirect methods is to compare their  BH distributions with that derived using dynamical methods.This is illustrated in Fig. 2, where the histograms yielded by the - ★ , FP, and X-ray scaling methods are compared to the histogram of the restricted dynamical sample (top panels), and the distributions of  BH produced by the SE and the X-ray variability method (using both the Akylas et al. (2022) prescription and the correlation derived in this work) are compared to the histogram of the extended dynamical sample (bottom panels).
From this figure, one can conclude that there is a reasonably good overlap between the dynamically based  BH distributions and the distributions obtained with most of the indirect methods tested here, with the notable exception of the distribution derived from the - ★ correlation, which appears to be shifted to the right toward larger values (top left panel).
In order to quantify the agreement (or lack thereof) between the  BH distributions obtained with the various indirect methods and the one inferred from dynamical measurements, we carried out Kolgomorov-Smirnov (K-S) and Student's  tests, which compare the whole distributions and their means, respectively.Large values of the probabilities associated with these tests indicate that the distributions (and means) are statistically consistent with the dynamical values.The results of these tests are reported in the first four columns    Columns. 1 = object name. 2 = full width at half maximum for the broad H α line.3 = log of the luminosity of the broad H α line.4 =  BH computed with the single epoch method.5 = stellar velocity dispersion in km s −1 .6 =  BH computed with the - ★ method of Kormendy & Ho (2013).7 = radio luminosity at 5 GHz in erg s −1 .8 =  BH computed with the FP method of Gültekin et al. (2019).9 = normalized excess variance computed with 80 ks segments in the 10-20 keV band.10 =  BH computed with the X-ray variability method.11 =  BH computed with the X-ray scaling method.
of Table 3 and are in general agreement with the qualitative conclusions inferred from the visual inspection of Fig. 2.
Although comparing the distributions of  BH offers a simple way to reveal whether an indirect method is inconsistent with dynamical measurements (as in the case of the - ★ method), it may be misleading in suggesting that other methods are fully consistent with the dynamical estimates.For example, if the same indirect method happens to yield for some objects  BH measurements that are substantially larger and for others substantially smaller than the accepted dynamical ones, the overall distribution may appear to be consistent with the dynamical one (as in the case of the FP method).In other words, the formal consistency between  BH distributions should be considered as a necessary but not sufficient condition for an indirect method to be considered reliable.
A more direct way to illustrate whether an indirect method is accurate is to plot the  BH obtained from indirect methods versus the corresponding dynamical ones, as shown in Fig. 3, where the continuous line indicates the perfect match  BH,indirect / BH,dyn = 1 and the dashed lines represent ratios of 3 and 1/3, which are of the order of the typical uncertainties associated to  BH estimates.This figure qualitatively confirms that there is a good agreement between dynamically estimated  BH and values obtained with the X-ray scaling method (top right panel), the SE method (bottom left panel), as well as with the values obtained with the X-ray normalized excess variance (bottom middle and right panels).In the latter case, for completeness, we have shown both the values obtained with the best fit correlation obtained by Akylas et al. (2022) using a larger RM sample (bottom middle panel) and the best fit obtained in this work using the extended dynamical sample (bottom right panel).Both choices provide values in general agreement with the dynamical ones.A quantitative way to assess the agreement (or lack thereof) illustrated in Fig. 3 is based on the Spearman's and Kendall's rank correlation coefficients  and  and their associated probabilities, with a positive correlation expected for the indirect methods that yield  BH estimates consistent with the dynamical ones.The results of these statistical tests are reported in the fifth and sixth columns of Table 3 and corroborate the previous conclusions with negative correlation coefficients obtained for the - ★ and FP methods and positive correlation coefficients for the SE and the two X-ray based methods.
An additional straightforward way to directly compare the  BH estimates from indirect methods with the corresponding dynamical ones is obtained by computing their ratio.However, since the same method can yield for some objects estimates that are much larger and for others much smaller than the dynamical ones, the differences may cancel out when computing the average of the ratios for the whole sample.To circumvent this problem, for each source we computed the maximum value between  BH,dyn / BH,indirect and  BH,indirect / BH,dyn (that is, a value that by definition cannot be smaller than 1) and then calculated its mean, ⟨max( BH ratio)⟩, for the dynamical samples.The results of this analysis, reported in the last column of Table 3, confirm our previous conclusions with smaller values of ⟨max( BH ratio)⟩ associated with the X-ray scaling method, the X-ray variability one, and the SE.
A synthesis of the various statistical comparisons of indirect methods with the dynamical ones is illustrated in Fig. 4, where circles indicate indirect methods tested using the restricted dynamical sample and squares refer to methods tested using the extended dynamical sample.In this diagram, the x-axis represents ⟨max( BH ratio)⟩ and the y-axis the probability associated with the K-S test, whereas the auxiliary color-coded axis indicates the Spearman's correlation rank coefficient .Based on these criteria, darker symbols located in the top left corner imply a good agreement between indirect and dynamical methods; the opposite is true of the lighter-colored symbols located in the bottom right corner.
Based on this diagram, we can conclude that the X-ray scaling method shows the best agreement with the dynamical measurements according to the three criteria illustrated in Fig. 4 (the symbol is dark, indicating the positive correlation, and located in the top left corner).Also the SE method (which is only applicable to AGN with broad lines) shows a good agreement, as well as the X-ray variability method.In the latter case, Fig. 4 reveals a clear distinction between the method based on the correlation obtained by Akylas et al. (2022), which shows a moderately good agreement with the dynamical values, and the one obtained in this work, which is more consistent with the dynamical one.Finally, the same figure reinforces the conclusion that neither the - ★ method nor the fundamental plane of BH activity yields  BH measurements consistent with the dynamical values.

Comparison for the whole sample
According to all statistical tests performed so far, the X-ray scaling method appears to yield  BH estimates that are consistently in good agreement with the dynamical values.We will therefore consider this method as an accurate proxy of the dynamical  BH estimate and carry out a similar comparative analysis for the whole sample of AGN.The direct comparison with - ★ , FP, SE, and X-ray excess variance methods is shown in Fig. 5, where the  BH estimates obtained with these four indirect methods are plotted versus the values from the X-ray scaling method.Once more the - ★ method confirms the tendency to systematically overestimate the  BH , whereas FP often yields values that differ by orders of magnitude in both directions from those considered more reliable.The SE method shows a generally good agreement with the exception of a handful of objects whose mass appears to be underestimated by this method.Finally, it is interesting to note the surprisingly good agreement between the two X-ray based methods, which estimate  BH starting from very different assumptions.In Fig. 6 we illustrate the results of this comparison showing the diagram  K−S vs. ⟨max( BH ratio Xscal )⟩, which is defined as the average of the maximum between  BH,Xscal / BH,indirect and  BH,indirect / BH,Xscal .Based on this figure we conclude that, using all the available sources in this volume-limited sample, the X-ray variability method, based on the correlation determined in this work, shows the best agreement with the X-ray scaling method.A reasonably good agreement is achieved also by the SE method and the X-ray variability one based on the Akylas et al. (2022) correlation.Conversely, the - ★ and FP methods appear to be inconsistent with the X-ray scaling method.

DISCUSSION
The primary goal of this work is to compare different indirect methods to determine the  BH in AGN, with particular emphasis on methods that could be applied to heavily obscured objects, where the most frequently used optically based methods are not applicable.To this end, we have used a volume-limited sample of 32 AGN (D < 40 Mpc) selected from the 70-month BAT catalog, which can be considered representative of the AGN population in the local universe, since the detection at very hard X-rays strongly reduces the bias towards unobscured AGN.We have taken advantage of the second data release of the Swift BAT AGN Spectroscopic Survey (BASS DR2), which provides a wealth of useful information, including reliable redshift-independent distances, which play a crucial role in the determination of  BH in some methods, including the direct dynamical ones.
We started from two subsamples of 11 (the restricted dynamical sample) and 17 AGN (the extended dynamical sample) respectively, with  BH determined via dynamical methods to assess the reliability of different indirect methods applicable to all AGN regardless of their level of obscuration, such as the X-ray scaling method, the X-ray variability one based on the normalized excess variance, the - ★ correlation, and the fundamental plane of BH activity.The first important result is that the  BH values derived from X-ray based methods are in good agreement with the dynamical values.On the contrary, the values obtained with the - ★ correlation are systematically and substantially overestimated with respect to the dynamical values, whereas the FP method often yields values that differ from the dynamical ones by one or two order of magnitudes in both directions (see Fig. 3).

X-ray based methods
The bottom right panel of Fig. 5 illustrates the remarkable agreement between the X-ray scaling method and the X-ray variability one.Note that the  BH values derived from the normalized excess variance method are based on a best fit anti-correlation, log  BH = −1.04× log( 2 NXV /0.005) + 7.14, obtained from the extended dynamical sample, which by construction will yield values fully consistent with the dynamical ones.This anti-correlation is slightly different from that obtained by Akylas et al. (2022), based on a large sample of AGN with RM masses, log  BH = −1.13×log( 2 NXV /0.005)+7.38,which still yields  BH in broad agreement with the dynamical values but at a lower statistical level (see Figs. 4 and 6).
These two X-ray methods are based on completely different assumptions and methodologies.More specifically, the X-ray scaling method scales up the dynamical mass of a stellar-mass BH under the assumptions that X-rays in XRBs and AGN are produced by the same Comptonization process and that the photon index Γ is a reliable indicator of the accretion rate in Eddington units.On the other hand, the method based on the  2 NXV is completely model-independent and exploits the universal anti-correlation observed in AGN between  BH and X-ray variability.Also their respective limitations are different: for example, the X-ray scaling method by construction cannot be applied to very low-accreting BHs because in that regime Γ is not positively correlated with  Edd (Jang et al. 2014), but for moderate and highly accreting objects this method constrains  BH with uncertainties ranging between 15 per cent and 40 per cent depending on the quality of X-ray spectra and the reference source utilized (see, e.g.Williams et al. 2023).The X-ray variability method instead has no specific limitations for very low-accreting objects; however, being a statistical method, it yields  BH values with an intrinsic uncertainty of the order of 0.25 dex (the average scatter around the best fit) and can be currently applied to a fairly limited number of AGN since objects with sufficiently long, good-quality hard X-ray light curves are still a minority compared to those with good-quality spectra.Considering that the X-ray variability method by construction must yield values in agreement with the dynamical ones, the full consistency with the X-ray scaling method confirms that  BH estimates derived with the latter method can be considered as reliable surrogates of dynamical values.We therefore utilize the values obtained with the X-ray scaling method as a reference for the comparison with those from other indirect methods over the entire sample.

SE method
Although we are mainly interested in constraining  BH in heavily obscured AGN, we can use the broad-lined AGN contained in our volume-limited sample (14 objects out of 32) and their X-ray scaling  BH estimates to assess the reliability of the SE method, which is the most utilized indirect method to determine  BH in type 1 AGN.This test is illustrated in the bottom left panel of Fig. 5, where the  BH obtained with the SE technique are plotted vs. those from the X-ray scaling method.This figure suggests a reasonably good agreement between the two methods, which is confirmed by the same statistical tests carried out systematically throughout this work:  K−S = 0.54,  t = 0.24,  = 0.48 ( S = 0.08),  = 0.38 ( K = 0.06), and ⟨max( BH ratio)⟩ = 4.8 and are summarized in Fig. 6.
However, the same figure also reveals a handful of objects apparently underestimated by the SE by nearly one order of magnitude with respect to the corresponding X-ray scaling values.Based on our X-ray spectral analysis, all these outliers possess  H along the line of sight in excess of 2×10 22 cm −2 with three of them (NGC 2992, NGC 1365, andNGC 5728) of the order of a few units in 10 23 cm −2 and one (NGC 5506) in the CT regime, suggesting that the discrepancy may be ascribed to the partial obscuration of the BLR that reduces the observed line width leading to  BH that are underestimated.This result appears to be in qualitative agreement with the conclusions of Caglar et al. (2020) and Mejía-Restrepo et al. (2022), who suggest the need of a correction to the  BH values obtained with the SE in heavily absorbed AGN.However, their conclusions were based on the comparison with  BH from the - ★ correlation, under the assumption that the latter yields reliable values, which is at odds with the findings of this work.The bottom left panel of Fig. 5 also shows that other sources that are heavily absorbed yield SE-based  BH that are fully consistent with the X-ray scaling estimates.This may suggest that there is not a simple corrective factor or a monotonic function of  H that can be systematically applied to all SE estimates in heavily obscured AGN.

M-sigma correlation
When good-quality X-ray data exist (that is, high signal-to-noise broadband spectra or light curves with sufficiently long exposures), then X-ray based methods appear to be the most reliable way to constrain  BH .Unfortunately, the current paucity of high-quality X-ray data limits the application of these methods to a few hundreds of AGN.On the other hand, since reliable optical data of host galaxy bulges exist for a much larger number of AGN, it is important to investigate the possibility of applying a correction factor to make the  BH obtained with the - ★ correlation more consistent with those from the X-ray scaling method (and hence with the dynamical ones).The top panel of Fig. 5, showing that the values obtained with the - ★ correlation are systematically overestimated compared to the X-ray scaling estimates for the whole sample, confirms and reinforces the conclusion derived from the comparison with the dynamical  BH estimates.This conclusion appears to be in agreement with the findings of Greene et al. (2016) and Ricci et al. (2017b), who compared the  BH estimates from the - ★ correlation of heavily obscured AGN to the values obtained either from maser measurements or from an IR-based SE method, and found that the - ★ values are substantially overestimated.
The reason why in this work we utilized the - ★ correlation of Kormendy & Ho (2013) is twofold: 1) Koss et al. (2022b) utilize the correlation of Kormendy & Ho (2013) to estimate the  BH for several hundreds of AGN in the BASS DR2.Therefore, determining a suitable correction factor would provide the simplest way to quickly obtain more reliable  BH (and, consequently, the Eddington ratio  Edd ) values for hundreds of AGN. 2) Based on our previous work, where we carried out a systematic comparison of several indirect methods using NGC 4151 as a reference, only the correlation of Kormendy & Ho (2013) yielded an estimate consistent with the  BH obtained from different dynamical methods (Williams et al. 2023).
To investigate whether a different - ★ correlation may provide better estimates of  BH in our volume-limited sample, in Fig. 7 we plotted the  BH obtained with the X-ray scaling method versus the bulge stellar velocity dispersion  ★ and overlaid the - ★ correlation of Kormendy & Ho (2013)  For completeness, we also independently computed the best linear fit with three different methods: we first used ladfit, a robust least absolute deviation method, which does not take into account the errors but is unaffected by outliers, and obtained log  BH M ⊙ = (7.14)+ (2.36) × log  ★ 200 km s −1 .Then, we utilized the linear regression routine linmix_err that accounts for the uncertainties on both variables using a Bayesian approach (Kelly 2007), obtaining the linear regression log  BH M ⊙ = (6.97± 0.03) + (1.79 ± 0.13) × log  ★ 200 km s −1 with a standard deviation of 0.75 on the posterior distribution of the slope values (this fit is represented by the short-dashed magenta line in Fig. 7).Finally, using fitexy, another linear least-squares routine that accounts for the the uncertainties on both x and y, we obtained log  BH M ⊙ = (7.18±0.06)+ (2.66±0.27)×log  ★ 200 km s −1 (represented by the short-dashed orange line in Fig. 7).
Regardless of the routine used for the linear fit, the slope obtained is fairly shallow, which is marginally consistent with the correlation of Woo et al. (2013) but not with that of Kormendy & Ho (2013).In the remainder of the paper, for our independently derived correlation we will use the results from fitexy, which accounts for both uncertainties in x and y, is fully consistent with the results from ladfit, and has a slope slightly steeper than that derived with linmix_err.However, the same conclusions can be derived from the latter linear fit as well.
Importantly, we note that our results, based on a volume-limited unbiased sample (D < 40 Mpc which roughly corresponds to z < 0.01) of type 1 and type 2 AGN selected from the BAT survey, are in full agreement with those recently obtained by Caglar et al. (2023), where an - ★ correlation with a slope of 3.09 ± 0.39 was derived using a larger volume limited sample (z < 0.08) of 154 type 1 AGN from the 105-month BAT sample with reliable  BH derived from the SE method.This may suggest that the - ★ correlations derived from inactive or nearly quiescent galaxies, such that of Kormendy & Ho (2013), are not appropriate to determine the  BH in AGN.
From Fig. 7 it is clear that the Woo et al. (2013) - ★ correlation, as well as the ones derived in this work, provide a better fit of the data while the correlation from Kormendy & Ho (2013) overestimates most of the  BH .Since the discrepancy appears to be related mostly to the intercept, in principle one could still obtain a reasonable estimate of  BH simply renormalizing the Kormendy & Ho (2013) correlation.This process is illustrated in Fig. 8, which shows that the  BH estimates, obtained from the - ★ correlation of Kormendy & Ho (2013) (grey circles), become broadly consistent with the corresponding X-ray scaling values once they are scaled down by a factor of 10 (red circles), and largely overlap the estimates derived from the - ★ correlation of Woo et al. (2013) (illustrated by the blue triangles), and the ones obtained from our best-fitting model (orange open squares).
A more comprehensive comparison of the various - ★ correlations is shown in Fig. 9, where we utilize the same diagnostic diagram introduced in Section 5 to illustrate the statistical comparisons of the different indirect methods.In this figure, M-sigma represents the correlation derived in this work, M-sigma W indicates the correlation of Woo et al. (2013), M-sigma KH the correlation of Kormendy & Ho (2013), and M-sigma KH*0.1 the same correlation scaled down by a factor of 10; finally, M-sigma S represents the correlation of Shankar et al. (2016).Fig. 9 reveals that  BH estimates obtained with the correlation of Kormendy & Ho (2013) scaled down by a factor of 10 (which yields values similar to those obtained from Woo et al. (2013)) have a considerably better agreement with the  BH values from the X-ray scaling method (and hence the dynamical ones) than the original Kormendy & Ho (2013) correlation.This figure also reveals that, based on the ⟨max( BH ratio)⟩ value, which is the most stringent of the three criteria depicted in the diagram, the - ★ correlation derived in this work shows the best agreement with the accepted  BH estimates.

CONCLUSION
Here, we summarize our main results, obtained using a volumelimited sample of AGN selected from the 70-month BAT catalog: • Of all indirect methods used to estimate  BH in obscured AGN, which were tested in this work, only the X-ray based ones (the scaling method and the variability one based on the normalized excess variance) yield values that are consistently in agreement with the dynamical estimates.For this reason, in the absence of a more direct dynamical way to constrain  BH and when good-quality X-ray data with sufficiently long exposure are available, X-ray based estimates should be considered as the most reliable and used as a reference to calibrate other indirect methods.
• BH values obtained with the - ★ correlation proposed by Kormendy & Ho (2013), which was also used to estimate the vast majority of  BH in the BASS DR2 (Koss et al. 2022b), appear to be systematically overestimated when compared to dynamical measurements (and hence to X-ray based ones).Scaling down these estimates by a factor of 10 substantially improves the agreement with the corresponding dynamical values (and X-ray based estimates) and represents the fastest way to obtain reasonable  BH for hundreds of AGN.
• In general, AGN  BH values appear to be more consistent with the estimates derived from the - ★ correlation obtained by Woo et al. (2013) for an RM AGN sample or the one derived in this work, log  BH M ⊙ = (7.18± 0.06) + (2.66 ± 0.27) × log  ★ 200 km s −1 .These correlations with shallower slopes are consistent with the one recently derived by Caglar et al. (2023) using a larger unbiased volume-limited sample of type 1 AGN.This may suggest that AGN follow a different - ★ correlation from the one inferred using inactive or nearly quiescent nearby galaxies.
• For the broad-lined AGN of our sample, a statistical comparison between SE and X-ray scaling estimates of  BH indicates a good agreement, but also confirms the tendency of the SE technique to underestimate the  BH in sources affected by substantial obscuration.It is therefore important to properly correct the SE values of heavily obscured sources, as suggested by Mejía-Restrepo et al. (2022).However, the corrections should not be based on the comparison with the - ★ correlation from inactive galaxies, since the latter has a tendency to overestimate the  BH .Our analysis also reveals that some heavily absorbed sources yield SE-based  BH that are fully consistent with the X-ray scaling estimates, suggesting that there is not a simple corrective factor that can be uniformly applied to all SE estimates in heavily obscured AGN.
• As expected from previous studies, the  BH estimates obtained with the FP are fairly inconsistent and unreliable, with values that can differ from the dynamical ones by two orders of magnitude in both directions and with only a small fraction of objects with  BH consistent with dynamical values (and X-ray based ones).Therefore, one should use considerable caution when deriving conclusions based on  BH estimates derived from this method, which should not be used in isolation.In the Analysis column, spec denotes objects we conducted spectral analysis on, and for those we calculated masses with the X-ray scaling method.In the same column, var denotes objects we conducted temporal analysis on, and for those we calculated masses with the X-ray variability method.

Figure 1 .
Figure 1. BH plotted vs.  2 NXV in the 10-20 keV band.The symbols are color coded based on the intrinsic  H along the line of sight; circles indicate objects of the dynamical restricted sample, whereas squares are the added objects in the extended sample.Compton-thick (CT) sources (NGC 1068, NGC 4945, Circinus) have been excluded from the fit but are shown in the plot.The green continuous line represents the best fit, log  BH = −1.04× log(  2 NXV /0.005) + 7.14, obtained from the extended dynamical sample without the changing-look AGN NGC 2992, whose low variability may be affected by its peculiar state and hence unreliable.The short-dashed black line represents the best fit, log  BH = −1.13× log(  2 NXV /0.005) + 7.38, from Akylas et al. (2022) using reverberation mapping data only.For completeness, we also show the best fit, log  BH = −0.53× log(  2 NXV /0.005) + 7.23, obtained using only the restricted dynamical sample (long-dashed red line), which is limited to six objects (after excluding CT AGN and objects without reliable  2 NXV ), and the best fit, log  BH = −0.29 × log(  2 NXV /0.005) + 6.97, obtained using the extended dynamical sample (ten objects; blue longdashed line).

Figure 2 .
Figure 2. Top panels: comparison between the restricted dynamical  BH distribution and those obtained with the - ★ correlation (left), with the fundamental plane of black hole activity (middle), and with the X-ray scaling method (right).Bottom panels: comparison between the extended dynamical  BH distribution and the ones obtained with single epoch H α measurements (left), with the X-ray variability method based on the correlation derived by Akylas et al. (2022) using the RM sample (middle), and the variability method based on the correlation derived in this work (right).

Figure 3 .
Figure 3. Top left:  BH obtained from the - ★ correlation plotted vs. those obtained from dynamical methods (circles indicates objects from the restricted sample, whereas squares represent the sources added to form the extended dynamical sample).Top middle:  BH from the fundamental plane of black hole activity vs. dynamical values.Top right: values from the X-ray scaling method vs. the dynamical ones.Bottom left:  BH obtained with single epoch H α measurements plotted vs. those obtained from dynamical methods.Bottom middle and right:  BH obtained from X-ray variability using the correlation of Akylas et al. (2022) and that obtained in this work, respectively, plotted vs. the dynamical measurements.The symbols are color coded based on the intrinsic  H along the line of sight.

Figure 4 .
Figure 4.The quantity along the x-axis is ⟨max(  BH ratio) ⟩, the average of the maximum between  BH,dyn / BH,indirect and  BH,indirect / BH,dyn , computed for each source: the smaller the x value the better the agreement with the dynamical value.The quantity along the y-axis is the probability associated with the K-S test: higher values indicate that the distribution of the  BH obtained with an indirect method is consistent with that of dynamical values.Finally, the auxiliary color-coded axis represents the Spearman's rank correlation coefficient : darker colors indicate a positive correlation, light colors indicate a negative or nonexistent linear correlation between  BH indirect values and the corresponding dynamical ones.Circles represent methods that were applied to the restricted dynamical sample, whereas squares refer to the indirect methods applied to the extended dynamical sample.Dark symbols in the top left corner indicate the best agreement between indirect and dynamical methods; the opposite is true for the bottom right corner and lighter colors.

Fig. 3
Fig.3also confirms that a substantial number of  BH values derived with the - ★ correlation method are overestimated compared to the corresponding dynamical ones (top left panel) and indicates that in some cases the FP estimates may exceed the dynamical values by orders of magnitude and in others yield significantly lower values (top middle panel).A quantitative way to assess the agreement (or lack thereof) illustrated in Fig.3is based on the Spearman's and Kendall's rank correlation coefficients  and  and their associated probabilities, with a positive correlation expected for the indirect methods that yield  BH estimates consistent with the dynamical ones.The results of these statistical tests are reported in the fifth and sixth columns of Table3and corroborate the previous conclusions with negative correlation coefficients obtained for the - ★ and FP methods and positive correlation coefficients for the SE and the two X-ray based methods.An additional straightforward way to directly compare the  BH estimates from indirect methods with the corresponding dynamical ones is obtained by computing their ratio.However, since the same method can yield for some objects estimates that are much larger and for others much smaller than the dynamical ones, the differences may cancel out when computing the average of the ratios for the whole sample.To circumvent this problem, for each source we computed the maximum value between  BH,dyn / BH,indirect and  BH,indirect / BH,dyn (that is, a value that by definition cannot be smaller than 1) and then calculated its mean, ⟨max( BH ratio)⟩, for the dynamical samples.The results of this analysis, reported in the last column of Table3, confirm our previous conclusions with smaller values of ⟨max( BH ratio)⟩ associated with the X-ray scaling method, the X-ray variability one, and the SE.A synthesis of the various statistical comparisons of indirect methods with the dynamical ones is illustrated in Fig.4, where circles indicate indirect methods tested using the restricted dynamical sample and squares refer to methods tested using the extended dynamical sample.In this diagram, the x-axis represents ⟨max( BH ratio)⟩ and the y-axis the probability associated with the K-S test, whereas the auxiliary color-coded axis indicates the Spearman's correlation rank coefficient .Based on these criteria, darker symbols located in the top left corner imply a good agreement between indirect and dynamical methods; the opposite is true of the lighter-colored symbols located in the bottom right corner.Based on this diagram, we can conclude that the X-ray scaling method shows the best agreement with the dynamical measurements according to the three criteria illustrated in Fig.4(the symbol is dark, indicating the positive correlation, and located in the top left

Figure 5 .
Figure 5. Top left panel:  BH obtained with the - ★ correlation plotted vs. those from the X-ray scaling method.Top right panel:  BH values from the fundamental plane of black hole activity vs. those from the X-ray scaling method.Bottom left panel:  BH obtained with the SE method plotted vs. the values derived from the X-ray scaling method.Bottom right panel:  BH from X-ray variability of 80 ks segments vs. those from the X-ray scaling method.The symbols are color coded based on the intrinsic  H along the line of sight.

Figure 6 .Figure 7 .
Figure 6.Probability associated with the K-S test plotted vs. ⟨max(  BH ratio Xscal ) ⟩, the average of the maximum between  BH,Xscal / BH,indirect and  BH,indirect / BH,Xscal .The auxiliary colorcoded axis represents the Spearman's rank coefficient .Dark symbols in the top left corner indicate good agreement between X-ray scaling methods and the other indirect methods; the opposite is true for the bottom right corner.
(black continuous line), as well as the correlation of Woo et al. (2013) (the blue long-dashed line).

Figure 8 .
Figure 8.  BH obtained with different - ★ correlations plotted vs. those from the X-ray scaling method.The grey circles indicate the values obtained from the original correlation from Kormendy & Ho (2013), whereas the red circles show the same values scaled down by a factor of 10.The blue triangles represent the  BH derived from the correlation of Woo et al. (2013), and the open orange squares are the values obtained from the - ★ correlation derived in this work.

Figure 9 .
Figure 9. Probability associated with the K-S test plotted vs. ⟨max(  BH ratio) ⟩, the average of the maximum between  BH,X−scal / BH,indirect and  BH,indirect / BH,X−scal using the whole sample.The auxiliary color-coded axis represents the Spearman's rank coefficient .In this plot we show the different locations in the diagram of the various - ★ correlations, along with the location of the X-ray variability and FP values.M-sigma W indicates the correlation of Woo et al. (2013), M-sigma KH the correlation of Kormendy & Ho (2013), and M-sigma KH*0.1 the same correlation scaled down by a factor of 10; M-sigma S represents the correlation of Shankar et al. (2016); finally, M-sigma is the one obtained in this work.

Table 1 .
Sample properties

Table 3 .
Statistical comparisons of indirect methods = dynamical sample (R means restricted and E extended) used for the comparison with the number of sources available in parentheses, 3 = probability of the Kolmogorov-Smirnov test (large values indicate that the  BH distribution obtained with the indirect method is indistinguishable from that obtained from direct dynamical methods), 4 = probability of the Student's  test (small values indicate that the mean of the distributions compared are significantly different), 5 = Spearman's rank correlation coefficient  and relative probability, 6 = Kendall's rank correlation coefficient  and relative probability, 7 = average of the maximum between  BH,dyn / BH,indirect and  BH,indirect / BH,dyn with the standard deviation divided by the square root of the number of objects in parentheses.