BASS. XXX. Distribution Functions of DR2 Eddington Ratios, Black Hole Masses, and X-Ray Luminosities

2.1.1. Input Catalog and Sample

2.1.2. Excluded Sources

2.1.3. Final, Redshift-restricted AGN Sample

After applying all the cuts mentioned so far, our base sample of nonbeamed, non-Galactic-plane AGNs includes 678 sources, and is used for (parts of) our XLF analysis. We next introduce a number of additional criteria to retain only those AGNs that lie in regions of parameter space in which the BAT and BASS selection functions are well understood, highly complete, and highly reproducible.

We first restrict our analysis to sources in the 0.01 ≤ z ≤ 0.3 range. This excludes the 53 nearest BASS/DR2 AGNs, for which redshift-based distance determinations may be affected by peculiar velocities, and/or which may be outliers in terms of their location in the L–z plane. Altogether, the 0.01 ≤ z ≤ 0.3 sample contains 619 objects (after excluding six more AGNs at z > 0.3).

We next limit our sample to have reliable measurements of L_bol, M_BH, and λ_E, within ranges that are reasonable for persistent, radiatively efficient accretion onto SMBHs, and that can be probed within BASS/DR2 with a high degree of completeness. These basic measurements are described in Section 2.2 below, while the chosen ranges are listed in Table 1. Specifically, we include only those AGNs with BH masses in the range $6.5\leqslant \mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\leqslant 10.5$ and with Eddington ratios in the range $-3\leqslant \mathrm{log}{\lambda }_{{\rm{E}}}\leqslant 1$.

Table 1. Overview of Sample Selection and Analysis Parameters ^a

Quantity/Variable	Symbol/Value
Minimum redshift considered	${z}_{\min ,{\rm{s}}}=0.01$
Maximum redshift considered	${z}_{\max ,{\rm{s}}}=0.3$
Minimum BH mass considered	$\mathrm{log}({M}_{\mathrm{BH},\min ,{\rm{s}}}/{M}_{\odot })=6.5$
Maximum BH mass considered	$\mathrm{log}({M}_{\mathrm{BH},\max ,{\rm{s}}}/{M}_{\odot })=10.5$
Minimum Eddington ratio considered	$\mathrm{log}{\lambda }_{{\rm{E}},\min ,{\rm{s}}}=-3$
Maximum Eddington ratio considered	$\mathrm{log}{\lambda }_{{\rm{E}},\max ,{\rm{s}}}=1$
Luminosity bin size for Vmax method	$d\mathrm{log}{L}_{{\rm{X}}}=0.3$
BH mass bin size for Vmax method	$d\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })=0.3$
Edd. ratio bin size for Vmax method	$d\mathrm{log}{\lambda }_{{\rm{E}}}=0.3$
Assumed uncertainty on $\mathrm{log}{M}_{\mathrm{BH}}$	${\sigma }_{{M}_{\mathrm{BH}}}=$ 0.3, 0.5 dex
Assumed uncertainty on $\mathrm{log}{\lambda }_{{\rm{E}}}$	${\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}=$ 0.3, b 0.5 dex
Galactic plane exclusion c	∣b∣ ≥ 5°
Other cuts	Beamed AGNs are excluded

Notes.

^a No explicit flux limits were imposed during our XLF/BHMF/ERDF analysis; instead, we used the full flux–area curve of the 70 month Swift/BAT survey (Baumgartner et al. 2013; see the main text). ^b We also investigate the effect of an uncertainty of 0.2 dex in luminosity, due to measurement uncertainty and the X-ray bolometric correction (Equation (2)). Variable bolometric corrections are explored in Appendix E. ^c Meant to guarantee a high completeness of optical counterpart identification and spectroscopic coverage.

Download table as:  ASCIITypeset image

Throughout our analysis, we further classify sources as being either broad-line ("Type 1," hereafter) or narrow-line ("Type 2," hereafter) AGNs, based on the presence of any broad Balmer lines (as described in Koss et al. 2022a, 2022b), and on how their most reliable measures of M_BH and λ_E were derived (see Section 2.2 below), which in turn has consequences for uncertainty and incompleteness estimations. Thus, our sample of Type 1 AGNs (366 sources)—sources that have at least one broad Balmer emission line—also includes so-called "Seyfert 1.9" AGNs, with broad Hα but no broad Hβ. Our Type 2 AGNs (220 sources) have only narrow Balmer emission lines. Type 1 AGNs are relatively unobscured [e.g., $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\leqslant 22]$, whereas Type 2 AGNs tend to be more heavily obscured $[\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\gt 22$]; although, as shown in Section 3, these $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})$ limits do not apply strictly (see also Oh et al. 2022).

For our BHMF and ERDF analysis, we exclude the four Seyfert 1.9 sources that have mass estimates only from broad Hα lines, as such mass estimates were shown to be highly uncertain for heavily obscured Seyfert 1.9s in the companion BASS/DR2 paper by Mejía-Restrepo et al. (2022). One of these four sources (ID 476) is the only AGN in our sample that has an estimated Eddington ratio formally greater than 10 ($\mathrm{log}{\lambda }_{{\rm{E}}}\simeq 1.4$). Therefore, the λ_E ratio upper limit we impose does not exclude objects that would otherwise have been included in the analysis. Similarly, there are no sources in the BASS/DR2 sample with BH masses greater than the upper limit on M_BH. These two upper limits are important for computational reasons, and are discussed in Section 3.4.

The Venn diagrams in Figure 1 show the number of nonbeamed BASS/DR2 AGNs that meet each of our selection criteria (in z, M_BH, and/or λ_E), as specified in Table 1, for both Type 1 and Type 2 AGNs. We also show the number of sources that fall outside all these criteria. Our main BHMF and ERDF analysis is done using the 586 AGNs that meet all criteria (366 Type 1 and 220 Type 2 AGNs), as shown in the central regions of the Venn diagrams.

Zoom In Zoom Out Reset image size

The present analysis is based on four basic, interlinked measurements that are available for the BASS/DR2 AGNs thanks to the extensive X-ray and optical spectroscopy that is at the heart of the BASS project: bolometric luminosities (L_bol), BH masses (M_BH), Eddington ratios (λ_E), and line-of-sight hydrogen column densities (N_H). The determination of these quantities from basic observables is described in detail in other BASS publications (see below). Here, we provide only a brief summary of the aspects most relevant for the present analysis. Two key considerations for deriving these quantities are to adopt a uniform approach for all BASS AGNs, whenever possible (e.g., for L_bol), and to revert to differential methodologies only when absolutely needed (e.g., for M_BH determination).

First, as a primary probe of AGN luminosity, we use the integrated, intrinsic luminosity between 14–195 keV (L_{14–195 kev} or L_X, hereafter), as determined through an elaborate spectral fitting of all the available X-ray data for each BASS source, as described in detail in Ricci et al. (2017a). This X-ray spectral decomposition also yields the N_H measurements we use here, and we stress that the intrinsic L_{14–195 kev} already accounts for the line-of-sight obscuration (as probed by N_H). The measurement uncertainties on the X-ray luminosities we use are rather small, not exceeding ∼0.05–0.1 dex for unobscured sources. Thus, whenever uncertainties on luminosities are invoked throughout our analysis, these relate to bolometric luminosities (unless otherwise noted), and are dominated by the systematics on the X-ray to bolometric luminosity conversion (see below).

Bolometric luminosities, L_bol, are then derived directly from L_{14–195 kev}, by using a simple, universal scaling of

$$\begin{eqnarray}&&{L}_{\mathrm{bol}}={\kappa }_{14-195\,\mathrm{keV}}\times {L}_{14-195\,\mathrm{keV}}.\end{eqnarray} \tag{ 2 }$$

As explained in other key BASS publications (e.g., Koss et al. 2022b), this bolometric correction of κ_{14−195 keV} = 7.4 corresponds to a lower-energy bolometric correction of κ_{2−10 keV} = 20, which is the average correction found for the BASS sample following the L_{2–10 kev}-dependent prescription of Marconi et al. (2004), and further assumes the median X-ray power-law index found for the BASS sample, Γ = 1.8 (e.g., Lanzuisi et al. 2013). This bolometric correction is also in agreement with Vasudevan et al. (2009). Our choice to use L_{14–195 kev}, and not other (X-ray) luminosity probes, is motivated by (1) the need to work as closely as possible with the Swift/BAT selection functions, (2) our desire to have the most reliable determinations of L_bol for even the most obscured AGNs [i.e., Compton-thick sources with $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\gtrsim 24$], and (3) our desire to be consistent with previous studies of the XLF of Swift/BAT-selected AGNs (e.g., Sazonov et al. 2007; Tueller et al. 2008; Ajello et al. 2012). We acknowledge that several other, higher-order bolometric correction prescriptions have been suggested and used in the literature, including corrections that depend on luminosity, λ_E, and/or other AGN properties (e.g., Marconi et al. 2004; Vasudevan & Fabian 2007, 2009; Jin et al. 2012; Lusso et al. 2012; Brightman et al. 2017; Netzer 2019; Duras et al. 2020, and references therein). In the main part of the text, we prefer to use a constant bolometric correction to simplify our already complicated decomposition of the XLF, BHMF, and ERDF, and to be consistent with the rest of the BASS (DR2) analyses. However, we report how our conclusions change with a luminosity-dependent bolometric correction in Appendix E.

BH masses, M_BH, are determined using two different approaches for AGNs specifically with broad Hα line emission, and with or without broad Balmer-line emission, which allows for a certain M_BH estimation procedure (see immediately below). As noted above, for simplicity we refer to such sources as Type 1 and Type 2 sources, respectively, but we note that their detailed classification may be more nuanced (see Koss et al. 2022a; Mejía-Restrepo et al. 2022).

For broad-line (Type 1) AGNs, we rely on detailed spectral decomposition of the Hα spectral complex and "virial" BH mass prescriptions that are calibrated against reverberation mapping experiments, as described in detail in Mejía-Restrepo et al. (2022). Specifically, Mejía-Restrepo et al. followed the prescription provided by Greene & Ho (2005, Equation (6)), but adjust the virial factor to f = 1, yielding:

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=2.67\times {10}^{6}{\left(\displaystyle \frac{L[\mathrm{bH}\alpha ]}{{10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)}^{0.55}{\left(\displaystyle \frac{\mathrm{FWHM}[\mathrm{bH}\alpha ]}{1000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{2.06}{M}_{\odot },\end{eqnarray} \tag{ 3 }$$

where L[bHα] and FWHM[bHα] are the luminosity and width of the broad part of the Hα emission line, respectively. For narrow-line (Type 2) AGNs, Koss et al. (2017) and Koss et al. (2022a) rely on measurements of the stellar velocity dispersion (σ_⋆) in the AGN host galaxies and the well-known M_BH − σ_⋆ relation. Specifically, σ_⋆ was measured from the Ca H + K, Mg i, and/or the Ca ii triplet spectral complexes, as described in detail in Koss et al. (2022c). We then used the relation determined by Kormendy & Ho (2013, Equation (3)):

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=3.09\times {10}^{8}{\left(\displaystyle \frac{{\sigma }_{\star }}{200\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{4.38}{M}_{\odot }.\end{eqnarray} \tag{ 4 }$$

We stress that these two types of M_BH prescriptions are consistently calibrated. Specifically, the broad-line M_BH prescription in Equation (3) is derived by assuming that low-redshift, broad-line AGNs (particularly those with reverberation mapping measurements) lie on the same M_BH − σ_⋆ relation as do narrow-line AGNs and inactive galaxies (that is, Equation (4); see, e.g., Onken et al. 2004 and Woo et al. 2013 for detailed discussions).

The uncertainties in both types of M_BH estimates are completely dominated by systematics and are of the order 0.3–0.5 dex (Gültekin et al. 2009; Shen 2013; Peterson 2014; Shankar et al. 2019). The lower end of this range is consistent with the scatter seen in the M_BH − σ_⋆ (or M_BH − M_host) relation (e.g., Sani et al. 2011; Kormendy & Ho 2013), while the upper end also includes the other key ingredients of "virial" M_BH determinations in broad-line AGNs (e.g., Shen 2013; Peterson 2014). The uncertainties on λ_E are naturally dominated by the systematic uncertainties on M_BH, and are thus also of the order 0.3–0.5 dex, although the ranges and trends seen for bolometric corrections (e.g., Marconi et al. 2004; Vasudevan & Fabian 2009) suggest that the uncertainties on λ_E may be yet higher. In comparison, our measurement uncertainties on L(bHα), FWHM(bHα), and σ_⋆ are typically of order 10%, which would add up to ∼0.1 dex uncertainty for M_BH estimates in broad-line AGNs and <0.2 dex in narrow-line AGNs. Our analysis takes into account the large (systematic) uncertainties on M_BH determinations for all BASS AGNs, as detailed in Section 3.

Finally, λ_E is determined by combining the L_bol and M_BH estimates mentioned above, using the relation

$$\begin{eqnarray}&&\mathrm{log}{\lambda }_{{\rm{E}}}=\mathrm{log}({L}_{\mathrm{bol}}/\mathrm{erg}\,{{\rm{s}}}^{-1})-\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })-38.18,\end{eqnarray} \tag{ 5 }$$

which is appropriate for solar metallicity gas. As the uncertainties in both luminosity and mass must be taken into account to calculate λ_E, we present one case in our main analysis where we assume a higher error in $\mathrm{log}\,{\lambda }_{{\rm{E}}}$ to take the uncertainties in luminosity measurement/bolometric correction into account, and find that our results converge on the same solution. Larger uncertainties in $\mathrm{log}\,{\lambda }_{{\rm{E}}}$ and results due to variable bolometric correction are presented in Appendix E.

There are alternative galaxy–BH scaling relationships suggested by Bernardi et al. (2007), Shankar et al. (2017), and Shankar et al. (2020), which correct for the selection biases that may affect the samples used for calibrating these mass measurement relationships. While recalculating the masses computed using Equation (4) (i.e., velocity dispersion) would be trivial, consistently recalibrating all the other masses, calculated using various other methods, to account for these selection effects is beyond the scope of the present analysis.

Figure 2 shows the cumulative and differential source counts (the number density per square degree and its differential form, respectively) for the various samples relevant for the present study, specifically: (1) the input 70 month catalog (z ≤ 0.3;672 sources); (2) our redshift-restricted AGN sample (0.01 ≤z ≤ 0.3; 619 sources); and (3) our final BASS/DR2 BHMF/ERDF sample (i.e., with M_BH and λ_E cuts; 586 sources). All three samples exclude the sources described in Section 2.1.2. In all cases, the uncertainties are derived assuming that the source counts follow Poisson distributions. Figure 2 also shows the best-fit curves corresponding to the various samples, which use the following functional form for the differential number counts of the three samples:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dS}}=A\times {\left(S/{10}^{-11}\right)}^{-\alpha }.\end{eqnarray} \tag{ 6 }$$

We limit our fits to the F_{14−195, obs} ≥ 10^−11.1 erg s⁻¹ cm⁻² flux regime, since the flux–area curve for the Swift/BAT 70 month survey (discussed in more detail in Section 3) is sparsely sampled at lower flux levels, making it difficult to accurately calculate the surface density of faint BAT sources (see the left panel of Figure 2). Our fits are derived using orthogonal distance regression (ODR), so as to properly account for uncertainties on both axes (e.g., $\mathrm{log}S$ and $\mathrm{log}N$). One limitation of the ODR method is the inherent assumption of symmetric uncertainties (which, in our case, are averages of the upper and lower errors). For the Poisson uncertainties relevant for our analysis, this introduces only a minor change compared with the real, asymmetric uncertainties (e.g., for a bin with 25 objects, the errors are +6.07 and −4.97; see Gehrels 1986). To verify that our results are not significantly affected by the choice of fitting method (and the associated treatment of uncertainties), we have carried out an additional maximum likelihood estimator (MLE) fitting, using a Fechner distribution (see Wallis 2014 and references therein).

Zoom In Zoom Out Reset image size

The best-fit slopes derived using ODR are: α = 2.50 ±0.04, 2.62 ± 0.07, and 2.67 ± 0.05, for the 70 month AGN, the redshift-restricted, and the final BHMF/ERDF AGN samples, respectively. The best-fit normalizations are $\mathrm{log}A=9.43\pm 0.03$, 9.44 ± 0.03, and 9.43 ± 0.02, respectively. The MLE-based best-fit results are in agreement with the ORD ones, within ±1σ errors. Our best-fit curve for the full sample is consistent with the expected slope for a fully uniform Euclidean distribution (α = 2.5), and also with the results of several previous studies of ultra-hard X-ray-selected AGNs (Tueller et al. 2008; Cusumano et al. 2010; Krivonos et al. 2010; Ajello et al. 2012; Harrison et al. 2016).

For the two limited samples, the slope is not perfectly Euclidean, but that is to be expected, as objects have been removed from the total sample. We conclude that our BASS/DR2 AGN samples, which are used to determine the XLF, BHMF, and ERDF, do not show any significant biases compared to the all-sky Swift/BAT 70 month catalog and/or other samples of ultra-hard X-ray-selected AGNs. We are thus confident that we can rely on the same selection criteria, and specifically the sky coverage curves, as derived for the parent Swift/BAT 70 month catalog.

We take advantage of the well-constrained flux–area curve of the Swift/BAT survey (Baumgartner et al. 2013). The flux–area curve ${{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{F}_{{\rm{X}}})$, shown in Figure 3, accounts for the fact that the effective area of the BAT survey is a function of the observed 14–195 keV X-ray flux. For bright, ultra-hard X-ray sources with $\mathrm{log}({F}_{14-195,\mathrm{obs}}/\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1})\gt -10.7$, the BAT survey is complete over the entire sky, i.e., Ω_sel = 1. Sources with $\mathrm{log}({{\rm{F}}}_{{\rm{X}},\ \mathrm{obs}}/\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1})\lt -11.3$ are completely missed by the BAT 70 month survey and catalog, and thus Ω_sel = 0. Our analysis makes explicit use of the complete flux–area curve, including the intermediate values for sources with $-11.3\lesssim \mathrm{log}({F}_{14-195,\mathrm{obs}}/\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1})\lesssim -10.7$. Note that this is the cumulative flux–area curve over the entire sky, and the sensitivity is somewhat nonuniform (as shown in Figure 1 of Baumgartner et al. 2013). Taking the slight positional differences in sensitivity into account is beyond the scope of this work.

Zoom In Zoom Out Reset image size

Additionally, we account for the fact that the XLF describes the distribution of AGNs according to their intrinsic ultra-hard X-ray luminosity, while their detectability is a function of their observed (ultra-hard) X-ray flux, which depends on the amount of obscuration along the line of sight (in addition to distance, obviously). To correct the intrinsic luminosities for obscuration, and to compute the observed ultra-hard X-ray luminosities, we use the Swift/BAT attenuation curve shown in Figure 4. This observed luminosity is a sum of the transmission of the intrinsic power-law radiation and reflection from the accretion disk through the obscuring torus. This 14−195 keV attenuation curve is similar to the attenuation curve of Ricci et al. (2015), which was calculated using the torus model of Brightman & Nandra (2011) and based on the spectral models used in Ueda et al. (2014). We recalculated the attenuation curve using the borus02 model (Baloković et al. 2018), which updated the Brightman & Nandra (2011) torus, further assuming a photon index of Γ = 1.8, E_cutoff = 200 keV, a torus opening angle of 60°, and an inclination angle of 72° (Gilli et al. 2007; Lanzuisi et al. 2013; Ueda et al. 2014; Aird et al. 2015; Ananna et al. 2019, 2020a). We have tested how our results vary with changes in the torus opening angle and the (line-of-sight) inclination angle. Following Ricci et al. (2015), we also consider another attenuation curve for a model with an opening angle of 35° and an inclination angle of 72° (shown in Figure 4). We discuss the impacts of varying our template spectra (and other model-dependent parts of our approach) on our final results in Section 3.4.2.

Zoom In Zoom Out Reset image size

While the BAT 14−195 keV energy range is not strongly affected by attenuation up to $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\simeq 23.5$, Figure 4 shows that for Compton-thick sources $[\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\gtrsim 24]$, at least 30% and up to ∼98% of the intrinsic ultra-hard X-ray emission is lost [i.e., only the scattered component is detectable as $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})$ approaches 26], and the luminosities of such sources may have been drastically underestimated. Our analysis makes explicit use of the complete attenuation curve, thus properly linking (limiting) observed fluxes and intrinsic luminosities, given the measured N_H of each AGN (Ricci et al. 2017a, Koss et al. 2022a).

We note that these links between observed fluxes, intrinsic luminosities, and surveyed angles (and volumes) affect not only the XLF, but also the BHMF and the ERDF, as the three key properties (L_{14–195 kev}, M_BH, and λ_E) are closely related, through Equations (2) and (5).

The 1/V_max method (Schmidt 1968) provides a way to correct for the incompleteness (i.e., lower-luminosity/mass objects falling below the survey sensitivity at larger distances). When counting the sources within a given luminosity (or mass, etc.) bin, each source i is weighted by ${V}_{\max ,{\rm{i}}}$—the maximum volume within which source i could have been detected, given the survey properties. In the case of an LF, low- and high-luminosity objects are weighted by small and large volumes, respectively. This increases and decreases their relative contribution to the respective space densities.

To measure the distribution of the AGN luminosities, BH masses, and Eddington ratios of the BASS sample, we bin in $\mathrm{log}{M}_{\mathrm{BH}}$, $\mathrm{log}{\lambda }_{{\rm{E}}}$, and $\mathrm{log}{L}_{{\rm{X}}}$. ²⁸ We then use the 1/V_max method to determine the corresponding space densities ${{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})$, $\xi (\mathrm{log}{\lambda }_{{\rm{E}}})$, and ${{\rm{\Phi }}}_{{\rm{L}}}(\mathrm{log}{L}_{{\rm{X}}})$. We compute the V_max values corresponding to the intrinsic, ultra-hard X-ray luminosity of each source by considering the respective observed flux and the survey completeness, as detailed below. As this 1/V_max-based calculation does not include a robust correction either for sample truncation or for the uncertainties on the relevant key quantities, it only allows us to gain initial guesses for the BHMF and the ERDF.

For the XLF, the space density in the luminosity bin j is given by the sum over all N_bin-weighted objects within this bin:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{L}},\ {\rm{j}}}\,d\mathrm{log}{L}_{{\rm{X}}}=\sum _{i}^{{N}_{\mathrm{bin}}}\displaystyle \frac{1}{{V}_{\max ,i}}.\end{eqnarray} \tag{ 7 }$$

We compute ${{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})$ and $\xi (\mathrm{log}{\lambda }_{{\rm{E}}})$ similarly, and use the bin sizes $d\mathrm{log}{L}_{{\rm{X}}}$, $d\mathrm{log}{M}_{\mathrm{BH}}$, and $d\mathrm{log}{\lambda }_{{\rm{E}}}$, given in Table 1.

To compute the V_max values for each AGN, we express the completeness of the BAT survey as a function of the intrinsic, ultra-hard X-ray luminosity, redshift, and column density: ${{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{L}_{{\rm{X}}},{N}_{{\rm{H}}},z)$. For each object i at redshift z_i with intrinsic luminosity $\mathrm{log}{L}_{{\rm{X}},i}$ and column density N_H,i , the maximum comoving volume ${V}_{\max ,i}$ is then given by (e.g., Hogg 1999; Hiroi et al. 2012):

$$\begin{eqnarray}&&\begin{array}{l}{V}_{\max ,i}(\mathrm{log}{L}_{{\rm{X}},i},{N}_{{\rm{H}},i},{z}_{i})\\ =\ \ 4\pi \times \displaystyle \frac{c}{{H}_{0}}\times {\displaystyle \int }_{{z}_{\min ,{\rm{s}}}}^{{z}_{\max ,{\rm{s}}}}{{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{L}_{{\rm{X}},i},{N}_{{\rm{H}},i},z)\displaystyle \frac{{D}_{C}{\left(z\right)}^{2}}{E(z)}{dz}.\end{array}\end{eqnarray} \tag{ 8 }$$

In this expression, D_C (z) corresponds to the comoving distance to redshift z, and E(z) is given by $\sqrt{{{\rm{\Omega }}}_{{\rm{M}}}{(1+z)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}$. The integration limits ${z}_{\min ,{\rm{s}}}$ and ${z}_{\max ,{\rm{s}}}$ correspond to the minimum and maximum redshifts of our BASS/DR2 (refined) sample, respectively (see Table 1). Note that, unlike what is done in some studies, our calculation does not explicitly introduce a rigid flux limit, which in turn would impose a maximal redshift up to which each source could be observed, ${z}_{\max ,i}$. Instead, this information is encoded in the flux–area curve, which ultimately provides the same outcome as Ω_sel = 0 when the source's observed flux becomes too faint [i.e., $\mathrm{log}({F}_{14-195,\mathrm{obs}}/\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1})\lesssim -11.3$].

Note that the 1/V_max method is very sensitive to the flux–area curve, which is nonuniform (as described in Section 3.1). As a specific example, the AGN in NGC 5283 (BAT ID 684), with z = 0.01036 and $\mathrm{log}({F}_{14-195,\mathrm{obs}}/\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1})=-11.14$, falls within the criteria used to select the redshift-restricted sample (Table 1), and it is close both to the low-z cut we use and to the regime where the flux–area curve drops to zero. As a result, the V_max within which this object can be detected above z > 0.1 is very low, and this drives up the corresponding 1/V_max value and any related contribution to the population-wide distributions under study. However, the method used to calculate the bias-corrected intrinsic BHMF and ERDF (described in Section 3.4) is not as dramatically affected by a single object as the direct 1/V_max approach. We show the effect of including this single object on the 1/V_max values in Appendix F, whereas the 1/V_max values shown in the plots in the main body of the text exclude this object. We include this object in the bias-corrected part of our analysis.

There are several published statistical methods that address (some of) the limitations of the 1/V_max approach, such as the Lynden-Bell–Woodroofe–Wang estimator (Lynden-Bell 1971; Woodroofe 1985; Choloniewski 1987; Wang 1989; Efron & Petrosian 1992), which corrects for sample truncation. Unlike the 1/V_max method, this estimator does not depend on the assumed bin size and does not assume a constant space density over every given bin. Another flexible Bayesian parametric framework for estimating LFs while correcting for sample truncation has been suggested by Kelly et al. (2008). However, most previous (X)LF works have used the 1/V_max approach. Motivated by the desire for our XLF results to be directly comparable to previous studies, and by the fact that our core analysis methodology does ultimately correct for both sample truncation and the uncertainties on key quantities (in all three distribution functions), we chose to present the 1/V_max results, if only as a first-order (albeit, somewhat biased) estimate of the XLF.

To estimate the random uncertainties on ${{\rm{\Phi }}}_{{\rm{L}}}(\mathrm{log}{L}_{{\rm{X}}})$, ${{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})$, and $\xi (\mathrm{log}{\lambda }_{{\rm{E}}})$, we follow the approach by Weigel et al. (2016; and see also Zhu et al. 2009; Gilbank et al. 2010). These errors are essentially Poisson errors on the number of sources in each bin, with effective weights applied based on the corresponding 1/V_max estimates. The exact prescriptions used to calculate these uncertainties are provided in Appendix A.

Having determined the XLF and having determined an initial guess for the BHMF and the ERDF via the 1/V_max method, we assume specific functional forms for all three distributions and fit the corresponding (differential) space densities. For the XLF, we assume a double power law of the following form:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{L}}}(\mathrm{log}{L}_{{\rm{X}}})=\displaystyle \frac{{dN}}{d\mathrm{log}L}={{\rm{\Phi }}}_{{\rm{L}}}^{* }\times {\left[{\left(\displaystyle \frac{{L}_{{\rm{X}}}}{{L}_{{\rm{X}}}^{* }}\right)}^{{\gamma }_{1}}+{\left(\displaystyle \frac{{L}_{{\rm{X}}}}{{L}_{{\rm{X}}}^{* }}\right)}^{{\gamma }_{2}}\right]}^{-1}.\end{eqnarray} \tag{ 9 }$$

For the fitting procedure, we parameterize the second power-law slope as γ₂ = γ₁ + _γ , with > 0 to avoid degeneracy between the two exponents. Given what is known about the rather universal nature of AGN LF (e.g., Shen et al. 2020, and references therein), in practice this means that our γ₁ and γ₂ correspond to what are often referred to as the faint- and bright-end LF slopes, respectively.

For the BHMF, ${{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})\equiv \tfrac{{dN}}{d\mathrm{log}{M}_{\mathrm{BH}}}$, we use a modified Schechter function, defined as:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})\propto \mathrm{ln}(10){{\rm{\Phi }}}^{* }{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\mathrm{BH}}^{* }}\right)}^{\alpha +1}\exp \left[-{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\mathrm{BH}}^{* }}\right)}^{\beta }\right].\end{eqnarray} \tag{ 10 }$$

This choice is motivated by the general shape of the galaxy stellar mass function (e.g., Baldry et al. 2012; Weigel et al. 2016; Davidzon et al. 2017, and references therein) and the close relations between SMBH and galaxy mass (Kormendy & Ho 2013).

For the ERDF, we use a broken power law, which we define as:

$$\begin{eqnarray}&&\xi (\mathrm{log}{\lambda }_{{\rm{E}}})=\displaystyle \frac{{dN}}{d\mathrm{log}{\lambda }_{{\rm{E}}}}\propto {\xi }^{* }\times {\left[{\left(\displaystyle \frac{{\lambda }_{{\rm{E}}}}{{\lambda }_{{\rm{E}}}^{* }}\right)}^{{\delta }_{1}}+{\left(\displaystyle \frac{{\lambda }_{{\rm{E}}}}{{\lambda }_{{\rm{E}}}^{* }}\right)}^{{\delta }_{2}}\right]}^{-1}.\end{eqnarray} \tag{ 11 }$$

This functional form is motivated by the fact that, qualitatively, the convolution of the BHMF with the ERDF should reproduce the LF, which directly links the bright-end slope of the XLF and the high-λ_E slope of the ERDF (e.g., Caplar et al. 2015; Weigel et al. 2017). Similar to what was done with the XLF, we parameterize δ₂ as δ₁ + _λ , with _λ > 0, thus linking δ₁ and δ₂ with the low- and high-λ_E ERDF slopes, respectively.

In the initial part of our analysis, we fit the 1/V_max measurements of the BHMF and the ERDF with these functional forms simply to obtain initial guesses for the more elaborate recovery of the intrinsic BHMF and ERDF. These same functional forms, however, are also relevant in various other parts of our analysis and interpretation.

To find the best-fitting parameters for the 1/V_max-based XLF, BHMF, and ERDF, we use the Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013). In Appendix B, we outline how we find the best-fitting functional form for the XLF. The BHMF and the ERDF are fit accordingly. We fit all three distributions independently. Note that this MCMC analysis only fits the functions given in Equations (9), (10), and (11) to 1/V_max estimates, and not directly to the observed astrophysical quantities (luminosities, masses, and/or Eddington ratios). In the next section, we describe a more elaborate method to correct for the limitations of the 1/V_max approach (discussed in Section 3.2), where we fit the aforementioned functions to the data directly.

3.4.1. The Basic Principle

To correct for sample truncation—i.e., the bias against low-mass and low–Eddington ratio AGNs due to the flux-limited nature of the Swift/BAT and BASS samples—we follow the approach of SW10 and S15. The technique requires the assumption of (intrinsic) functional forms for the BHMF and the ERDF, and thus it should be considered as a parametric maximum likelihood approach.

Our aim is to constrain the intrinsic bivariate distribution function ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},z)$, which also provides the BHMF and the ERDF, when integrated over $\mathrm{log}{\lambda }_{{\rm{E}}}$ and $\mathrm{log}{M}_{\mathrm{BH}}$, respectively:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}},z)={\displaystyle \int }_{\mathrm{log}{\lambda }_{{\rm{E}},\min ,{\rm{s}}}}^{\mathrm{log}{\lambda }_{{\rm{E}},\max ,{\rm{s}}}}{\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},z)d\mathrm{log}{\lambda }_{{\rm{E}}},\\ \xi (\mathrm{log}{\lambda }_{{\rm{E}}},z)={\displaystyle \int }_{\mathrm{log}{M}_{\mathrm{BH},\min ,{\rm{s}}}}^{\mathrm{log}{M}_{\mathrm{BH},\max ,{\rm{s}}}}{\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},z)d\mathrm{log}{M}_{\mathrm{BH}}.\end{array}\end{eqnarray} \tag{ 12 }$$

The integration limits $\mathrm{log}{\lambda }_{{\rm{E}},\min ,{\rm{s}}}$, $\mathrm{log}{\lambda }_{{\rm{E}},\max ,{\rm{s}}}$, $\mathrm{log}{M}_{\mathrm{BH},\min ,{\rm{s}}}$, and $\mathrm{log}{M}_{\mathrm{BH},\max ,{\rm{s}}}$ correspond to the minimum and maximum Eddington ratios and BH mass values that we are considering in our analysis (see Table 1). The bivariate distribution ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$ has physical units of space density, namely h³ Mpc⁻³ dex⁻² (i.e., per dex in M_BH and per dex in λ_E). As the redshift range is very small, we consider only a single redshift bin and we do not model the redshift evolution of the XLF, the BHMF, and/or the ERDF; this is why our bivariate distribution function ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$ does not depend directly on z. As noted above, we further assume that the ERDF is mass-independent. Under all these assumptions, the calculation of the BHMF and the ERDF from ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$ (the expressions in Equation (12)) reduces to:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})=\tilde{{\rm{\Phi }}}(\mathrm{log}{M}_{\mathrm{BH}})\times {\displaystyle \int }_{\mathrm{log}{\lambda }_{{\rm{E}},\ \min ,\ {\rm{s}}}}^{\mathrm{log}{\lambda }_{{\rm{E}},\max ,{\rm{s}}}}\tilde{\xi }(\mathrm{log}{\lambda }_{{\rm{E}}})d\mathrm{log}{\lambda }_{{\rm{E}}}\\ \xi (\mathrm{log}{\lambda }_{{\rm{E}}})=\tilde{\xi }(\mathrm{log}{\lambda }_{{\rm{E}}})\times {\displaystyle \int }_{\mathrm{log}{M}_{\mathrm{BH},\min ,{\rm{s}}}}^{\mathrm{log}{M}_{\mathrm{BH},\max ,{\rm{s}}}}\tilde{{\rm{\Phi }}}(\mathrm{log}{M}_{\mathrm{BH}})d\mathrm{log}{M}_{\mathrm{BH}}.\end{array}\end{eqnarray} \tag{ 13 }$$

The assumed functional forms for $\tilde{{\rm{\Phi }}}(\mathrm{log}{M}_{\mathrm{BH}})$ and $\tilde{\xi }(\mathrm{log}{\lambda }_{{\rm{E}}})$ are given by the right sides of Equations (10) and (11), respectively. These two functions are proportional to the actual BHMF and ERDF. To determine ${{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})$ and $\xi (\mathrm{log}{\lambda }_{{\rm{E}}})$, we have to define the M_BH and λ_E range that is being considered and marginalize over the other distribution. As Equation (13) shows, $\tilde{{\rm{\Phi }}}(\mathrm{log}{M}_{\mathrm{BH}})$ and $\tilde{\xi }(\mathrm{log}{\lambda }_{{\rm{E}}})$ retain their assumed shapes, but their normalizations are adjusted. In Section 3.4.4, we explain how these normalizations are determined after constraining the functional forms.

To estimate ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$, we compute the log-likelihood of observing our main BASS sample, which is given by

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=\sum _{i}^{{N}_{\mathrm{obs}}}\mathrm{ln}{p}_{i}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{N}_{{\rm{H}},i},\mathrm{log}{\lambda }_{{\rm{E}},i},{z}_{i}),\end{eqnarray} \tag{ 14 }$$

where p_i is the probability of observing object i with BH mass $\mathrm{log}{M}_{\mathrm{BH},i}$, Eddington ratio $\mathrm{log}{\lambda }_{{\rm{E}},i}$, and redshift z_i . The log-likelihood thus represents the sum of such (log) probabilities, over the total number of observed sources, N_obs. For an assumed ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$, p_i is given by the expected number of objects with similar properties to object i, relative to the total number of sources that are predicted to be observed; p_i thus encodes all the relevant selection effects, such as the flux limit of the survey, as we explain in detail in Section 3.4.2 immediately below. To estimate ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$, and thus the intrinsic BHMF and ERDF, we maximize the likelihood ${ \mathcal L }$, i.e., the probability of observing our ensemble of N_obs sources.

3.4.2. The Probability of Observing a Given Source

For a given bivariate distribution function Ψ, we express the probability of observing a specific object (index i) in the following way:

$$\begin{eqnarray}&&\begin{array}{l}{p}_{i}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},i},{z}_{i})\\ =\,\displaystyle \frac{{N}_{i}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},i},{z}_{i})}{{N}_{\mathrm{tot}}}\\ =\,\displaystyle \frac{1}{{N}_{\mathrm{tot}}}{\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i})\\ \,{{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},i},{z}_{i})\\ \,\,p(\mathrm{log}{N}_{{\rm{H}},i})\,p({z}_{i})\,\displaystyle \frac{{{dV}}_{C}({z}_{i})}{{dz}}.\end{array}\end{eqnarray} \tag{ 15 }$$

Note that this expression by itself does not correct for measurement uncertainties. We explain the mechanism that does ultimately account for these uncertainties in Section 3.4.3 below. For simplicity, we assume that p(z_i ) and $p(\mathrm{log}{N}_{{\rm{H}},i})$ are independent of Ψ in this expression. The specifics of how these functions are treated for our survey and redshift range are discussed in more detail below.

Here, we describe each of the the terms used for the calculation of p_i :

• 1.  
  Ψ, the intrinsic bivariate distribution function of the BH mass and Eddington ratio: for a given set of $\mathrm{log}{M}_{\mathrm{BH},i}$ and $\mathrm{log}{\lambda }_{{\rm{E}},i}$ values, the bivariate distribution function returns the space density of the objects with those properties, i.e., the intrinsic number of such AGNs per unit of comoving volume, per dex of BH mass and Eddington ratio.
• 2.  
  dV_C /dz, the comoving volume element: by multiplying the space density of the objects with $\mathrm{log}{M}_{\mathrm{BH},i}$ and $\mathrm{log}{\lambda }_{{\rm{E}},i}$ with the comoving volume element at z_i , the space density is converted to an absolute number of sources. The comoving volume element for the entire sky is defined as (e.g., Hogg 1999):

  $$\begin{eqnarray}&&\displaystyle \frac{{{dV}}_{C}(z)}{{dz}}=4\pi \displaystyle \frac{c}{{H}_{0}}\displaystyle \frac{{D}_{C}{\left(z\right)}^{2}}{E(z)},\end{eqnarray} \tag{ 16 }$$

  where D_C and E(z) were defined in Section 3.2.
• 3.  
  N_tot, the normalization: to obtain a probability, the number of sources with properties similar to object i is normalized by the total number of sources that are expected to be part of the sample after the selection effects have been taken into account. N_tot is given by an integral that accounts for all the quantities that affect the observability of our AGNs within the survey:

  $$\begin{eqnarray}&&\begin{array}{l}{N}_{\mathrm{tot}}=\iiint \displaystyle \int {{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},\mathrm{log}{N}_{{\rm{H}}},z)\\ \,{\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})\,p(\mathrm{log}{N}_{{\rm{H}}})\,p(z)\\ \,\displaystyle \frac{{{dV}}_{C}(z)}{{dz}}d\mathrm{log}{M}_{\mathrm{BH}}\,d\mathrm{log}{N}_{{\rm{H}}}\,d\mathrm{log}{\lambda }_{{\rm{E}}}\,{dz}.\end{array}\end{eqnarray} \tag{ 17 }$$

  Here, $d\mathrm{log}{M}_{\mathrm{BH}}$, $d\mathrm{log}{\lambda }_{{\rm{E}}}$, and dz are computed over the corresponding minimum and maximum values considered in the sample (see Table 1).
• 4.  
  $p(\mathrm{log}{N}_{{\rm{H}}})$: a bias-corrected intrinsic absorption function needs to be assumed, from which the $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})$ values of the intrinsic AGN population are drawn. As explained above, while the intrinsic luminosities of AGNs are the product of their M_BH and λ_E, any selection or distribution function that depends on their observed fluxes would also have to depend on (the distribution of the) line-of-sight obscuration, which we generally associate with circumnuclear (torodial) dusty gas. Ricci et al. (2015) derived an intrinsic $\mathrm{log}{N}_{{\rm{H}}}$ distribution specifically for Swift/BAT-detected AGNs, considering X-ray reflection from a torus. This model assumed an opening angle of 60° and considered two luminosity bins—above and below $\mathrm{log}({L}_{14-195\,\mathrm{keV}}/\mathrm{erg}\,{{\rm{s}}}^{-1})=43.7$. This intrinsic $\mathrm{log}{N}_{{\rm{H}}}$ distribution is shown in Figure 5.To test the impacts of different models, we also calculate the results for an absorption function with a torus opening angle of 35° (also from Ricci et al. 2015), the results of which are reported in Section 4 and shown in Appendix F. For consistency, the attenuation curve used when testing the effect of this absorption function also assumes an opening angle of 35° (see Figure 4). While one may expect that our results would have changed significantly if the input absorption function and attenuation curve were substantially different, we find that this small change in the assumed opening angle did not alter our results substantially. This indicates that our overall conclusions are robust against reasonably motivated changes to the model-dependent components.We argue that since the intrinsic absorption function derived in Ricci et al. (2015; in particular, with 60°) is calculated specifically for the BAT sample, this is the appropriate function for our purposes. We note that this absorption function is defined over the range $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=20\mbox{--}25$ in discrete bins of 1 dex width. While calculating N_tot for all AGNs, we assign $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})$ to each object from the underlying distribution by drawing values from this absorption function, taking the luminosity dependence into account. In this work, we choose to incorporate the luminosity-dependent absorption function as provided by Ricci et al. (2015), and do not impose any direct λ_E dependence on the absorption function. While a λ_E dependence is supported by some studies (e.g., Ricci et al. 2017b, and references therein), imposing such a dependence a priori would limit our ability to reveal the differences in the BHMF and the ERDF of obscured and unobscured sources. The distributions used in our calculations are shown in the bottom panel of Figure 5. When we compute the BHMF and the ERDF for Type 1 AGNs, we assume the intrinsic absorption function to be identical to the observed one, which is justified given the negligible effects of absorption in such sources. For Type 2 AGNs, we subtract the distribution of Type 1 AGNs from the bias-corrected overall distribution in each luminosity bin. These distributions are shown using the red dashed and solid lines (for high- and low-luminosity bins, respectively) in the bottom panel of Figure 5. Note that as the Ricci et al. (2015) absorption function is luminosity-dependent, the $p(\mathrm{log}{N}_{{\rm{H}}})$ used in this work is in fact $p(\mathrm{log}{N}_{{\rm{H}}}| {L}_{{\rm{X}}}({M}_{\mathrm{BH}},{\lambda }_{{\rm{E}}}))$.
• 5.  
  p(z): this is the redshift dependence term in the general expression, which in our case is considered constant (i.e., set to 1 in both the numerator and denominator of p_i ).
• 6.  
  Ω_sel: in Equation (15), the term ${{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},{\rm{i}}},{z}_{i})$ corresponds to the selection function of the survey. In its simplest form, Ω_sel returns the value of the flux–area curve, ${{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{F}_{{\rm{X}}})$, for each source (see Section 3.2). To predict the X-ray flux $\mathrm{log}{F}_{{\rm{X}},i}$ for source i with BH mass $\mathrm{log}{M}_{\mathrm{BH},i}$, an Eddington ratio $\mathrm{log}{\lambda }_{{\rm{E}},i}$, a column density $\mathrm{log}{N}_{{\rm{H}},i}$, and a redshift z_i , we perform the following sequence of calculations.First, we compute the corresponding bolometric luminosity, L_bol, using Equation (5) (in Appendix E, we experiment with a variable, luminosity-dependent bolometric correction); from L_bol, we then calculate the intrinsic, ultra-hard X-ray luminosity over the 14–195 keV range, L_{14–195 kev}, using Equation (2); from that we deduce the luminosity and flux as measured by BAT, taking into account the column density (using the curve shown in Figure 4) and the luminosity distance to the source (using z_i ). The calculated $\mathrm{log}{F}_{{\rm{X}},i}$ is then used to obtain the selection function, based on the curve shown in Figure 3.

Zoom In Zoom Out Reset image size

3.4.3. The Observed Bivariate Distribution Function

As presented in Equation (15), p_i is a four-dimensional normalized probability distribution function, and in order to properly account for the various selection functions in play, it needs to be convolved with the uncertainties in each of the four underlying parameters. Specifically, the selection effect–corrected p_i is as follows:

$$\begin{eqnarray}&&\begin{array}{l}{p}_{{\rm{i}},\ \mathrm{conv}}(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},i},{z}_{i})=\displaystyle \frac{{N}_{{\rm{i}},\ \mathrm{conv}}}{{N}_{\mathrm{tot}}}\\ =\,\displaystyle \frac{1}{{N}_{\mathrm{tot}}}\iiint \int {\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},\mathrm{log}{N}_{{\rm{H}}},z)\\ \ \times \,{{\rm{\Omega }}}_{\mathrm{sel}}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},\mathrm{log}{N}_{{\rm{H}}},z)\\ \,\times \,p(\mathrm{log}{N}_{{\rm{H}}})\,p(z)\,\displaystyle \frac{{{dV}}_{C}(z)}{{dz}}\\ \,\times \,\omega (\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},i},{z}_{i}| \mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},\mathrm{log}{N}_{{\rm{H}}},z)\\ \,\,d\mathrm{log}{M}_{\mathrm{BH}}\,\,d\mathrm{log}{\lambda }_{{\rm{E}}}\,\,d\mathrm{log}{N}_{{\rm{H}}}\,\,{dz}.\end{array}\end{eqnarray} \tag{ 18 }$$

Here, the denominator N_tot is the normalization constant of the original probability distribution p_i , and retains the value presented in Equation (17), while ω is the four-dimensional Gaussian distribution function that reflects the uncertainties related to the ith object, for all four directly or indirectly measured quantities. That is:

$$\begin{eqnarray}&&\begin{array}{l}\omega \left(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{{\rm{H}},i},{z}_{i}| \mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}},\mathrm{log}{N}_{{\rm{H}}},z\right)\\ =\,\displaystyle \frac{1}{4{\pi }^{2}{\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}\,{\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}\,{\sigma }_{\mathrm{log}{N}_{{\rm{H}}}}\,{\sigma }_{z}}\\ \times \,\exp \left(-\displaystyle \frac{{\left(\mathrm{log}{M}_{\mathrm{BH},i}-\mathrm{log}{M}_{\mathrm{BH}}\right)}^{2}}{2{\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}^{2}}-\displaystyle \frac{{\left(\mathrm{log}{\lambda }_{{\rm{E}},i}-\mathrm{log}{\lambda }_{{\rm{E}}}\right)}^{2}}{2{\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}^{2}}\right)\\ \times \,\exp \left(-\displaystyle \frac{{\left(\mathrm{log}{N}_{{\rm{H}},i}-\mathrm{log}{N}_{{\rm{H}}}\right)}^{2}}{2{\sigma }_{\mathrm{log}{N}_{{\rm{H}}}}^{2}}-\displaystyle \frac{{\left({z}_{i}-z\right)}^{2}}{2{\sigma }_{z}^{2}}\right),\end{array}\end{eqnarray} \tag{ 19 }$$

where ${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}$, ${\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}$, ${\sigma }_{\mathrm{log}{N}_{{\rm{H}}}}$, and σ_z correspond to the assumed uncertainties on $\mathrm{log}{M}_{\mathrm{BH}}$, $\mathrm{log}{\lambda }_{{\rm{E}}}$, $\mathrm{log}{N}_{{\rm{H}}}$, and z, respectively (see Table 1). Note that—at least for $\mathrm{log}\,{M}_{\mathrm{BH}}$ and $\mathrm{log}\,{\lambda }_{{\rm{E}}}$—these are dominated by systematic uncertainties inherent to the BH mass estimation methods we rely on.

Since we have spectroscopic redshifts for all sources in our sample, and the obscuring column densities $\mathrm{log}{N}_{{\rm{H}}}$ are derived from extensive spectral decomposition of multimission X-ray data (combining Swift/BAT and ancillary spectra at E < 10 keV), the convolution over $\mathrm{log}\,{N}_{{\rm{H}}}$ and z (with ${\sigma }_{\mathrm{log}\,{N}_{{\rm{H}}}}=0.25$ and σ_z = 0.005) does not change our results, but makes the process significantly more computationally expensive. Thus, in practice, we choose to only convolve over two dimensions:

$$\begin{eqnarray}&&\begin{array}{l}\omega (\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i}| \mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})\\ =\,\displaystyle \frac{1}{2\pi {\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}{\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}}\\ \times \exp \left(-\displaystyle \frac{{\left(\mathrm{log}{M}_{\mathrm{BH},i}-\mathrm{log}{M}_{\mathrm{BH}}\right)}^{2}}{2{\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}^{2}}-\displaystyle \frac{{\left(\mathrm{log}{\lambda }_{{\rm{E}},i}-\mathrm{log}{\lambda }_{{\rm{E}}}\right)}^{2}}{2{\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}^{2}}\right).\end{array}\end{eqnarray} \tag{ 20 }$$

We note that convolving over z and/or $\mathrm{log}{N}_{{\rm{H}}}$ may be, in principle, essential and beneficial for the analysis of samples that have photometric redshifts and/or more uncertain column density measurements.

In our main analysis, we have assumed ${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}={\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}\,=0.3$ dex (see Table 1). The assumption of no uncertainties (i.e., ${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}={\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}=0$) is used only to test and demonstrate our analysis framework in Appendix D. A higher value of ${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}={\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}=0.5$ dex is required for Type 2 AGNs, to take into account the rotation and aperture effects in the optical spectroscopy that is used to measure σ_⋆ (see, e.g., Gültekin et al. 2009 and Shankar et al. 2019). Additionally, we report results for a scenario where ${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}=0.3$ dex, along with a total uncertainty of 0.2 dex in bolometric luminosity (${\sigma }_{\mathrm{log}L,\mathrm{scatt}}$), which reflects both measurement and systematic uncertainties (dominated by the latter). As $\mathrm{log}{\lambda }_{{\rm{E}}}$ is calculated using two observed quantities (luminosity and mass), the uncertainty in the luminosity measurement also contributes to the total uncertainty in $\mathrm{log}{\lambda }_{{\rm{E}}}$. Therefore, the uncertainty in $\mathrm{log}{\lambda }_{{\rm{E}}}$ is ${\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}\simeq 0.36$(i.e., the log uncertainties in mass and luminosity added in quadrature). We note that this approach to calculating ${\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}$ is conservative, since the actual measurement errors on (X-ray) luminosities are much smaller than 0.2 dex, and since the intertwined nature of L, M_BH, and λ_E means that any systematic errors on $\mathrm{log}{M}_{\mathrm{BH}}$ and $\mathrm{log}{\lambda }_{{\rm{E}}}$ would be anticorrelated, rather than independent. As shown in Section 4, the various (conservative) levels of uncertainty we assume do not significantly affect our key results, attesting to the robustness of our conclusions.

In Appendix E, we report the results of an even more complex scenario, assuming a luminosity-dependent bolometric correction (from Duras et al. 2020), and an even larger error on $\mathrm{log}{\lambda }_{{\rm{E}}}$. Our framework has also been tested by using mock catalogs for this scenario, and the results are shown in Appendix D. In the main part of the analysis, we present the results of the simpler scenario with constant bolometric correction.

We stress again that we have not imposed any dependence of obscuration on luminosity, M_BH, or λ_E in the main analysis. Thus, any difference between the XLF, BHMF, and/or ERDF derived for the obscured and unobscured subsamples (Type 1 and 2 AGNs) would occur independent of our model assumptions. In Appendix E, where we present results for variable bolometric correction, some more complex assumptions are introduced. For example, Type 1 and Type 2 AGNs have different bolometric corrections. This might artificially introduce differences in the results between the two populations, therefore we keep such assumptions to a minimum in our main analysis.

3.4.4. Constraining Ψ

To determine the intrinsic bivariate distribution function, ${\rm{\Psi }}(\mathrm{log}{M}_{\mathrm{BH}},\mathrm{log}{\lambda }_{{\rm{E}}})$, we maximize the likelihood of observing our main input sample. Expressing the log-likelihood (Equation (14)) using Equation (15) gives:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}{ \mathcal L }=\displaystyle \sum _{i}^{{N}_{\mathrm{obs}}}\mathrm{ln}\left[{N}_{{\rm{i}},\ \mathrm{conv}}\right.\\ \left.\,\times \,(\mathrm{log}{M}_{\mathrm{BH},i},\mathrm{log}{\lambda }_{{\rm{E}},i},\mathrm{log}{N}_{\ {\rm{H}},i},{z}_{i})\right]-\mathrm{ln}{N}_{\mathrm{tot}}.\end{array}\end{eqnarray} \tag{ 21 }$$

Similar to the fitting procedure for the 1/V_max values, we use the MCMC python package emcee to maximize $\mathrm{ln}{ \mathcal L }$ (Foreman-Mackey et al. 2013). The number of free parameters in these fits, which is six (γ₁, _λ , ${\lambda }_{{\rm{E}}}^{* }$, α, β, and ${M}_{\mathrm{BH}}^{* }$), reflects the functional forms assumed for the BHMF and the ERDF. The initial guesses are based on the MCMC fits to the 1/V_max values, and 50 walkers are allowed to take 3000 steps (the chains usually converge within 2500 steps).

Equation (18) shows that the normalization of the intrinsic bivariate distribution function, Ψ*, does not affect the probability of observing source i, and thus also has no impact on the log-likelihood $\mathrm{ln}{ \mathcal L }$. When applying the MCMC procedure, we thus use constant normalizations for both the BHMF (${{\rm{\Phi }}}_{\mathrm{init}}^{* }$) and the ERDF (${\xi }_{\mathrm{init}}^{* }$). Having determined the best-fitting parameters, we renormalize Ψ and determine Ψ^* as follows:

$$\begin{eqnarray}&&{{\rm{\Psi }}}^{\ast }={a}_{{\rm{r}}{\rm{e}}{\rm{s}}{\rm{c}}{\rm{a}}{\rm{l}}{\rm{e}}}\times \displaystyle \frac{{N}_{{\rm{o}}{\rm{b}}{\rm{s}}}}{{N}_{{\rm{t}}{\rm{o}}{\rm{t}},\text{best-fit}}({{\rm{\Phi }}}_{{\rm{i}}{\rm{n}}{\rm{i}}{\rm{t}}}^{\ast },{\xi }_{{\rm{i}}{\rm{n}}{\rm{i}}{\rm{t}}}^{\ast })}.\end{eqnarray} \tag{ 22 }$$

Here, N_obs corresponds to the total number of objects observed in our sample, while N_{tot, best−fit} is determined by using Equation (17), the best-fitting parameters for Ψ, and the initial normalizations ${{\rm{\Phi }}}_{\mathrm{init}}^{* }$ and ${\xi }_{\mathrm{init}}^{* }$. The parameter a_rescale corresponds to an additional rescaling factor that can be used to correct for only partially available data (i.e., z and/or M_BH measurements). As our sample is spectroscopically complete (excluding the nonbeamed, non-Galactic sources, along with the others discussed in Section 2.1.2), we assume a_rescale = 1 for this work.

Finally, we use Equation (13) to determine the bias-corrected intrinsic BHMF and ERDF:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Phi }}}_{{\rm{M}}}(\mathrm{log}{M}_{\mathrm{BH}})={{\rm{\Psi }}}^{* }\tilde{{\rm{\Phi }}}(\mathrm{log}{M}_{\mathrm{BH}},{{\rm{\Phi }}}_{\mathrm{init}}^{* })\\ {\int }_{\mathrm{log}{\lambda }_{{\rm{E}},\min ,{\rm{s}}}}^{\mathrm{log}{\lambda }_{{\rm{E}},\max ,{\rm{s}}}}\tilde{\xi }(\mathrm{log}{\lambda }_{{\rm{E}}},{\xi }_{\mathrm{init}}^{* })d\mathrm{log}{\lambda }_{{\rm{E}}}\\ \xi (\mathrm{log}{\lambda }_{{\rm{E}}})={{\rm{\Psi }}}^{* }\tilde{\xi }(\mathrm{log}{\lambda }_{{\rm{E}}},{\xi }_{\mathrm{init}}^{* })\\ {\int }_{\mathrm{log}{M}_{\mathrm{BH},\min ,{\rm{s}}}}^{\mathrm{log}{M}_{\mathrm{BH},\max ,{\rm{s}}}}\tilde{{\rm{\Phi }}}(\mathrm{log}{M}_{\mathrm{BH}},{{\rm{\Phi }}}_{\mathrm{init}}^{* })d\mathrm{log}{M}_{\mathrm{BH}}.\end{array}\end{eqnarray} \tag{ 23 }$$

As Ψ^* has units of space density, namely h³ Mpc⁻³ dex⁻², so do Φ* and ξ*.

We have tested our methodology end to end to verify that we can indeed uncover the intrinsic distributions of BHMF and ERDF using this approach. To this end, we created two mock catalogs for both Type 1 and Type 2 AGNs, and tested our method using three different levels of uncertainty on the (mock) $\mathrm{log}{M}_{\mathrm{BH}}$ and $\mathrm{log}{\lambda }_{{\rm{E}}}$ values. These tests and mock catalogs are described in detail in Appendix D.

To determine the AGN XLF of the ultra-hard X-ray-selected AGNs, we use the 1/V_max approach. As discussed in Section 3.2, we bin our sources according to their intrinsic 14−195 keV luminosities. When determining the corresponding space densities, we take the flux–area curve and attenuation by high column densities into account (Figures 3, 4). Once we have determined Φ_L, we fit the XLF with a double power law (Equation (9)). The 1/V_max space densities are given in Appendix F (specifically, Table F5) and the best-fitting double-power-law parameters of the various subsamples we consider are given in Table 2. The analysis steps are further discussed in what follows.

Table 2. Best-fitting XLF Parameters ^a

Selection	$\mathrm{log}({L}_{{\rm{X}}}^{* }/\mathrm{erg}\,{{\rm{s}}}^{-1})$	$\mathrm{log}({{\rm{\Phi }}}_{{\rm{L}}}^{* }/{{h}}^{3}\ {\mathrm{Mpc}}^{-3})$	γ1	γ
BASS AGNs at z ≤ 0.3 (672; no obscuration correction)	${44.17}_{-0.18}^{+0.13}$	$-{4.71}_{-0.19}^{+0.25}$	${0.75}_{+0.12}^{-0.11}$	${1.64}_{-0.10}^{+0.11}$
0.01 ≤ z ≤ 0.3 (619; no obscuration correction)	${44.23}_{-0.24}^{+0.18}$	$-{4.79}_{-0.30}^{+0.34}$	${0.87}_{-0.21}^{+0.19}$	${1.59}_{-0.13}^{+0.14}$
0.01 ≤ z ≤ 0.3 (619; obscuration-corrected)	${44.24}_{-0.31}^{+0.24}$	$-{4.78}_{-0.39}^{+0.35}$	${0.92}_{-0.25}^{+0.21}$	${1.63}_{-0.19}^{+0.10}$
0.01 ≤ z ≤ 0.3, $\mathrm{log}\,{{\rm{N}}}_{{\rm{H}}}\lt 22$ (316)	${44.17}_{-0.20}^{+0.19}$	$-{5.04}_{-0.31}^{+0.36}$	${0.77}_{-0.26}^{+0.16}$	${1.53}_{-0.13}^{+0.21}$
0.01 ≤ z ≤ 0.3, $22\leqslant \mathrm{log}\,{{\rm{N}}}_{{\rm{H}}}\lt 24$ (263)	${44.30}_{-0.19}^{+0.06}$	$-{5.35}_{-0.04}^{+0.40}$	${0.99}_{-0.16}^{+0.13}$	${1.81}_{-0.18}^{+0.17}$
0.01 ≤ z ≤ 0.3, $\mathrm{log}\,{{\rm{N}}}_{{\rm{H}}}\geqslant 24$ (40)	${43.59}_{-0.38}^{+0.23}$	$-{4.51}_{-0.44}^{+0.50}$	${0.33}_{-0.40}^{+0.44}$	${1.81}_{-0.57}^{+0.13}$
Type 1 AGNs only (366; all Table 1 cuts)	${44.12}_{-0.22}^{+0.21}$	$-{4.75}_{-0.36}^{+0.27}$	${0.69}_{-0.26}^{+0.22}$	${1.60}_{-0.14}^{+0.22}$
Type 2 AGNs only (220; all Table 1 cuts)	${44.14}_{-0.33}^{+0.13}$	$-{4.88}_{-0.49}^{+0.29}$	${0.93}_{-0.26}^{+0.21}$	${1.75}_{-0.21}^{+0.20}$
All Table 1 Selection AGNs (586)	${44.10}_{-0.19}^{+0.20}$	$-{4.46}_{-0.44}^{+0.19}$	${0.85}_{-0.25}^{+0.17}$	${1.60}_{-0.12}^{+0.19}$

Note.

^a For the XLF, we assume a double-power-law shape (see Equation (9)). All samples described here exclude sources falling within the Galactic plane (−5 <latitude <5), weak X-ray associations, dual sources falling below the flux limit of the survey, beamed sources, and sources with z > 0.3.

Download table as:  ASCIITypeset image

In Figure 6, we present the XLF determined for the BASS/DR2 AGNs, regardless of their classification (i.e., both Type 1 and Type 2 sources). First, we use the redshifts, intrinsic 14–195 keV X-ray luminosities (L_{14–195 kev}), and column densities (N_H) of all 672 nonbeamed, high–Galactic latitude BASS/DR2 AGNs at z ≤ 0.3; their XLF is shown in the left panel with the blue open squares and dotted blue line. Second, the left panel of Figure 6 shows the XLF of the 619 sources that meet our redshift restrictions (0.01 ≤ z ≤ 0.3), without any corrections due to obscuration being applied to their selection function. The only noticeable difference between the XLF of all BASS/DR2 AGNs and the redshift-restricted subsample is in the low-L end. The redshift-restricted subsample lacks the 15 lowest-L AGNs (i.e., in the $\mathrm{log}[{L}_{{\rm{X}}}/\mathrm{erg}\,{{\rm{s}}}^{-1}]\lt 42.25$ bins), which are also some of the lowest-redshift AGNs in BASS/DR2 (i.e., z < 0.01). At higher luminosities, the number of AGNs in each luminosity bin for the second sample is nonzero, and usually equal to the number of AGNs in the supersample. As shown in Figure 6 and Table F5, in a few bins, the redshift-restricted sample has slightly higher space densities than the supersample, even when the latter has more sources in those bins. This is possible because the lower redshift cut leads to smaller V_max (and therefore higher 1/V_max) values in some cases.

Zoom In Zoom Out Reset image size

We further analyze the XLF of our samples in the right panel of Figure 6. In addition to the XLF of the redshift-restricted subsample, we show how this XLF changes when obscuration corrections are applied to the selection function. The obscuration corrections result in slightly higher values of 1/V_max at all luminosities (as the V_max value of the obscured object is smaller). The 1/V_max values and the corresponding errors are given in Tables F2 and F5. The best-fitting parameters for each XLF are shown in Table 2.

The correction due to obscuration is generally very small, which is to be expected, as the radiation in the 14−195 keV regime only starts to get attenuated beyond $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\gtrsim 23$ in the local universe (as shown in Figure 4). Therefore, taking the effect of obscuration into account for ultra-hard X-rays may not lead to significant changes, unless the sample at hand preferentially selects heavily obscured sources. For example, the lowest-luminosity bin of the redshift-restricted sample contains a single object with $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=23.15$, and therefore the effect of obscuration in that bin is noticeable, as shown in the right panel of Figure 6. As our sample includes both unobscured and heavily obscured sources, we apply an obscuration correction to all other 1/V_max calculations, except for the two cases shown in Figure 6. In the figure, we show the 1/V_max values and XLF fit to the sample used for the BHMF/ERDF analysis (after applying all cuts from Table 1).

Figure 6 also shows previous determinations of the ultra-hard XLF of low-redshift AGNs by Sazonov et al. (2007), Tueller et al. (2008), and Ajello et al. (2012). We also show the binned 1/V_max XLF measurements of the latter study (black symbols). To convert these previous results to the 14–195 keV range we use here, we have assumed Γ = 1.8. For the results by Ajello et al. (2012), this implies L_{14−195 keV} = 2.31 ×L_{15−55 keV}. To convert from L_{20−40 keV} and L_{17−60 keV} to L_{14−195 keV}, we used multiplicative factors of 4.34 and 2.33, respectively.

The left panel in Figure 6 shows that our parent sample of (unbeamed) BASS/DR2 AGNs reaches down to hard X-ray luminosities that are ≈1 dex fainter than the faintest luminosities covered by Ajello et al. (2012), which relied on the 60 month Swift/BAT all-sky catalog. As a natural result of the lower redshift cut at z = 0.01, our redshift-restricted sample (0.01 ≤ z ≤ 0.3) only goes down to $\mathrm{log}({L}_{{\rm{X}}}/\mathrm{erg}\,{{\rm{s}}}^{-1})=42.48$, which is higher than the lowest luminosities reached by Ajello et al. (2012).

Overall, Figure 6 shows that the 1/V_max-based Φ_L values for our BASS/DR2 AGNs are in excellent agreement with previous results.

The size of the BASS sample allows us to determine the XLF for various subsets of AGNs. In Figure 7, we split the AGNs by $\mathrm{log}{N}_{{\rm{H}}}$, and calculate the binned XLF separately for unabsorbed Compton-thin $[\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\lt 22$], absorbed Compton-thin $[22\leqslant \mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\lt 24]$, and Compton-thick $[\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\geqslant 24]$ sources. For comparison, the blue solid lines in the top panels of Figure 7 illustrate the XLF for our redshift-restricted unbeamed AGN sample. For reference, in the top right panel of Figure 7, we also show the XLF of Compton-thick AGNs reported by Akylas et al. (2016). They used 53 Compton-thick AGNs selected from the the same parent catalog as the one we use (the Swift/BAT 70 month catalog), and determined a relatively low "break" luminosity in the 14−195 keV range, $\mathrm{log}({L}_{{\rm{X}}}^{* }/\mathrm{erg}\,{{\rm{s}}}^{-1})\simeq 42.8$, compared to our value of ${43.59}_{-0.23}^{+0.38}$. Note that we also have 53 Compton-thick objects in our parent sample, 11 of which are at z < 0.01, one of which falls above z > 0.3, and another one of which is a faint dual source. Our individual 1/V_max Φ_L values are generally consistent with the fit by Akylas et al. (2016), within errors, in all but the highest luminosity bin, although at $\mathrm{log}({L}_{14\mbox{--}195\,\mathrm{kev}}/\mathrm{erg}\,{{\rm{s}}}^{-1})\gt 44$ the Akylas et al. (2016) Compton-thick XLF lies above our data points. This may be caused by the different volumes considered, and/or by the fact that Akylas et al. (2016) uses a Poissonian MLE that considers individual sources, rather than a fit to 1/V_max (see also Loredo 2004).

Zoom In Zoom Out Reset image size

To further demonstrate the fractional densities of the $\mathrm{log}{N}_{{\rm{H}}}$ subsamples, the bottom panels of Figure 7 show the ratios between the space densities of each subsample compared to the entire sample, as a function of L_{14–195 keV}. These ratios are computed using the 1/V_max Φ_L values, and the errors on the ratios are calculated using the Wilson Score Interval method (Wilson 1927), which provides binomial confidence intervals without any assumption about the symmetry of the error bars.

We stress again that the XLFs and the related ratios shown in Figure 7 are estimated using the 1/V_max method, and thus are for the observed sample of sources only. While these XLFs account for the effects of absorption on the sources observed in our survey, they do not account for absorbed populations that are completely missing from the sample due to selection biases (that are accounted for in more sophisticated studies, e.g., Ueda et al. 2014 and Ananna et al. 2019). Deriving such intrinsic N_H-dependent XLFs and ratios requires a much more involved approach than the 1/V_max estimates we use here, which is beyond the scope of the present work.

We discuss the insights from Figure 7 in more detail in Section 5.

To determine the BHMF and the ERDF for local AGNs, we use the M_BH and λ_E measurements of 586 ultra-hard X-ray-selected BASS/DR2 AGNs, as described in Section 2.2. Since the M_BH (and thus λ_E) of Type 1 (broad-line) and Type 2 (narrow-line) BASS/DR2 AGNs are determined through two different approaches, with different systematics and potentially different selection effects in play, we first treat these two subsets separately, and only then address the BHMF and the ERDF of the total BASS/DR2 AGN sample.

To gain initial guesses for the two distribution functions, we use the 1/V_max approach, assume functional forms, and fit the individual Φ and ξ values independently. For the BHMF and the ERDF, we assume a modified Schechter function and a double power law, respectively (see Equations (10) and (11) in Section 3.3). To correct for sample truncation, i.e., the bias against low-mass and low–Eddington ratio AGNs, we use the parametric maximum likelihood approach outlined in Section 3.4. The 1/V_max values for both the BHMF and the ERDF are given in Appendix F.

Figure 8 shows the BHMFs (top panels) and the ERDFs (center panels) for Type 1 and Type 2 BASS AGNs (the left and right columns), respectively. The two bottom panels of Figure 8 show the XLFs reproduced from these BHMFs and ERDFs, as explained below. Figure 9 presents the BHMF, ERDF, and XLF for the overall sample (Type 1 + Type 2) in a similar format. We recall that the normalizations are kept constant during the bias correction, and the method of computing them is described in Section 3.4.1. We are thus left with six free parameters: the break and two slopes of the BHMF (${M}_{\mathrm{BH}}^{* }$, α, and β), and the break and two slopes of the ERDF (${\lambda }_{{\rm{E}}}^{* }$, δ₁, and δ₂). We note that the high-λ_E slope, δ₂, is parameterized as δ₁ + _λ , with _λ > 0. We vary δ₁ and _λ in the MCMC. We also recall that for the bias correction we assume uncertainties of 0.3 dex and 0.5 dex on both $\mathrm{log}{M}_{\mathrm{BH}}$ and $\mathrm{log}{\lambda }_{{\rm{E}}}$ (see Equations (18) and (20)), as well as an additional scenario with ${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}=0.3$ and ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.3$ (i.e., ${\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}=0.36$). The results from all these different σ values reassuringly converge on the same solution, as shown in the figure. In Appendix E, we present the BHMF/ERDF calculated by assuming a luminosity-dependent bolometric correction from Duras et al. (2020). Note that we prefer the constant bolometric correction because it minimizes the number of assumptions we have to make, and because there is some conflict between different prescriptions of luminosity-dependent bolometric corrections (shown in Figure 6 of Duras et al. 2020), and it is unclear which prescription is the most accurate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The contour plots presenting the likelihoods resulting from our MCMC analysis are shown in Figure F1 in Appendix F. The best-fitting BHMF and ERDF parameters for all three samples are given in Tables 3 and 4, respectively. Note that we also report the results assuming an attenuation curve and absorption function calculated using a torus opening angle of 35° (as discussed in Section 3.4.2) in these tables. As shown in Figure F2, our final results for both torus geometries are consistent with each other.

Table 3. Sample Truncation-corrected BHMFs ^a and Fits to /V_max Values for the BHMFs for Type 1 AGNs, Type 2 AGNs, and Both Samples Together ^b

	$\mathrm{log}({M}_{\mathrm{BH}}^{* }/{M}_{\odot })$	$\mathrm{log}({{\rm{\Psi }}}^{* }/{{h}}^{3}\ {\mathrm{Mpc}}^{-3})$	α	β
All
Intrinsic (σ = 0.3)	${7.88}_{-0.22}^{+0.21}$	−3.52	$-{1.576}_{-0.078}^{+0.147}$	${0.593}_{-0.069}^{+0.078}$
Intrinsic (σ = 0.3; ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.2$)	${7.92}_{-0.22}^{+0.13}$	−3.67	$-{1.530}_{-0.094}^{+0.114}$	${0.612}_{-0.063}^{+0.053}$
Intrinsic (σ = 0.5)	${7.67}_{-0.20}^{+0.25}$	−3.37	$-{1.26}_{-0.11}^{+0.19}$	${0.630}_{-0.086}^{+0.065}$
Intrinsic (σ = 0.3; OA = 35°)	${7.92}_{-0.27}^{+0.18}$	−3.49	$-{1.576}_{-0.157}^{+0.057}$	${0.600}_{-0.082}^{+0.066}$
1/Vmax	${8.12}_{-0.11}^{+0.14}$	$-{4.33}_{-0.37}^{+0.22}$	$-{1.06}_{-0.30}^{+0.12}$	${0.574}_{-0.040}^{+0.112}$
Type 1
Intrinsic (σ = 0.3)	${7.97}_{-0.30}^{+0.17}$	−4.19	$-{1.753}_{-0.088}^{+0.137}$	${0.561}_{-0.073}^{+0.062}$
Intrinsic (σ = 0.3; ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.2$)	${7.93}_{-0.22}^{+0.25}$	−4.27	-1.73 ± 0.11	${0.566}_{-0.071}^{+0.085}$
Intrinsic (σ = 0.5)	${7.91}_{-0.27}^{+0.16}$	−4.27	-1.56 ± 0.13	${0.590}_{-0.050}^{+0.104}$
1/Vmax	${8.73}_{-0.31}^{+0.26}$	$-{5.10}_{-0.50}^{+0.26}$	$-{1.35}_{-0.23}^{+0.18}$	${0.681}_{-0.114}^{+0.087}$
Type 2
Intrinsic(σ = 0.3)	${7.82}_{-0.26}^{+0.15}$	−3.6	$-{1.16}_{-0.20}^{+0.14}$	${0.637}_{-0.075}^{+0.067}$
Intrinsic (σ = 0.3; ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.2$)	${7.79}_{-0.22}^{+0.17}$	−3.64	$-{1.18}_{-0.17}^{+0.16}$	${0.617}_{-0.044}^{+0.100}$
Intrinsic (σ = 0.5)	${7.76}_{-0.15}^{+0.19}$	−3.6	$-{0.99}_{-0.19}^{+0.21}$	${0.703}_{-0.069}^{+0.086}$
Intrinsic (σ = 0.3; OA = 35°)	${7.73}_{-0.12}^{+0.25}$	−3.44	$-{1.26}_{-0.16}^{+0.17}$	${0.635}_{-0.053}^{+0.085}$
1/Vmax	${8.102}_{-0.095}^{+0.172}$	$-{4.33}_{-0.23}^{+0.19}$	$-{1.04}_{-0.29}^{+0.30}$	${0.732}_{-0.050}^{+0.074}$

Notes.

^a We assume a modified Schechter function (see Equation (10)) for the BHMF. ^b We use a constant bolometric correction (Equation (2)) to compute these results.

Download table as:  ASCIITypeset image

Table 4. Sample Truncation-corrected ERDFs ^a and Fits to /V_max Values for the ERDFs for Type 1 AGNs, Type 2 AGNs, and the Full AGN Sample ^b

	$\mathrm{log}{\lambda }_{{\rm{E}}}^{* }$	$\mathrm{log}({\xi }^{* }/{{h}}^{3}\,{\mathrm{Mpc}}^{-3})$	δ1	λ
All
Intrinsic (σ = 0.3)	−1.338 ± 0.065	−3.64	${0.38}_{-0.10}^{+0.12}$	${2.260}_{-0.082}^{+0.121}$
Intrinsic (σ = 0.3; ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.2$)	$-{1.286}_{-0.076}^{+0.059}$	−3.76	0.40 ± 0.12	${2.322}_{-0.105}^{+0.095}$
Intrinsic (σ = 0.3; OA = 35°)	$-{1.332}_{-0.068}^{+0.077}$	−3.68	${0.484}_{-0.133}^{+0.091}$	${2.210}_{-0.106}^{+0.066}$
Intrinsic (σ = 0.5)	$-{1.249}_{-0.061}^{+0.059}$	−3.8	${0.28}_{-0.12}^{+0.13}$	${2.72}_{-0.14}^{+0.13}$
1/Vmax	$-{1.19}_{-0.21}^{+0.23}$	$-{3.76}_{-0.30}^{+0.18}$	$-{0.02}_{-0.39}^{+0.29}$	${2.06}_{-0.27}^{+0.30}$
Type 1
Intrinsic (σ = 0.3)	−${1.152}_{-0.089}^{+0.053}$	−4.08	0.30 ± 0.14	${2.51}_{-0.15}^{+0.11}$
Intrinsic (σ = 0.3; ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.2$)	$-{1.138}_{-0.058}^{+0.070}$	−4.09	${0.27}_{-0.12}^{+0.13}$	${2.57}_{-0.16}^{+0.12}$
Intrinsic (σ = 0.5)	$-{1.103}_{-0.054}^{+0.065}$	−4.23	${0.13}_{-0.14}^{+0.15}$	${2.97}_{-0.19}^{+0.18}$
1/Vmax	$-{1.06}_{-0.25}^{+0.28}$	$-{4.02}_{-0.32}^{+0.22}$	$-{0.51}_{-0.41}^{+0.53}$	${2.57}_{-0.45}^{+0.33}$
Type 2
Intrinsic (σ = 0.3)	−${1.657}_{-0.087}^{+0.064}$	−3.82	${0.376}_{-0.253}^{+0.099}$	2.50 ± 0.17
Intrinsic (σ = 0.3; ${\sigma }_{\mathrm{log}L,\mathrm{scatt}}=0.2$)	$-{1.628}_{-0.084}^{+0.079}$	−3.84	${0.32}_{-0.20}^{+0.18}$	2.50 ± 0.17
Intrinsic (σ = 0.3; OA = 35°)	$-{1.675}_{-0.067}^{+0.091}$	−3.8	${0.33}_{-0.20}^{+0.15}$	${2.51}_{-0.16}^{+0.14}$
Intrinsic (σ = 0.5)	$-{1.593}_{-0.096}^{+0.080}$	−3.92	${0.30}_{-0.20}^{+0.21}$	${2.53}_{-0.11}^{+0.26}$
1/Vmax	$-{1.87}_{-0.40}^{+0.38}$	$-{3.74}_{-0.43}^{+0.35}$	$-{0.50}_{-0.56}^{+1.08}$	${2.30}_{-0.69}^{+0.81}$

Notes.

^a We assume a double-power-law shape for the ERDF (see Equation (11)). ^b We use a constant bolometric correction (Equation (2)) to compute these results.

Download table as:  ASCIITypeset image

As discussed in Section 1, and shown in detail in Appendix D, the bolometric LF corresponds to the convolution of the BHMF and the ERDF. By "reversing" the bolometric correction, the XLF can thus be predicted from the best-fit BHMF and ERDF. We test if the bias-corrected BHMFs and ERDFs for each subset of AGNs (i.e., Type 1, Type 2, and the overall sample) allow us to predict the corresponding observed XLF. Note that for the convolution, we use the normalized ERDF (see Equation (32)). The normalization of the predicted XLF is thus driven by the normalization of the BHMF.

Figures 8 and 9 illustrate the importance of bias correction. For the BHMF, relying solely on the model fit to 1/V_max values would have led us to underestimate the space density of low-mass AGNs at the lowest-mass bin ($\mathrm{log}[{M}_{\mathrm{BH}}/{M}_{\odot }]=6.65$) by ⪆ 1 dex for both Type 1 and Type 2 sources. At low M_BH and/or λ_E, the 1/V_max method underestimates the intrinsic space density of AGNs, due to the survey's incompleteness. At high M_BH and/or λ_E, the 1/V_max method overestimates the intrinsic AGN space densities, due to the uncertainties associated with the key AGN parameters. As shown by the mock catalogs (Figures D1 and D2), and Figure 17 of S15, as measurement uncertainty (${\sigma }_{\mathrm{log}{M}_{\mathrm{BH}}}$, ${\sigma }_{\mathrm{log}{\lambda }_{{\rm{E}}}}$) increases, this overestimation increases at the high-M_BH and/or high-λ_E end. This effect is a manifestation of the so-called Eddington bias (Eddington 1913): the uncertainty causes objects from lower-M_BH (or lower-λ_E) bins to scatter into higher-M_BH (or higher-λ_E) bins, and vice versa. Intrinsically, there are always fewer objects with higher M_BH (or λ_E), therefore the scattering from the lower to the higher bins causes significant overestimation of the space densities in the higher bins. The bottom panels in Figures 8 and 9 demonstrate that our convolution-based reconstruction of the XLF matches what we measure directly from observations.

Figure 8 shows that the normalization of the bias-corrected BHMF of Type 2 AGNs is higher than that of Type 1 AGNs at most masses. The bias-corrected ERDF of Type 2 AGNs has higher space densities at $\mathrm{log}{\lambda }_{{\rm{E}}}\leqslant -1.7$. Beyond this point, the ERDF of Type 2 AGNs drops off rapidly below the ERDF of Type 1 AGNs. This is because the break in the ERDF of Type 2 AGNs is at log λ_E = −${1.657}_{-0.064}^{+0.087}$, which is significantly below the break of the Type 1 AGNs at log λ_E = −${1.152}_{-0.053}^{+0.089}$. Figure 9 shows the Type 1 and Type 2 BHMFs, ERDFs, and XLFs, along with the overall sample results.

Figure 10 illustrates how we can use the bivariate distribution (i.e., Ψ) to reproduce the (univariate) BHMF and ERDF, following Equation (23). In the top panels, we show how the bivariate distribution varies as a function of $\mathrm{log}{M}_{\mathrm{BH}}$ and $\mathrm{log}{\lambda }_{{\rm{E}}}$ for Type 1 (left panels) and Type 2 (right panels) AGNs. We also indicate the intervals over which we integrate to produce the BHMFs and ERDFs, as shown in the lower panels. In the middle panels, we show the reconstructed BHMFs for two bins of log λ_E (bin widths of 0.3 dex), as well as the integrated BHMF (over all λ_E we consider). Similarly, the bottom panels show the reconstructed ERDFs for two log M_BH bins (bin widths of 0.3 dex), as well as the integrated ERDF over all M_BH. We note that, in a graphical sense, the reconstructed XLFs would correspond to integrating the bivariate distribution along antidiagonal stripes (see the top left panel in Figure 10), illustrating that the bivariate distribution function fully captures the statistical properties of AGNs.

Zoom In Zoom Out Reset image size

Note that there is some degeneracy between several pairs of parameters describing the fitting functions, as shown in our MCMC chain contour plots in Appendix F (Figure F1). For all three AGN populations (Type 1, Type 2, and overall), the parameters of the BHMF show significant correlation. The pairs of parameters (α, β) and $(\alpha ,\mathrm{log}{M}_{\mathrm{BH}}^{* })$ seem to be anticorrelated, while the pair $(\beta ,\mathrm{log}{M}_{\mathrm{BH}}^{* })$ is positively correlated. These trends are expected: while fitting the same intrinsic population, if we fix the break in the BHMF at a higher mass, the slope at low mass (α) has to be shallower, and the slope at high mass (β) has to be steeper, to compensate for the higher mass break and to produce a good fit. For the ERDF, the slopes (δ₁–_λ ) are negatively correlated for all three samples. The break of the ERDF ($\mathrm{log}{\lambda }_{* }$) is weakly positively correlated with δ₁, and shows no significant correlation with _λ . The figure also shows that even when these functions converge to the same distribution for one parameter, it may occupy distinctly different locations in six-dimensional parameter space. These degeneracies should be kept in mind when one tries to directly compare individual fitting parameters within our own analysis (e.g., between Type 1 and Type 2 AGNs), or when comparing our best-fit parameters to those found in other studies. In what follows, instead of comparing the values of individual parameters, we often refer to the similarities (or lack thereof) in the shapes of certain fitting functions shown in the figure.

In Figure 11, we compare the BASS BHMFs and ERDFs for Type 1 AGNs (left panels) and Type 2 AGNs (right panels) to previously published determinations of these distributions. In the top panels, we show the 1/V_max Φ values and the bias-corrected BHMF in blue (left panels; Type 1) and red (right panels; Type 2). In the top left panels, the green data points show the Φ values determined by Greene & Ho (2007, 2009), based on a large sample of SDSS AGNs with broad Hα lines, which is not corrected for sample truncation. The yellow data points and dashed–dotted line illustrate the 1/V_max Φ values and the incompleteness-corrected distribution, respectively, as determined by SW10, based on the quasar sample of HES (Wisotzki et al. 2000). Note that SW10 does not correct for measurement uncertainties. The purple dashed lines in Figure 11 show the BHMF and ERDF determined by Schulze et al. (2015), based on a joint analysis of 1 < z < 2 AGNs from the SDSS, zCOSMOS, and VVDS samples (Schneider et al. 2010; Lilly et al. 2009 and Gavignaud et al. 2008, respectively). Note that the Schulze et al. (2015) BHMF and ERDF for Type 2 AGNs are predicted from the corresponding distributions for Type 1 AGNs, by assuming a luminosity-dependent fraction of obscured systems, which is in turn determined from X-ray surveys. The higher space densities of the Schulze et al. (2015) BHMFs and ERDFs are expected, given the well-known redshift evolution of AGN abundance (e.g., Aird et al. 2015; Caplar et al. 2015; Aird et al. 2018; Shen et al. 2020). The turquoise shaded region in the top panels of Figure 11 also shows the total local BHMF—that is, including both active and inactive SMBHs—determined by Shankar et al. (2009). The total BHMF naturally occupies a much higher space density than our BASS AGNs–only BHMF. The recent study by Shankar et al. (2020) used a different set of M_BH − σ relations (discussed in Section 2.2) to derive an updated total BHMF. However, as our analysis relies on the more widely used M_BH − σ relation of Kormendy & Ho (2013), we do not compare our results with the Shankar et al. (2020) BHMF.

Zoom In Zoom Out Reset image size

We used our BASS/DR2 sample to directly verify that our assumption of a mass-independent ERDF is a reasonable assumption, at least for the data in hand. To this end, we divided each of the Type 1 and Type 2 AGN subsamples into two broad mass bins, $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\leqslant 7.8$ and $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\geqslant 8.2$. The 0.4 dex wide gap in $\mathrm{log}{M}_{\mathrm{BH}}$ was imposed to minimize the "mixing" between the mass bins, due to uncertainties on the M_BH estimation. We then derived the 1/V_max measurements (which are susceptible to selection biases, as discussed in detail in Section 3), the associated functional fits of these measurements, and the bias-corrected intrinsic functional forms for each of the four subsamples (Type 1s and Type 2s, low-mass and high-mass), following the same methodology as used for our main analysis.

The results of this analysis are shown in Figure 12. We find that the shapes of the bias-corrected intrinsic distributions of the low-mass and high-mass bins—that is, the power-law exponents and the locations of the breaks—are in excellent agreement (for each of the AGN subtypes), as can also be seen from the best-fit parameters (Table F1). The low-mass subsets naturally have higher number densities (i.e., normalizations), as expected, given the generally decreasing nature of the BHMF with increasing M_BH. This analysis indicates that the shape of the ERDF is indeed mass-independent, for both obscured and unobscured AGNs.

Zoom In Zoom Out Reset image size

Our analysis of the BASS/DR2 sample offers several insights concerning the role of obscuration in the distributions describing the (low-redshift) AGN population. First, we have fully considered the effect of line-of-sight obscuration on the derived intrinsic luminosity of every AGN in our sample, and thus on the derived XLF (as well as the BHMF and the ERDF). This is motivated by the significant attenuation expected for highly obscured sources $[\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\gt 23$], even in the ultra-hard 14−195 keV band (i.e., Figure 4).

Second, the BASS sample allowed us to construct and explore the XLF for AGNs of several N_H regimes. The top panels of Figure 7 show the XLFs in each of the three N_H bins, while the bottom panels show the fractions of objects in each N_H bin relative to all AGNs (in a luminosity-resolved way; see the details in Section 4.2). We find that (1) the fraction of unabsorbed AGNs increases with luminosity; (2) the fraction of Compton-thin sources decreases with luminosity; and (3) the fraction of Compton-thick objects remains roughly constant with luminosity.

The observation that unobscured sources dominate the high-L end of the AGN LF is often interpreted as evidence for the so-called "receding torus model" (e.g., Lawrence 1991; Simpson 2005). In this model, the covering factor of the dusty torus, which obscures the nuclear region, decreases with increasing AGN luminosity—a trend that is observed in many (X-ray) AGN samples and surveys (e.g., Lawrence & Elvis 1982; Steffen et al. 2003; Barger et al. 2005; Simpson 2005; Hasinger et al. 2005; La Franca et al. 2005; Treister & Urry 2006; Maiolino et al. 2007; Brusa et al. 2010; Burlon et al. 2011). As the opening angle increases, the probability of detecting an obscured source decreases. Some studies have suggested that unobscured AGNs with higher bolometric luminosities may have less (dusty) torus material, based on the mid-IR to bolometric ratio, again supporting the receding torus model (e.g., Treister et al. 2008). Some alternative explanations have also been suggested. For example, Akylas & Georgantopoulos (2008) showed that the photoionization of the obscuring screen around AGNs roughly reproduces the relationship between the fraction of obscured AGNs and X-ray luminosity, as observed by XMM-Newton and Chandra. Hönig & Beckert (2007) suggested that the Eddington limit on a clumpy torus may also cause the obscured fraction to decrease with luminosity. It is important to note that the apparent decrease in the fraction of obscured AGNs at high luminosities could be partially due to selection effects, as suggested by, e.g., Treister & Urry (2006)—further emphasizing the need for large, highly complete samples, drawn from homogeneous input catalogs (such as BASS).

Many of these studies have often made the (pragmatic) assumption that Compton-thin AGNs trace all obscured objects (e.g., Ueda et al. 2014; see also Hickox & Alexander 2018). Specifically, several analyses of AGN LFs have assumed that the abundance of AGNs stays constant throughout the $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=24-26$ regime, thus predicting significant space densities of $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\gt 25$ AGNs (e.g., Ueda et al. 2014; Aird et al. 2015; Buchner et al. 2015; Ananna et al. 2019). In contrast, only a handful of objects (three) with $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\geqslant 25$ are observed in the (70 month) all-sky Swift/BAT survey, conducted in the ultra-hard 14–195 keV band (Ricci et al. 2015, 2017b). To see how our BASS-based results relate to this issue, we consider the Ananna et al. (2019) XLF, which suggests that Compton-thick AGNs would comprise ≈50% of all AGNs in the local universe. If we apply the appropriate BAT flux limits to that XLF (i.e., the 70 month survey), and limit the XLF to $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\lt 25$, we find that <10% of the observable sample would in fact be Compton-thick, which is consistent with the results shown in Figure 7. We also note that the Ricci et al. (2015) study, from which the intrinsic N_H distribution shown in Figure 5 was derived, also presents the observed N_H distribution of the Swift/BAT AGN sample used in that study (their Figure 4, bottom panel). The observed Compton-thick fraction in that study is ≲10%, again consistent with the result we show in Figure 7 here (note that the Ricci et al. (2015) analysis also includes AGNs at z < 0.01).

As the attenuation curve in Figure 4 shows, the observed luminosity of an object with $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})\geqslant 25$ is ≤ 5% of its intrinsic luminosity in the 14–195 keV band, therefore we have to probe very faint fluxes to be able to detect such heavily obscured AGNs. Vito et al. (2018), Yan et al. (2019), and Carroll et al. (2021) find evidence for IR-selected luminous AGNs that completely lack X-ray detection, even in the NuSTAR 3–79 keV band. ²⁹ It is therefore possible that even high-energy X-ray observations, such as the Swift/BAT survey, do not individually identify the most obscured Compton-thick objects, and thus the fractions reported here for the high-N_H subsample should be treated with some caution.

Considering the physical driver for the trends linking (Compton-thin) obscuration and accretion power, our results in Figure 7 seem at face value to be consistent with the receding torus scenario. However, note that there is no clear drop in the (relative) space densities of Compton-thick AGNs with increasingly high luminosities (i.e., there is no downward trend in the bottom right panel of Figure 7). Indeed, the small sample size and the correspondingly large errors mean we cannot robustly rule out the possibility that the luminosity dependence of the Compton-thick space densities is consistent with that of Compton-thin sources. With this caveat in mind, if the Compton-thick fraction indeed remains (roughly) constant with luminosity, it may lead to some important insights regarding the distribution of circumnuclear matter in AGNs.

Fabian et al. (2006, 2008) and Fabian et al. (2009) suggested that the radiation pressure exerted on the dusty torus gas is crucial for understanding the links between AGN accretion power and obscuration. This was corroborated by the BASS/DR1-based study by Ricci et al. (2017b), which tied radiation pressure, gravity, and orientation angle together, and suggested an explanation for the relation between the fraction of obscured sources and luminosity. In this scenario, the fraction of obscured (Compton-thin) sources fundamentally depends on the Eddington ratio, which dictates the effective radiation pressure on the dusty torus gas. As the luminosity exceeds the effective Eddington limit for the dusty gas, the torus material is pushed away from the central engines, thus significantly decreasing the abundance of high-λ_E, high-N_H sources. The observed luminosity dependence is then simply a consequence of the λ_E dependence, dictated by the (limited) range of M_BH probed in AGN surveys. Most importantly, Ricci et al. (2017b) showed that the fraction of Compton-thick AGNs is independent of luminosity, indicating that these clouds are apparently unaffected by radiation pressure, or that the effective λ_E threshold for Compton-thick clouds is too high and is seldom exceeded. ³⁰

Our results provide further support for this radiation pressure–driven scenario. The ERDFs we constructed for Type 1 and Type 2 AGNs (Figure 8) clearly show that obscured AGNs (Type 2 sources) are indeed increasingly rare beyond $\mathrm{log}{\lambda }_{{\rm{E}}}^{* }\simeq -1.7$, which is remarkably consistent with the effective Eddington limit for dusty gas $[\mathrm{log}{\lambda }_{{\rm{E}}}\approx -1.7$ for $\mathrm{log}({N}_{{\rm{H}}}/{\mathrm{cm}}^{-2})=22;$ see Ricci et al. (2017b) and references therein]. The downturn in the space densities of unobscured (Type 1) AGNs occur at higher accretion rates, $\mathrm{log}{\lambda }_{{\rm{E}}}^{* }\,=$ −${1.152}_{-0.053}^{+0.089}$, and may instead be linked to the physics of accretion disks and/or the (circumnuclear-scale) fueling mechanisms. We therefore propose that the much higher space densities of unobscured AGNs (relative to obscured AGNs) at high λ_E could be naturally explained through the effect of radiation pressure on the dusty obscuring material.

Indeed, looking more closely at the different shapes of the ERDFs for Type 1 and Type 2 AGNs, another line of interpretation suggests itself, relating small-scale physics to large-scale population statistics. First, high accretion rates (i.e., high λ_E) can be triggered by major galaxy mergers, and theoretical models (e.g., Hopkins et al. 2006a, 2006b; Blecha et al. 2018) and observations (Treister et al. 2012; Glikman et al. 2012; Banerji et al. 2015; Glikman et al. 2015; Ricci et al. 2017c; Glikman et al. 2018; Koss et al. 2018; Banerji et al. 2021) suggest that luminous, merger-triggered AGNs start from a highly obscured state (i.e., Type 2), but eventually blow away the obscuring material to become unobscured AGNs (i.e., Type 1). The ratio of Type 2 to Type 1 AGN number densities at high λ_E could therefore reflect the ratio of the short duration of the obscured phase to a much longer, unobscured phase. In contrast, at low λ_E, where accretion is less violent or disruptive, the ratio of Type 2 to Type 1 AGNs likely reflects the geometry of circumnuclear obscuration, as in the traditional unification scheme that is well established locally (Barthel 1989; Antonucci 1993; Urry & Padovani 1995). At intermediate λ_E, a transition occurs, where a mix of AGN types coexist, and where the break in the ERDF of Type 2 AGNs likely reflects the onset of significant radiative feedback, as explained above. If the picture of obscured AGNs transitioning into unobscured AGNs during a merger-driven accretion episode is correct, then we would expect that, at higher redshifts, more high-λ_E obscured AGNs would be transitioning into unobscured AGNs, since both merger rates and gas fractions generally increase with redshift. Therefore, we predict that the break in the Type 2 AGN ERDF at higher redshift should occur at higher λ_E (e.g., Jun et al. 2021). Quantitative exploration of these ideas will be deferred to a future work.

We caution that our sample still has relatively few high-luminosity obscured sources. Still, the very existence of such sources indicates that the receding torus scenario is at the very least incomplete (see also Bär et al. 2019), and our extensive bias corrections suggest that the dearth of high-λ_E obscured AGNs is real, and cannot be easily explained by obvious observational biases.

In Figure 8, the similarity in shape of the BHMFs for Type 1 and Type 2 AGNs, and the difference in the shape of the ERDFs, could potentially imply that the observed difference between the two populations is mainly due to the different distributions in the Eddington ratio. It could also highlight the possibility of a fundamental difference between the two AGN populations. As explained in the preceding section, this difference may be related to AGN fueling mechanisms, including galaxy mergers and large-scale environments (e.g., as seen in the BASS-based clustering analysis of Powell et al. 2018); to smaller-scale feedback mechanisms, including the radiation-regulated unification scheme (as supported by the BASS/DR1-based analysis of Ricci et al. 2017b); or perhaps to some other mechanisms. In any case, any comparison of our results to other studies should take into account these differences between Type 1 and Type 2 AGNs.

In Figure 11, we compare the BHMFs and ERDFs for BASS Type 1 and Type 2 AGNs with other BHMFs and ERDFs reported in the literature. We compare our BHMF for the Type 1 AGN subsample to the 1/V_max-based BHMF of the local broad-line SDSS AGNs derived by Greene & Ho (2007) and Greene & Ho (2009). We find that the BASS observed space densities are higher, by ≥ 0.5 dex, in the range $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\gt 7.0$, which may have been caused by differences in the selection method. We note that the Greene & Ho (2007) and Greene & Ho (2009) studies focused on a broad (optical) selection function, which included host-dominated continuum sources with broad Hα lines, and not just quasar-like AGNs. Our higher space densities thus indicate that the ultra-hard X-ray selection of AGNs provides larger, likely more complete samples even of (high-luminosity) broad-line AGNs, compared to SDSS. Some of this discrepancy may be related to the bright flux limit imposed as part of the SDSS spectroscopy (i > 15 mag), which BASS does not impose. Specifically, the fact that our 1/V_max space densities extend to higher masses implies that the X-ray selection is probing a different M_BH (and possibly λ_E) range.

SW10 assumed a modified Schechter function–shaped BHMF and a Schechter function–shaped ERDF, evaluated over essentially the same redshift interval as the present work (z < 0.3). While SW10 does correct for incompleteness at the low-M_BH/λ_E range, unlike the 1/V_max method, it does not account for measurement uncertainty. For the BHMF, SW10 report the following best-fitting parameters: $\mathrm{log}({M}_{\mathrm{BH}}^{* }/{M}_{\odot })=8.11$, $\mathrm{log}({{\rm{\Phi }}}^{* }/{{h}}^{3}\,{\mathrm{Mpc}}^{-3})=-5.10$, α = − 2.11, and β = 0.5. The 1/V_max values of the BHMF of SW10 are ≃ 0.5 dex lower at $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\leqslant 9$ than this work. These differences in 1/V_max between SW10 and this work are also reflected in the discrepancy between the intrinsic BHMFs, particularly at intermediate masses. As both SW10 and the present work focus on the same redshift regime, the differences between the two observed populations likely arise from the distinct selections by the optical and ultra-hard X-ray bands.

The differences in the two samples become clearer when we look at the ERDF. There is a noticeable discrepancy between SW10 and our Type 1 results—in both the observed data points and the intrinsic functions. Again, the discrepancy in the observed samples might be due to the different selection methods. Optical surveys are more likely to be biased toward high-luminosity AGNs, because in order for an object to be identified as an AGN, it has to dominate over host galaxy emission. Therefore, this sample would be skewed toward high luminosity (and λ_E). Additionally, our sample of Type 1 AGNs covers all objects with broad Hα lines, including intermediate types (such as 1.5 and 1.9). These intermediate types tend to be somewhat obscured, even in the Compton-thin regime, which means that if the radiation-regulated unification model is indeed the dominant mechanism, a sample of luminous quasars (such as that of SW10) would be skewed toward higher λ_E, compared to the more complete BASS sample. As the SW10 study does not account for measurement uncertainty, that could also lead to part of the discrepancy.

S15 presented the BHMF and ERDF for Type 1 AGNs as being in the redshift range 1 < z < 2. Additionally, S15 predicted the BHMF and ERDF of Type 2 AGNs as being in the same redshift range, by using a luminosity-dependent obscured fraction function (from Merloni et al. 2014). By comparing our local BASS results with the 1 < z < 2 S15 results, it appears that, at higher redshifts, there are more high-mass AGNs of both types, whereas in the local universe the lower-mass AGNs become more abundant. This suggests that many high-mass SMBHs have become inactive between 1 ≲ z ≲ 2 and z ≲ 0.3, or that the average accretion rate has decreased over time (e.g., due to fewer mergers or because interstellar gas has been depleted). The S15 Type 1 AGN ERDF agrees well with our Type 1 ERDF, while both the normalization and the shape of the S15 Type 2 ERDF are significantly different from our Type 2 results. Specifically, the shapes of the two S15 ERDFs are consistent with each other, whereas the shapes of our Type 1 and Type 2 ERDFs are significantly different. Note that the predicted S15 Type 2 estimates are based on several assumptions. As suggested in Merloni et al. (2014), even though the obscuration fraction is reported as redshift-independent, some redshift dependence is still seen at high luminosities. Additionally, Merloni et al. (2014) reported some issues that could lead to incorrect estimations of intrinsic luminosities (i.e., due to incorrect assumptions about the complex geometry of the obscurer). As this fraction is constrained using ≤ 10 keV data, these uncertainties could be significant, due to the degeneracy of AGN spectral parameters (Gilli et al. 2007; Ricci et al. 2017a; Ananna et al. 2020a). If the z ≃ 1 − 2 ERDF represents the underlying Type 2 population at that redshift, it would imply that AGN activity has decreased over time. The overall normalizations of the BHMFs at 1 < z < 2 are also higher than the BHMFs of the local universe, which supports the decreased activity scenario.

A highly useful observational constraint on theoretical models of SMBH evolution is the AGN duty cycle—that is, the fraction of all SMBHs (including inactive SMBHs) that are actively accreting, above a certain Eddington ratio (and at any given cosmic epoch). The AGN duty cycle has been addressed by numerous observational and theoretical studies. Some hydrodynamical galaxy simulations have tried to quantify the AGN duty cycle by tracing the gas inflow onto the central SMBHs. Novak et al. (2011) simulated a single galaxy and found that the SMBH accretes at $\mathrm{log}{\lambda }_{{\rm{E}}}\gt -3$ for ≲ 30% of the time span (12 Gyr) covered by the simulation. Angles-Alcazar et al. (2021) recently reported a duty cycle of ∼0.25 at z < 1.1. Phenomenological AGN population models can constrain and/or deduce the AGN duty cycle by linking the observed (redshift-resolved) AGN LF with the local (active and inactive) BHMF, or indeed the (integrated) BH mass density, generally following the Soltan (1982) argument (e.g., Cavaliere & Padovani 1989; Marconi et al. 2004; Shankar et al. 2009). For example, Shankar et al. (2009) found that, in the local universe, less than 1% of all SMBHs should be considered as active (i.e., accreting at the fiducial accretion rate of their model; see their Figure 7). The active fraction or duty cycle is often defined as a ratio of LF to mass function, independent of Eddington ratio. Both simulations and population synthesis studies typically have to assume an AGN radiative efficiency (Shankar et al. 2009), and often also have to assume an ERDF, or even a universal λ_E (e.g., Shankar et al. 2013; Weigel et al. 2017).

Highly complete AGN surveys naturally provide the observational benchmark for the AGN duty cycle. For example, Goulding et al. (2010) reported an active fraction of ≈0.27, based on a volume-limited, mid-IR-selected sample of D < 15 Mpc galaxies—although their definition of "active" includes AGNs with Eddington ratios as low as $\mathrm{log}{\lambda }_{{\rm{E}}}\geqslant -5$.

There are theoretical, phenomenological, and observational lines of argument for the AGN duty cycle to depend on galaxy and/or BH mass, and perhaps on other properties as well (e.g., galaxy environment and/or clustering; Haiman & Hui 2001; Martini & Weinberg 2001; Shen et al. 2007; White et al. 2008; Shen et al. 2009; Shankar et al. 2010b, 2010a). In Figure 13, we show the AGN duty cycle in the local universe, based on our BASS/DR2 AGN sample. In the left panel of Figure 13, we show the λ_E-dependent duty cycle, expressed as the cumulative probability of having $\mathrm{log}{\lambda }_{{\rm{E}}}$ greater than a given value, $P(\gt \mathrm{log}{\lambda }_{{\rm{E}}})$. We calculate this probability by integrating over all BH mass bins and $\mathrm{log}{\lambda }_{{\rm{E}}}$ bins above a given value, then dividing by the integrated total local BHMF (i.e., including inactive SMBHs), taken from Shankar et al. (2009). The Shankar et al. (2009) BHMF was compiled by taking into account the dispersions in all the local SMBHs that were available at the time (e.g., McLure & Dunlop 2002; Marconi & Hunt 2003; Tundo et al. 2007). As discussed in Section 2.2, Shankar et al. (2020) report a BHMF ∼4 lower than the Shankar et al. (2009) BHMF (in terms of space densities at all masses), as the later study uses a recalibrated M–σ relationship. As Koss et al. (2022c) use the canonical Kormendy & Ho (2013) prescription to estimate the masses for BASS/DR2 objects, we compare our results to the Shankar et al. (2009) BHMF in Figure 13.

Zoom In Zoom Out Reset image size

The M_BH-integrated AGN duty cycles for the entire BASS/DR2 sample relative to the Shankar et al. (2009) BHMFs at $\mathrm{log}{\lambda }_{{\rm{E}}}=-2,-1$, and 0 are about $P(\mathrm{log}{\lambda }_{{\rm{E}}}\gt -2)\,\simeq$ 2.85 ×10⁻², $P(\mathrm{log}{\lambda }_{{\rm{E}}}\gt -1)\,\simeq$ 6.45 × 10⁻⁴, and $P(\mathrm{log}{\lambda }_{{\rm{E}}}\gt 0)\,\simeq$ 1.61 × 10⁻⁶, respectively. At the lowest-λ_E threshold that is reasonable for radiatively efficient SMBH accretion, we obtain $P(\mathrm{log}{\lambda }_{{\rm{E}}}\gt -3)\simeq$ 0.1–0.16. According to the Shankar et al. (2009) SMBH mass function, the total mass density of all SMBHs in the local universe is ≈ (3 − 5) × 10⁻⁵ M_⊙ h³ Mpc⁻³, while the active SMBH mass density is ${3.58}_{-0.23}^{+0.41}\times {10}^{-4}$ M_⊙ h³ Mpc⁻³ (i.e., 6%–10% of the total SMBH mass density).

Considering the Type 1 and Type 2 AGNs in our samples, their duty cycles are essentially identical at the fiducial threshold corresponding to $P(\mathrm{log}\,{\lambda }_{{\rm{E}}}\simeq -2.5)$. For lower-threshold Eddington ratios ($\mathrm{log}{\lambda }_{{\rm{E}}}\lt -2.5$), the cumulative duty cycle is slightly higher for Type 2 sources than it is for Type 1 sources (but within the 1σ error budget), while for higher-threshold λ_E, the duty cycle of Type 1 AGNs is significantly higher. This means that a lower fraction of obscured AGNs have such high λ_E, which is not surprising given the differences between the Type 1 and Type 2 ERDFs (Figure 8).

The right panel of Figure 13 shows the fraction of active SMBHs (AGNs) with $\mathrm{log}{\lambda }_{{\rm{E}}}\gt -2$ among the total SMBH population (including inactive BHs), as a function of M_BH. The general trend for all BASS AGNs is that the AGN fraction decreases with increasing M_BH. This general trend is in agreement with what was found in several previous studies, including both direct observations (e.g., Greene & Ho 2007; Goulding et al. 2010) and population models (e.g., Marconi et al. 2004; Shankar et al. 2009, 2013). Our AGN fraction (10%–16%) is about an order of magnitude above what was found by Shankar et al. (2009) using their fiducial model (<1% active), which assumes a constant radiative efficiency of 0.065. Our results are in slightly better agreement with the redshift-dependent Eddington ratio model from the Shankar et al. (2009) study. This demonstrates the importance of independent, observational determinations of the AGN ERDFs, fractions, and duty cycles to (phenomenological) models of the cosmic evolution of SMBHs. We note that these trends may evolve with redshift.

Among the more nuanced trends in the right panel of Figure 13, we note that at lower masses, $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\lt 9$, the fraction of Type 2 AGNs with $\mathrm{log}{\lambda }_{{\rm{E}}}\gt -2$ is higher than that of Type 1 AGNs. This is driven by the fact that the space density of Type 2 AGNs is higher than that of Type 1 AGNs at $\mathrm{log}{\lambda }_{{\rm{E}}}\lt -1.7$, as shown in Figure 8. This trend flips completely if the threshold is moved to $\mathrm{log}{\lambda }_{{\rm{E}}}\geqslant -1$ (not shown in the figure), and Type 1 AGNs dominate the fraction of AGNs at all mass bins.

Assuming that the trends found here hold for other surveys (and out to higher redshifts), they highlight why different survey strategies may lead to ambiguous or contradictory conclusions about the obscured AGN fraction. Specifically, wide-field surveys that pick up the rarest, highest-M_BH systems are expected to be biased toward unobscured systems; conversely, deeper (and narrower) surveys that uncover the more abundant lower-M_BH population may be (mildly) biased toward obscured AGNs. Mateos et al. (2017) reported that by studying the torus structure and covering factor of X-ray-selected samples (at <10 keV), they find that a significant population of obscured objects should exist at high luminosities, and are missed by X-ray surveys. These heavily obscured high-luminosity AGNs have been identified in recent studies using IR, optical, and X-ray data (e.g., Yan et al. 2019; Carroll et al. 2021). Treister et al. (2010) also reached similar conclusions by analyzing IR-selected sources. As discussed in Section 5.1, and as implied by Figure 4 and these results, given Swift/BAT's current flux limits, many heavily obscured (i.e., $\mathrm{log}[{N}_{{\rm{H}}}/{\mathrm{cm}}^{-2}]\gt 25$), massive AGNs may still be undetected by BAT.

We finally use the BHMF we derive from the BASS/DR2 data to speculate about the BH to stellar mass ratio. The local universe provides ample evidence for a close relation between the mass of SMBHs and the stellar mass of their host galaxies (particularly their bulge components). The ratio of SMBHs to stellar mass lies in the range $-3.55\lt {\rm{l}}{\rm{o}}{\rm{g}}({M}_{{\rm{B}}{\rm{H}}}/{M}_{{\rm{g}}{\rm{a}}{\rm{l}}})\lt -2.31$ (e.g., Marconi & Hunt 2003; Häring & Rix 2004; Sani et al. 2011; Kormendy & Ho 2013; Marleau et al. 2013; McConnell & Ma 2013; Reines & Volonteri 2015, and references therein). If SMBH masses scale (roughly) linearly with host galaxy masses, then the corresponding mass functions should have similar shapes, with a horizontal shift that scales as M_BH/M_gal. ³¹ We explore the M_BH − M_gal scaling relationship by comparing the break in our BHMF (${M}_{\mathrm{BH}}^{* }$) with the break in the galaxy stellar mass function (${M}_{{\rm{g}}{\rm{a}}{\rm{l}}}^{\ast }$).

Some recent determinations of the galaxy stellar mass function, based on large optical low-redshift surveys, find breaks at $10.5\lesssim {\rm{l}}{\rm{o}}{\rm{g}}({M}_{{\rm{g}}{\rm{a}}{\rm{l}}}^{\ast }/{M}_{\odot })\lesssim 10.7$ (e.g., MacLeod et al. 2010; Baldry et al. 2012; Weigel et al. 2016), with relatively limited variance between different galaxy types (see the discussion in, e.g., Moffett et al. 2016; Davidzon et al. 2017, and references therein). Comparing these galaxy stellar mass function breaks directly with the BHMF break we find for all BASS/DR2 AGNs $[{\rm{l}}{\rm{o}}{\rm{g}}({M}_{{\rm{B}}{\rm{H}}}^{\ast }/{M}_{\odot })\,={7.88}_{-0.22}^{+0.21}$], we get −2.82 $\lesssim \,{\rm{l}}{\rm{o}}{\rm{g}}({M}_{{\rm{B}}{\rm{H}}}^{\ast }/{M}_{{\rm{g}}{\rm{a}}{\rm{l}}}^{\ast })\,\lesssim$ −2.62. This agrees well with the range of ${\rm{l}}{\rm{o}}{\rm{g}}({M}_{{\rm{B}}{\rm{H}}}^{\ast }/{M}_{{\rm{g}}{\rm{a}}{\rm{l}}}^{\ast })$ derived from direct measurements in individual systems.

We caution that these results are highly speculative, as they inherently link the active BHMF (i.e., the BHMF of AGNs) to a quantity of stellar mass functions that are dominated by inactive galaxies, thus assuming that the BASS-based ${M}_{\mathrm{BH}}^{* }$ is representative of all SMBHs and/or that BASS AGN hosts are indistinguishable from the general galaxy population. The actual analysis of BASS AGN hosts and the measurement of their stellar masses is beyond the scope of this work, but is expected to be pursued in a future BASS study.

Figure 14 shows the galaxy mass function (shifted left by 2.88 dex, for ease of comparison in the figure) and the AGN BHMF for our overall sample, as well as for the Type 1 and Type 2 subsamples. In this figure, we use a modified Schechter function to represent the AGN samples, and a double Schechter function to represent the galaxy stellar mass function (e.g., Weigel et al. 2016, 2017). The shape of the galaxy stellar mass function may differ from the BH mass functions because of the different functional forms used to fit the data. However, this might also mean that the ratio of BH to galaxy mass, or even bulge to galaxy mass, varies with galaxy mass (e.g., Bell et al. 2017). Fitting a modified Schechter function to galaxy stellar masses is beyond the scope of this work, but will also be explored in our future project.

Zoom In Zoom Out Reset image size

We stress again that these simplistic galaxy–BH scaling relationships provide only a limited view of the relations between (BASS) AGNs and their hosts, as some studies suggest that close SMBH–host links should only be applicable to the (true) bulge or spheroidal components of galaxies and/or to certain types of galaxies (see the detailed discussions in, e.g., Graham et al. 2011 and Kormendy & Ho 2013).