# BASS XXX: Distribution Functions of DR2 Eddington-ratios, Black Hole Masses,
and X-ray Luminosities

Tonima Tasnim Ananna Department of Physics and Astronomy, Dartmouth College,
6127 Wilder Laboratory, Hanover, NH 03755, USA  Anna K. Weigel Modulos AG,
Technoparkstrasse 1, CH-8005 Zurich, Switzerland  Benny Trakhtenbrot School of
Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel  Michael J.
Koss Eureka Scientific, 2452 Delmer Street Suite 100, Oakland, CA 94602-3017,
USA  Space Science Institute, 4750 Walnut Street, Suite 205, Boulder, Colorado
80301, USA  C. Megan Urry Department of Physics and Yale Center for Astronomy
& Astrophysics, Yale University, PO BOX 201820, New Haven, CT 06520, USA
Claudio Ricci Núcleo de Astronomía de la Facultad de Ingeniería, Universidad
Diego Portales, Av. Ejército Libertador 441, Santiago 22, Chile  Kavli
Institute for Astronomy and Astrophysics, Peking University, Beijing 100871,
People’s Republic of China  Ryan C. Hickox Department of Physics and
Astronomy, Dartmouth College, 6127 Wilder Laboratory, Hanover, NH 03755, USA
Ezequiel Treister Instituto de Astrofísica, Facultad de Física, Pontificia
Universidad Católica de Chile, Casilla 306, Santiago 22, Chile  Franz E. Bauer
Instituto de Astrofísica and Centro de Astroingeniería, Facultad de Física,
Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile
Millennium Institute of Astrophysics, Nuncio Monseñor Sótero Sanz 100, Of 104,
Providencia, Santiago, Chile  Space Science Institute, 4750 Walnut Street,
Suite 205, Boulder, Colorado 80301  Yoshihiro Ueda Department of Astronomy,
Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan
Richard Mushotzky Department of Astronomy, University of Maryland, College
Park, MD 20742, USA  Federica Ricci Dipartimento di Fisica e Astronomia,
Universitá di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy  INAF -
Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Gobetti,
93/3, I40129 Bologna, Italy  Kyuseok Oh Korea Astronomy & Space Science
institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku,
Kyoto 606-8502, Japan  JSPS Fellow  Julian E. Mejía-Restrepo European Southern
Observatory, Casilla 19001, Santiago 19, Chile  Jakob Den Brok Argelander-
Institut fur Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn,
Germany  Daniel Stern Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Drive, MS 169-224, Pasadena, CA 91109, USA
Meredith C. Powell Kavli Institute of Particle Astrophysics and Cosmology,
Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA  Turgay Caglar
Leiden Observatory, PO Box 9513, 2300 RA, Leiden, The Netherlands  Kohei
Ichikawa Frontier Research Institute for Interdisciplinary Sciences, Tohoku
University, Sendai 980-8578, Japan  O. Ivy Wong CSIRO Astronomy & Space
Science, PO Box 1130, Bentley, WA 6102, Australia  ICRAR-M468, University of
Western Australia, Crawley, WA 6009, Australia  Fiona A. Harrison  Cahill
Center for Astronomy and Astrophysics, California Institute of Technology,
Pasadena, CA 91125, USA  Kevin Schawinski Modulos AG, Technoparkstrasse 1,
CH-8005 Zurich, Switzerland

###### Abstract

We determine the low-redshift X-ray luminosity function (XLF), active black
hole mass function (BHMF), and Eddington-ratio distribution function (ERDF)
for both unobscured (Type 1) and obscured (Type 2) active galactic nuclei
(AGN) using the unprecedented spectroscopic completeness of the BAT AGN
Spectroscopic Survey (BASS) data release 2. In addition to a straightforward
1/$V_{\mathrm{max}}$ approach, we also compute the intrinsic distributions,
accounting for sample truncation by employing a forward modeling approach to
recover the observed BHMF and ERDF. As previous BHMFs and ERDFs have been
robustly determined only for samples of bright, broad-line (Type 1) AGNs
and/or quasars, ours is the first directly observationally constrained BHMF
and ERDF of Type 2 AGN. We find that after accounting for all observational
biases, the intrinsic ERDF of Type 2 AGN is significantly skewed towards lower
Eddington ratios than the intrinsic ERDF of Type 1 AGN. This result supports
the radiation-regulated unification scenario, in which radiation pressure
dictates the geometry of the dusty obscuring structure around an AGN.
Calculating the ERDFs in two separate mass bins, we verify that the derived
shape is consistent, validating the assumption that the ERDF (shape) is mass
independent. We report local AGN duty cycle as a function of mass and
Eddington ratio, by comparing the BASS active BHMF with local mass function
for all SMBH. We also present the $\log N{-}\log S$ of Swift-BAT 70-month
sources.

Active galactic nuclei (16), Supermassive black holes (1663), X-ray surveys
(1824), X-ray active galactic nuclei (2035), Luminosity function (942),
Accretion (14)

##  1 Introduction

Supermassive black holes (SMBHs) are found at the centers of nearly all
massive galaxies, and are understood to co-evolve with their host galaxies
(see Kormendy & Ho, 2013 for a review). Actively accreting SMBHs, identified
by their high luminosities or rates of accretion, are known as active galactic
nuclei (AGNs). The space density of AGN as a function of luminosity—i.e., the
AGN luminosity function (LF)— represents a key statistical measure for the AGN
population which allows us to constrain the abundance and growth history of
SMBHs (e.g., Soltan, 1982).

The redshift-resolved AGN LF and space density have had major impacts on our
understanding of the evolving SMBH population. For example, it is used to
determine the epoch of peak SMBH growth at around $z\sim 2$ (e.g., Barger et
al., 2001; Ueda et al., 2003; Hasinger et al., 2005; Croom et al., 2009; Ueda
et al., 2014; Ananna et al., 2020b), quite similar to the peak in cosmic star
formation activity (e.g., Lilly et al., 1996; Madau et al., 1998; Zheng et
al., 2009; Madau & Dickinson, 2014; Aird et al., 2015; Caplar et al., 2015).
It is also clear that the space densities of low-luminosity AGN peak at lower
redshifts compared to higher-luminosity systems (so-called “downsizing”, see,
e.g., Barger et al. 2001; Ueda et al. 2003; Miyaji et al. 2015; Brandt &
Alexander 2015; Ueda et al. 2014; Ananna et al. 2020b). At yet higher
redshifts, the AGN LF can help constrain the contribution of accreting SMBHs
to cosmic reionization (e.g., Willott et al., 2010a; Kashikawa et al., 2015;
Giallongo et al., 2015; Ricci et al., 2017d; Parsa et al., 2018; Matsuoka et
al., 2018; Ananna et al., 2020b). Indeed, when used as a key ingredient in
phenomenological population models, the evolving AGN LF is used to trace the
growth of SMBHs throughout cosmic history, ultimately accounting for the local
population of (relic) SMBHs, and even SMBH-host relations (e.g., Soltan, 1982;
Marconi et al., 2004; Shankar et al., 2009; Ueda et al., 2014; Aird et al.,
2015; Buchner et al., 2015; Caplar et al., 2018; Ananna et al., 2019). The AGN
LF is therefore a very useful statistical tool for understanding the AGN
population and its evolution (see, e.g., Brandt & Alexander, 2015 for a
review).

The AGN LF alone, however, cannot constrain the crucial characteristics of the
underlying SMBH population. This is because the AGN (bolometric) luminosity is
essentially the product of two more fundamental properties of a black hole—
its mass ($M_{\rm BH}$) and relative accretion rate, which we parameterize as
the dimensionless Eddington-ratio ($\lambda_{\rm E}\equiv L_{\rm bol}/L_{\rm
Edd}$), that is

$L_{\rm bol}\propto M_{\rm BH}\times(L_{\rm bol}/L_{\rm Edd})=M_{\rm
BH}\times\lambda_{\rm E}\,.$ (1)

Therefore, only after measuring the underlying BH mass and Eddington-ratio for
sizable, representative AGN samples can we decisively answer questions such as
when was the epoch during which the most massive BHs ($M_{\rm BH}\gtrsim
10^{9}\,M_{\odot}$) grew most of their mass. Several studies show that such
high-mass BHs accreted at maximal Eddington rates, reaching the Eddington
limit at $z\gtrsim 5$ (e.g., Willott et al., 2010b; Trakhtenbrot et al., 2011;
De Rosa et al., 2014). In the local Universe, on the other hand, it seems that
lower mass ($\sim 10^{6}-10^{8}\,M_{\odot}$) AGN with lower $\lambda_{\rm E}$
dominate space density distributions, even among the most luminous AGN (i.e.
quasars; e.g., McLure & Dunlop, 2004; Netzer & Trakhtenbrot, 2007; Schulze &
Wisotzki, 2010; Schulze et al., 2015). To obtain a complete census of the AGN
population, we thus have to consider three key distribution functions: the AGN
(bolometric) LF, the active black hole mass function (BHMF)111Note that in
this work we use BHMF to denote the mass function of the active SMBH
population alone, unless stated otherwise, whereas the total BHMF is the sum
of active and inactive BHMFs, and the Eddington-ratio distribution function
(ERDF), after correcting for obscuration, uncertainties on the
observationally-derived key quantities ($M_{\rm BH}$, $L_{\rm bol}$, and
$\lambda_{\rm E}$), and selection effects. These distributions are
fundamentally interlinked through the ensemble version of Eq. 1; that is, the
bolometric LF can be expressed as the convolution of the BHMF and the ERDF.

Compared to the AGN LF, determining the BHMF and the ERDF is much more
challenging. First, it requires reliable $M_{\rm BH}$ and $\lambda_{\rm E}$
measurements for large, unbiased AGN samples. In practice, beyond the local
Universe this is only possible for unobscured, broad-line AGNs, thanks to
“virial” mass prescriptions calibrated against reverberation mapping
experiments (e.g., Shen, 2013; Peterson, 2014). Furthermore, certain selection
effects have to be taken into account, which go beyond the more common biases
affecting the LF (i.e., the Eddington bias; Eddington, 1913). For flux limited
surveys, the dominant effect is a bias against low mass and low Eddington-
ratio AGN. To address this bias, and others, the selection function of the
sample at hand has to be well understood, and both the BHMF and ERDF have to
be constructed (and/or fitted) simultaneously.

With these challenges in mind, the BHMF (both active and total) and ERDF have
previously been constrained from large surveys, both indirectly and directly.
First, by assuming a universal relation between total stellar mass and BH mass
(see, e.g., Marconi et al., 2004; Sani et al., 2011a; Kormendy & Ho, 2013),
the shape of the total BHMF can be associated to the galaxy stellar mass
function and determined empirically. Since the AGN LF can be expressed as a
convolution of the BHMF and the ERDF, once the shape of the BHMF is known (or
assumed), then the LF can be used to constrain the ERDF indirectly (e.g.,
Caplar et al., 2015; Weigel et al., 2017). Other studies have used the ratio
between bolometric AGN luminosity and stellar mass as an indirect proxy for
$\lambda_{\rm E}$ and thus the ERDF (e.g., Aird et al., 2018a; Georgakakis et
al., 2017).

In contrast to such indirect approaches, the active BHMF and the ERDF can also
be determined directly from observations AGN for which it is possible to get
reliable SMBH mass measurements. Greene & Ho (2007, 2009) constrained a low-
redshift active BHMF for broad-line AGNs drawn from the Sloan Digital Sky
Survey (SDSS), focusing on relatively low-mass systems, and using the
1/$V_{\mathrm{max}}$ method to estimate the space densities (Schmidt, 1968).
On the other hand, Schulze & Wisotzki (2010) determined the active BHMF and
ERDF of highly luminous low-redshift quasars, drawn from the Hamburg-ESO
survey. The 1/$V_{\mathrm{max}}$ method was used to determine the BHMF of
$0.3{\lesssim}z{\lesssim}5$ quasars drawn from the SDSS/DR3 (Vestergaard et
al., 2008) and L/BQS surveys (Vestergaard & Osmer, 2009). The BHMF of quasars
of comparably high luminosities and redshifts from the BQS and SDSS was also
determined by Kelly et al. (2009) and Kelly & Shen (2013), respectively, with
the latter also constraining the ERDF. Nobuta et al. (2012) reported the BHMF
and ERDF of broad-line AGN at $z=1.4$ selected from the Subaru XMM-Newton Deep
Survey (SXDS) field. Schulze et al. (2015) constrained the BHMF and ERDF for
AGNs in the $1<z<2$ range, combining data from the SDSS, zCOSMOS and VVDS
optical surveys. Assuming a luminosity-dependent fraction of obscured AGN,
Schulze et al. (2015) also indirectly deduced the BHMF and ERDF of obscured
(Type 2) AGNs.

Several of these studies addressed some limitations of the
1/$V_{\mathrm{max}}$ method, and generally showed that the BHMF can be
described by a (modified) Schechter-like functional form (Aller & Richstone,
2002), resembling the shape of the galaxy stellar mass function. The ERDF, on
the other hand, is often described using a broken-power-law shape (Caplar et
al., 2015; Weigel et al., 2017), sharply decreasing towards high $\lambda_{\rm
E}$, and rarely exceeding the nominal Eddington limit.

High energy X-rays are understood to be more suited for probing large samples
of obscured, narrow-line (Type 2) AGNs, than the optical and UV bands due to
the penetrating power of X-ray photons through high column densities of gas
and dust. X-rays are also considered a higher-purity tracer of AGN than the
infrared (IR) band, as X-rays arising from SMBH accretion are less
contaminated by stellar and gas emission originating from the host galaxy.
Surveys such as SDSS are highly complete in optical bands, however they (as
well as soft X-ray surveys) are naturally biased towards unobscured AGN [i.e.,
$\log(N_{\rm H}/{\rm cm}^{-2})$$<22$]. This, combined with the inability to
measure $M_{\rm BH}$ and/or $\lambda_{\rm E}$ in narrow-line (mostly obscured)
AGNs, means that essentially all existing literature presents BHMFs and ERDFs
for only broad-line, unobscured (Type 1) AGNs. Nobuta et al. (2012) reported
BHMF/ERDF for an X-ray selected sample, but only for broad-line AGN. Aird et
al. (2018b) constrained the distribution of $\lambda_{\rm E}$, as probed
indirectly by the X-ray luminosity to host mass ratio, out to $z\simeq 4$,
using a large (Chandra) X-ray-selected sample, also concluding that the AGN
population tends towards higher $\lambda_{\rm E}$ with increasing redshift,
but not exceeding the Eddington limit. However, that study did not correct for
obscuration, as it was assumed that obscuration does not significantly affect
a hard X-ray-selected (2–7 keV) sample. Probing the key distributions of the
AGN population using a large sample selected in the ultra-hard X-ray regime
thus offers a crucial addition to our understanding of SMBH accretion and
triggering, even in the local Universe.

The Swift/BAT AGN Spectroscopic Survey (BASS, originally presented in Koss et
al. 2017) provides a large, highly complete sample of AGN selected in the
ultra-hard X-ray band (14–195 keV), along with reliable measurements of key
properties. This includes X-ray fluxes and column densities derived from
detailed X-ray spectral analysis (Ricci et al., 2017a), and—crucially—optical
counterpart matching, redshift measurements, and $M_{\rm BH}$ and
$\lambda_{\rm E}$ measurements, obtained through extensive optical
spectroscopic observations and analysis. As shown in Figure 2 of Koss et al.
(2016a), in terms of obscuration, BASS is the least biased of all X-ray
surveys to date, and largely unaffected by obscuring column densities below
$\log(N_{\rm H}/{\rm cm}^{-2})$$\simeq 23$. The 2nd data release of BASS
(BASS/DR2; Koss et al., subm.a) provides an essentially complete census of
luminosity, BH mass, Eddington-ratio, and obscuration towards all AGNs in the
70-month catalog of the Swift/BAT all-sky survey, with over 800 unbeamed AGNs,
mostly at $z\lesssim 0.3$. BASS/DR2 therefore allows us to determine the XLF,
BHMF and ERDF of low-redshift AGNs in an unprecedented way. It also presents
the first sample of ultra-hard X-ray selected AGN large enough to allow direct
measurements of the BHMF and the ERDF for both Type 1 and Type 2 AGN.
Importantly, the BH masses for Type 1 and Type 2 AGN are derived through
different, though consistent & inter-calibrated methods: broad Balmer line
measurements and “virial” prescriptions are used for the former, while stellar
velocity dispersion ($\sigma_{\star}$) measurements and the $M_{\rm
BH}-\sigma_{\star}$ relation is used for the latter. Highly complete (95% mass
measurement completeness at z $<0.3$) and with a well-understood selection
function, BASS thus allows us to gain a new understanding of the local AGN
population, while accounting for potential biases and statistical inference
limitations (e.g., see Schulze & Wisotzki, 2010; Kelly & Shen, 2013; Schulze
et al., 2015). With the intrinsic, bias-corrected XLF, BHMF, and ERDF at hand
for both Type 1 and Type 2 AGN, we can also investigate trends with column
density, AGN duty cycle, $M_{\rm BH}$ and $\lambda_{\rm E}$, as well as
whether the AGN BHMF is consistent with galaxy-SMBH scaling relations.

We present our work as follows: in §2 we describe the data and selection
criteria, in §3 we discuss the details of our method of calculating the BHMF,
ERDF and the XLF, in §4 we present the results of our analysis and in §5 we
discuss the physical implications of our results. Further details about our
methods, such as estimating errors and testing our code with mock catalogs are
described in the Appendices. A flat $\Lambda$CDM cosmology with $H_{\rm
0}=70\,{\rm km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{\rm M}=0.3$, and
$\Omega_{\rm\Lambda}=0.7$ is assumed throughout the paper. All uncertainties
reported in this paper are $\pm 1\sigma$ from the best-fit values.

##  2 Sample, Data and Basic Measurements

###  2.1 Sample selection

####  2.1.1 Input catalog and sample

Figure 1:  Venn diagram of all Swift/BAT 70-month sources below $z=0.3$, after
removing all beamed sources, Galactic-plane sources ($|b|<5^{\circ}$), as well
as more subtle cuts that exclude a total of 181 sources (see text for
details). The diagram shows how the final samples of 366 Type 1 objects and
220 Type 2 objects were selected. All cuts applied to the parent sample are
explicitly listed in Table 1. For each of the AGN sub-samples, the red circle
denotes the set of objects that fall within $0.01\leq z\leq 0.3$, the green
circle denotes the set that falls within $6.5\leq\log(M_{\rm
BH}/M_{\odot})\leq 10.5$ and the black circle denotes the set that falls
within $-3\leq\log\lambda_{\rm E}\leq 1$ . The numbers of sources in all sets
and intersections of sets are also noted.

Figure 2:  Swift-BAT cumulative number counts (left) and differential number
counts (right) from 70-month survey with beamed sources, low flux dual
sources, sources close to the Galactic plane and weak associations removed
(672 sources, at $z\leq 0.3$; turquoise data points), the same sample with an
additional redshift cut applied ($z\geq 0.01$, 619 sources; yellow data
points), and then with all Table 1 cuts applied (586 sources; red data
points). The latter sample is used in our BHMF and ERDF calculations, and as
seen in this figure, it is representative of the redshift restricted sample.
We produce fits to these data using Equation 6, and report the slope and
normalization of the fitted line in §2.3. To produce the fits, we applied a
flux cut at $F_{\rm 14-195,\leavevmode\nobreak\ obs}>10^{-11.1}\,{\rm
erg\,s^{-1}\,cm^{-2}}$ (shown in the left panel with black dashed line),
because the flux-area curve below that flux limit is not well constrained (§3
and Fig. 3), artificially affecting the slope of the number counts (shown in
the top left corner of the cumulative counts).

We base our analysis on the second data release (DR2) of BASS, which is
described in detail in Koss et al., (subm.a). BASS/DR2 combines extensive
optical spectroscopic measurements, X-ray spectral analysis, and derived
quantities, for essentially all 838 X-ray detected sources that are part of
the 70-month Swift-BAT catalog (Baumgartner et al., 2013; Ricci et al.,
2017a). We further discuss aspects of this catalog and the BASS/DR2 parent
sample that are crucial for our analysis, such as the flux limit (or rather,
the flux-area curve), in subsequent Sections.

For the X-ray related properties of the vast majority of these AGN, we rely on
the detailed spectral measurements described in Ricci et al. (2017a). BASS/DR2
also includes 14 AGN that were robustly detected as ultra-hard X-ray sources
in the 70-month Swift/BAT catalog, but not identified as AGN in the Ricci et
al. (2017a) compilation. The X-ray spectra of these source are analyzed as
part of BASS/DR2, following the same procedures as in Ricci et al. (2017a). A
detailed account of these newly identified AGN, as well as on a few more minor
changes to counterpart matches, can be found in Koss et al., (subm.b).

####  2.1.2 Excluded sources

First, we exclude all sources with Galactic latitudes $|b|<5^{\circ}$, as the
reliability of cross-matching BAT sources with optical counterparts, as well
as the completeness of the BASS optical spectroscopy efforts, drop
significantly for sources in the Galactic plane. All our survey area and/or
volume coverage calculations are adjusted to reflect the exclusion of this
region of the sky.

We also exclude from our analysis 105 beamed sources (i.e., blazars — BL-Lac-
like or FSRQ sources). More details about these sources are provided in Koss
et al., (subm.a). Here we briefly mention that this set constitutes of sources
identified through the BZ_flag in Ricci et al. (2017a), sources identified in
a dedicated BAT-Fermi analysis (Paliya et al., 2019), and a few additional
sources for which extensive multi-wavelength data suggests that they are most
likely beamed Paiano et al. (2020); Marcotulli+ (in prep.). After making these
two adjustments, $713$ non-beamed, non-Galactic-plane AGNs remains in our
sample.

We further exclude several dual AGN systems. These dual ultra-hard X-ray
emitting systems are identified in the Swift/BAT 70-month catalog thanks to
the combined emission of the dual sources but are too faint to individually
fall above the flux limit (see, e.g., the earlier Swift/BAT study by Koss et
al. 2012). Such sources should be removed from our analysis, as it is
fundamentally based on a well-defined, flux-limited sample. Details about each
of these dual AGN systems are provided in Koss et al., (subm.a). Of the 10
sources in dual systems in BASS/DR2, only NGC 6240S (BAT ID 841; Puccetti et
al., 2016) and MCG+04-48-002 (BAT ID 1077; Koss et al., 2016b) fall above our
nominal the flux cut.

An additional 26 sources that are detectable due to their flux being enhanced
through blending with a nearby brighter BAT source are also excluded. These
unassociated faint X-ray sources are described in detail in Ricci et al.
(2017a) and Koss et al., subm.a. Together, the exclusion of faint/weak
associations and of dual sources removes 37 objects from the sample, and
leaves us with of $678$ AGN.

####  2.1.3 Final, redshift-restricted AGN samples

After applying all the cuts mentioned so far, our base sample of non-beamed,
non-Galactic-plane AGN includes 678 sources, and is used for (parts) of our
XLF analysis. We next introduce a number of additional criteria to retain only
those AGN 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\leq z\leq 0.3$ range.
This excludes the 53 nearest BASS/DR2 AGN, 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\leq z\leq 0.3$ sample contains 619 objects (after excluding 6 more AGN
at at $z>0.3$).

We next limit our sample to have reliable measurements of $L_{\rm bol}$,
$M_{\rm BH}$, and $\lambda_{\rm 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 AGN with black hole masses in the
range $6.5\leq\log(M_{\rm BH}/M_{\odot})\leq 10.5$, and with Eddington ratios
in the range $-3\leq\log\lambda_{\rm E}\leq 1$.

Throughout our analysis, we further classify sources as being either broad-
line (“Type 1” hereafter) or narrow-line (“Type 2” hereafter) AGN, based on
the presence of any broad Balmer lines (as described in Koss et al.,, subm.a,
subm.b) and on how their most reliable measures of $M_{\rm BH}$ and
$\lambda_{\rm E}$ were derived (see Section 2.2 below), which in turn has
consequences for uncertainty and incompleteness estimation. 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$\alpha$ but no broad H$\beta$. Our Type 2 AGNs (220 sources) have only
narrow (Balmer) emission lines. Type 1 AGN are relatively unobscured [e.g.,
$\log(N_{\rm H}/{\rm cm}^{-2})\leq 22$] whereas Type 2 AGN tend to be more
heavily obscured [$\log(N_{\rm H}/{\rm cm}^{-2})>22$], although as shown in
Section 3, these $\log(N_{\rm H}/{\rm cm}^{-2})$ limits do not apply strictly
(see also Oh et al., submitted).

For our BHMF and ERDF analysis, we exclude the four Seyfert 1.9 sources that
have mass estimates only from broad H$\alpha$ 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., (submitted). One of these
four sources (ID 476) is the only AGN in our sample that has an estimated
Eddington ratio formally greater than 10 ($\log\lambda_{\rm E}\simeq 1.4$).
Therefore, the $\lambda_{\rm 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 mass greater than the
upper limit on $M_{\rm BH}$. These two upper limits are important for
computational reasons, and are discussed in § 3.4.

The Venn diagrams in Figure 1 show the number of non-beamed BASS/DR2 AGN that
meet each of our selection criteria (in $z$, $M_{\rm BH}$, and/or
$\lambda_{\rm E}$), as specified in Table 1, for both Type 1 and Type 2 AGN.
We also show the number of sources that fall outside all these criteria. Our
main BHMF and ERDF analysis is done using the 586 AGN that meet all criteria
(366 Type 1 and 220 Type 2 AGN), as shown in the central regions of the Venn
diagrams.

###  2.2 Data and Basic measurements

The present analysis is based on four basic, interlinked measurements that are
available for the BASS/DR2 AGN thanks to the extensive X-ray and optical
spectroscopy that is in the heart of the BASS project: bolometric luminosities
($L_{\rm bol}$), black hole masses ($M_{\rm BH}$), Eddington-ratios
($\lambda_{\rm E}$), and the line-of-sight hydrogen column densities ($N_{\rm
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.
One key consideration for deriving these quantities is to adopt a uniform
approach to all BASS AGN, whenever possible (e.g., for $L_{\rm bol}$), and to
revert to differential methodologies only when absolutely needed (e.g., for
$M_{\rm BH}$ determination).

First, as a primary probe of the AGN luminosity we use the integrated,
intrinsic luminosity between 14–195 keV ($L_{14-195\,{\text{keV}}}$ or $L_{\rm
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_{\rm H}$
measurements we use here, and we stress that the intrinsic
$L_{14-195\,{\text{keV}}}$ already accounts for the line-of-sight obscuration
(as probed by $N_{\rm H}$). The measurement uncertainties on the X-ray
luminosities we use are rather small, not exceeding ${\sim}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_{\rm bol}$, were then derived directly from
$L_{14-195\,{\text{keV}}}$, by using a simple, universal scaling of

$L_{\rm bol}=\kappa_{14{-}195\,{\text{keV}}}\times L_{\rm
14-195\,{\text{keV}}}\,.$ (2)

As explained in other key BASS publications (e.g., Koss et al.,, subm.b), this
bolometric correction of $\kappa_{14{-}195\,{\text{keV}}}=7.4$ corresponds to
a lower-energy bolometric correction of $\kappa_{2{-}10\,{\text{keV}}}=20$,
which is the average correction found for the BASS sample following the
$L_{2-10\,{\text{keV}}}$-dependent prescription of Marconi et al. (2004), and
further assumes the median X-ray power-law index found for the BASS sample,
$\Gamma=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\,{\text{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_{\rm bol}$ for even the most obscured AGN [i.e., Compton-thick sources
with $\log(N_{\rm H}/{\rm 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 were suggested and used in the literature, including corrections
that depend on luminosity, $\lambda_{\rm 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 therin). 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.

Black hole masses, $M_{\rm BH}$, are determined using two different approaches
for AGN with or without broad Balmer line emission, and specifically broad
H$\alpha$ line emission, which allows for a certain $M_{\rm BH}$ estimation
procedure (see immediately below). As noted above, for simplicity we refer to
such sources simply as Type 1 and Type 2 sources, respectively, however we
note that their detailed classification may be more nuanced (see Koss et al.,
subm.a; Mejía-Restrepo et al., submitted for a detailed discussion).

For broad line (Type 1) AGN, we rely on detailed spectral decomposition of the
H$\alpha$ spectral complex, and “virial” BH mass prescriptions which are
calibrated against reverberation mapping experiments, as described in detail
in Mejía-Restrepo et al., (submitted). Specifically, Mejía-Restrepo et al.,
(submitted) followed the prescription provided by Greene & Ho (2005, Eq. 6),
but adjust the virial factor to $f=1$, yielding:

$M_{\rm BH}=2.67\times 10^{6}\,\left(\frac{L[{\rm
b}{\text{H}}\alpha]}{10^{42}\,{\text{erg\,s}}^{-1}}\right)^{0.55}\left(\frac{{\rm
FWHM}[{\rm b}{\text{H}}\alpha]}{1000\,{\rm
km\,s}^{-1}}\right)^{2.06}M_{\odot}\,,$ (3)

where $L[{\rm b}{\text{H}}\alpha]$ and ${\rm FWHM}[{\rm b}{\text{H}}\alpha]$
are the luminosity and width of the broad part of the H$\alpha$ emission line,
respectively. For narrow line (Type 2) AGN, Koss et al. (2017) and Koss et
al., (subm.a) rely on measurements of the stellar velocity dispersion
($\sigma_{\star}$) in the AGN host galaxies and the well-known $M_{\rm
BH}-\sigma_{\star}$ relation. Specifically, $\sigma_{\star}$ was measured from
either the Ca H+K, Mg i, and/or the Ca ii triplet spectral complexes, as
described in detail in Koss et al., (in prep.). We then used the relation
determined by Kormendy & Ho (2013, Eq. 3):

$M_{\rm BH}=3.09\times 10^{8}\left(\frac{\sigma_{\star}}{200\,{\rm
km\,s}^{-1}}\right)^{4.38}M_{\odot}.$ (4)

We stress that these two types of $M_{\rm BH}$ prescriptions are consistently
calibrated. Specifically, the broad-line $M_{\rm BH}$ prescription in Eq. 3 is
derived by assuming that low-redshift, broad-line AGNs (particularly those
with reverberation mapping measurements) lie on the same $M_{\rm
BH}-\sigma_{\star}$ relation as do narrow-line AGNs and inactive galaxies
(that is - Eq. 4; see, e.g., Onken et al. 2004 and Woo et al. 2013 for
detailed discussion).

The uncertainties in both types of $M_{\rm BH}$ estimates are completely
dominated by systematics and are of the order $0.3-0.5$ dex (Gültekin et al.,
2009; Shankar et al., 2019; Shen, 2013; Peterson, 2014). The lower end of this
range is consistent with the scatter seen in the $M_{\rm BH}-\sigma_{\star}$
(or $M_{\rm BH}-M_{\rm host}$) relation (e.g., Sani et al., 2011b; Kormendy &
Ho, 2013), while the upper end includes also the other key ingredients of
“virial” $M_{\rm BH}$ determinations in broad-line AGNs (e.g., Shen, 2013;
Peterson, 2014). The uncertainties on $\lambda_{\rm E}$ are naturally
dominated by the systematic uncertainties on $M_{\rm BH}$, and are thus also
of the order $0.3-0.5$ dex, although the range and trends seen for bolometric
corrections (e.g., Marconi et al., 2004; Vasudevan & Fabian, 2009) suggest
that uncertainties on $\lambda_{\rm E}$ may be yet higher. In comparison, our
measurement uncertainties on $L({\rm b}{\text{H}}\alpha)$, FWHM(bH$\alpha$),
and $\sigma_{\star}$ are typically of order 10%, which would add up to $\sim$
0.1 dex uncertainty for $M_{\rm BH}$ estimates in broad-line AGN, and $<$ 0.2
dex for narrow-line AGN. Our analysis takes into account the large
(systematic) uncertainties on $M_{\rm BH}$ determinations for all BASS AGN, as
detailed in §3.

Finally, $\lambda_{\rm E}$ is determined by combining the $L_{\rm bol}$ and
$M_{\rm BH}$ estimates mentioned above, using the relation

$\log\lambda_{\rm E}=\log(L_{\rm bol}/{\text{erg\,s}}^{-1})-\log(M_{\rm
BH}/M_{\odot})-38.18\,,$ (5)

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

There are alternative galaxy-black hole scaling relationships suggested by
Bernardi et al. (2007), Shankar et al. (2017) and Shankar et al. (2020), that
correct for selection biases that may affect the samples used for calibrating
these mass measurement relationships. While recalculating 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.

###  2.3 Source number counts

Figure 2 shows the cumulative and differential source counts (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\leq 0.3$; 672 sources), (2) our redshift-restricted AGN sample
($0.01\leq z\leq 0.3$; 619 sources), and (3) our final BASS/DR2 BHMF/ERDF
sample (i.e., with $M_{\rm BH}$ and $\lambda_{\rm E}$ cuts; 586 sources). All
three samples exclude the sources described in §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:

$\frac{dN}{dS}=A\times(S/10^{-11})^{-\alpha}\,.$ (6)

We limit our fits to the $F_{\rm 14-195,\leavevmode\nobreak\ obs}\geq
10^{-11.1}\,{\rm erg\,s^{-1}\,cm^{-2}}$ flux regime, since the flux-area curve
for Swift-BAT 70-month survey (discussed in more detail in §3) is sparsely
sampled at lower flux level, making it difficult to accurately calculate the
surface density of faint BAT sources (see left panel of Fig. 2).  Our fits are
derived using orthogonal distance regression (ODR), so to properly account for
uncertainties on both axes (i.e., $\log S$ and $\log N$). One limitation of
the ODR method is the inherent assumption of symmetric uncertainties (which in
our case are an average 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 choice of fitting method
(and 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).

The best-fit slopes derived using ODR are: $\alpha=2.50\pm 0.04$, $2.62\pm
0.07$ and $2.67\pm 0.05$, for the 70-month AGN, the redshift-restricted, and
final BHMF/ERDF AGN samples, respectively. The best-fit normalizations are
$\log A=9.43\pm 0.03$, $9.44\pm 0.03$, and $9.43\pm 0.02$, respectively. The
MLE-based best fitting results are in agreement with the ORD ones, within $\pm
1\sigma$ errors. Our best-fit curve for the full sample is consistent with the
expected slope for a fully uniform, Euclidean distribution ($\alpha=2.5$), and
also with the results of several previous studies of ultra-hard X-ray selected
AGN 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
expected as objects have been removed from the total sample. We conclude that
our BASS/DR2 AGN sample(s), 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 AGN. We are
thus confident that we can rely on the same selection criteria, and
specifically sky coverage curves, as derived for the parent Swift/BAT 70-month
catalog.

##  3 Statistical Inference Methods

Table 1: Overview of sample selection and analysis parameters.a

Quantity/variable | Symbol/value
---|---
Minimum redshift considered | $z_{\rm min,s}=0.01$
Maximum redshift considered | $z_{\rm max,s}=0.3$
Minimum black hole mass considered | $\log(M_{\rm BH,min,s}/M_{\odot})=6.5$
Maximum black hole mass considered | $\log(M_{\rm BH,max,s}/M_{\odot})=10.5$
Minimum Eddington-ratio considered | $\log\lambda_{\rm E,min,s}=-3$
Maximum Eddington-ratio considered | $\log\lambda_{\rm E,max,s}=1$
Luminosity bin size for $V_{\mathrm{max}}$ method | $d\log L_{\rm X}=0.3$
BH mass bin size for $V_{\mathrm{max}}$ method | $d\log(M_{\rm BH}/M_{\odot})=0.3$
Edd. ratio bin size for $V_{\mathrm{max}}$ method | $d\log\lambda_{\rm E}=0.3$
Assumed uncertainty on $\log M_{\rm BH}$ |  $\sigma_{M_{\rm BH}}=$ 0.3, 0.5 dex
Assumed uncertainty on $\log\lambda_{\rm E}$ |  $\sigma_{\log\lambda_{\rm E}}=$ 0.3b, 0.5 dex
Galactic plane exclusionc | $|b|\geq 5^{\circ}$
Other cuts | Beamed AGN are excluded

aafootnotetext: 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 main text).

bbfootnotetext: We also investigate the effect of an uncertainty of 0.2 dex in
luminosity due to measurement uncertainty and the X-ray bolometric correction
(Eq. 2). Variable bolometric corrections are explored in Appendix E.

ccfootnotetext: Meant to guarantee a high completeness of optical counterpart
identification and spectroscopic coverage.

The main goal of the present study is to determine and interpret the intrinsic
distributions of X-ray luminosity ($L_{14-195\,{\text{keV}}}$), SMBH mass
($M_{\rm BH}$), and Eddington-ratio ($\lambda_{\rm E}$)—namely the X-ray
luminosity function (XLF), black hole mass function (BHMF), and Eddington-
ratio distribution function (ERDF)—for AGN in the present-day Universe, in the
most complete way possible. In what follows, we provide a detailed description
of the statistical inference methods we use to derive these distributions from
the basic measurements available for our sample of AGN drawn from the all-sky
Swift-BAT survey and the BASS project.

The main obstacle in deriving the statistical properties of any sample of
astrophysical sources is to account for the various factors of incompleteness
and bias that are encoded in the observed sample in hand. The first and most
obvious source of incompleteness for a flux-limited survey such as that of
Swift-BAT AGN is the Malmquist bias, where less luminous sources can only be
detected within a small volume at low redshift, whereas higher-luminosity ones
are detected even at high $z$. This leads to the severe underestimation of the
space densities of the former and an overestimation of the space densities of
the latter.

More complex forms of bias are introduced once the incompleteness in terms of
luminosity is translated into incompletenesses in the distributions of related
quantities, such as masses or growth rates (e.g., Marchesini et al., 2009;
Pozzetti et al., 2010; Weigel et al., 2016). Specifically for the present
study, the BASS AGN are selected based on their ultra-hard X-ray AGN
luminosity. For the BHMF, this results in a bias against low mass black holes
which will be too faint to lie above the flux limit, unless they have
exceptionally high Eddington-ratios. Similarly, low Eddington-ratio AGN will
not be part of the ERDF since they are too faint, even if they are massive.
Schulze & Wisotzki (2010, hereafter SW10) refer to this incompleteness as
“sample censorship”, however we refer to it as “sample truncation” hereafter.

Another potential source of bias may arise from the need to measure the
properties for which the statistics are to be surveyed. Specifically, our
estimates of BH mass depend on robustly measuring the luminosities and widths
of broad emission lines (Eq. 3) or the widths of stellar absorption features
(Eq. 4). This, in turn, requires well-calibrated, medium-to-high spectral
resolution observations, and also carries significant (systematic)
uncertainties.

Our statistical inference methodology accounts for all these possible biases.
We first constrain the XLF by using the classical 1/$V_{\mathrm{max}}$ method
(Schmidt, 1968; Avni et al., 1980). We apply the same method to the $M_{\rm
BH}$ and $\lambda_{\rm E}$ measurements to gain an initial guesses for the
BHMF and the ERDF. We assume functional forms and fit the resulting
distributions. We then use a parametric maximum likelihood approach and our
initial guess to simultaneously correct the BHMF and the ERDF for sample
truncation. For the bias correction we follow the approach by SW10 and Schulze
et al. (2015, hereafter S15). We test our approach by creating mock catalogs,
determining the corresponding distributions, and comparing those to the
assumed input.

Throughout this work, we implicitly assume that the ERDF is $M_{\rm
BH}$-independent. This assumption is further justified by the analysis we
present, based on splitting our sample to two $M_{\rm BH}$ regimes (Section
4). We also do not impose any Eddington-ratio ($\lambda_{\rm E}$) dependence
of the distribution of absorbing columns (i.e., $N_{\rm H}$). Although a
certain, complex, link between $\lambda_{\rm E}$ and $N_{\rm H}$ was suggested
by several studies (Ricci et al., 2017b, and references therein), our choice
allows us to obtain independent evidence for such a link using our Type 1
(mostly unobscured) and Type 2 (mostly obscured) AGN samples, instead of
assuming it a priori.

Figure 3: Flux-area curve of the Swift/BAT 70-month all-sky survey. To compute
the XLF and gain an initial guess for the BHMF and the ERDF, we use this flux-
area curve rather than a constant flux limit. For bright ultra-hard X-ray
sources, $\log(F_{\rm 14{-}195,obs}/{\rm erg\,cm}^{-2}\,{\rm s}^{-1})>-10.7$,
the survey is complete over the entire sky. For fainter sources, the
completeness decreases as a function of the observed 14$-$195 keV X-ray flux.
This Figure converts Figure 10 of Baumgartner et al. (2013) cumulative flux-
area curve from mCrab to erg cm-2 s-1units.

###  3.1 Survey specific considerations

Figure 4:  Hard X-ray (14$-$195 keV) attenuation curve for the BAT survey,
assuming two torus opening angles: 60∘ (solid blue line) and 35∘ (dotted blue
line). For extremely obscured sources, with $\log(N_{\rm H}/{\rm
cm}^{-2})>25$, the fraction of observed luminosity relative to the intrinsic
luminosity is less than 5%, making them very difficult to detect even with the
high energy X-ray window of BAT. This is an update to the Ricci et al. (2015)
attenuation curve with a newer torus model (see text for details).

We take advantage of the well-constrained flux-area curve of the Swift-BAT
survey (Baumgartner et al., 2013). The flux-area curve $\Omega_{\rm sel}(\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 $\log(F_{\rm 14{-}195,obs}/{\rm
erg\,cm}^{-2}\,{\rm s}^{-1})>-10.7$, the BAT survey is complete over the
entire sky, i.e. $\Omega_{\rm sel}=1$. Sources with $\log({\rm F}_{\rm
X,obs}/{\rm erg\,cm}^{-2}\,{\rm s}^{-1})<-11.3$ are completely missed by the
BAT 70-month survey & catalog, and thus $\Omega_{\rm sel}=0$. Our analysis
makes explicit use of the complete flux-area curve, including the intermediate
values for sources with $-11.3\lesssim\log(F_{\rm 14{-}195,obs}/{\rm
erg\,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 non-
uniform (as shown in Figure 1 of Baumgartner et al., 2013). Taking the slight
positional differences in sensitivity is beyond the scope of this work.

Additionally, we account for the fact that the XLF describes the distribution
of AGN 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 observed ultra-hard X-ray luminosities, we use the Swift/BAT
attenuation curve shown in Figure 4. This observed luminosity is a sum of
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), 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
$\Gamma=1.8$, ${\rm E}_{\rm cutoff}=200$ keV, a torus opening angle of
$60^{\circ}$ and an inclination angle of $72^{\circ}$ (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 torus
opening angle and (line-of-sight) inclination angle. Following Ricci et al.
(2015), we also consider another attenuation curve for a model with opening
angle of $35^{\circ}$ and an inclination angle of $72^{\circ}$ (shown in
Figure 4). We discuss the impact of varying our template spectra (and other
model dependent parts of our approach) on our final results in §3.4.2.

While the BAT 14$-$195 keV energy range is not strongly affected by
attenuation up to $\log(N_{\rm H}/{\rm cm}^{-2})\simeq 23.5$, Fig. 4 shows
that for Compton-thick sources [$\log(N_{\rm H}/{\rm cm}^{-2})\gtrsim 24$] at
least 30% and up to $\sim$98% of the intrinsic ultra-hard X-ray emission is
lost [i.e. only the scattered component is detectable as $\log(N_{\rm H}/{\rm
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_{\rm H}$ of each AGN (Ricci et
al., 2017a; Koss et al.,, subm.a).

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
ERDF, as the three key properties ($L_{14-195\,{\text{keV}}}$, $M_{\rm BH}$,
$\lambda_{\rm E}$) are closely related, through Eqs. 2 and 5.

###  3.2 The 1/$V_{\mathrm{max}}$ method and the XLF

The 1/$V_{\mathrm{max}}$ method (Schmidt, 1968) provides a way to correct for
the incompleteness (i.e., lower luminosity/mass objects fall below survey
sensitivity at larger distances). When counting the sources within a given
luminosity (or mass etc.) bin, each source $i$ is weighted by $V_{\rm max,i}$
– the maximum volume within which source $i$ could have been detected, given
the survey properties. In the case of a luminosity function, 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 AGN luminosities, black hole masses, and
Eddington-ratios of the BASS sample, we bin in $\log M_{\rm BH}$,
$\log\lambda_{\rm E}$, and $\log L_{\rm X}$.222Note that for simplicity, we
use $L_{\rm X}$ to denote $L_{14-195\,{\text{keV}}}$unless explicitly noted
otherwise. We then use the 1/$V_{\mathrm{max}}$ method to determine the
corresponding space densities $\Phi_{\rm M}(\log M_{\rm BH})$,
$\xi(\log\lambda_{\rm E})$, and $\Phi_{\rm L}(\log L_{\rm X})$. We compute the
$V_{\mathrm{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_{\mathrm{max}}$-based
calculation does not include a robust correction neither for sample truncation
nor for uncertainties on the relevant key quantities, it only allows us to
gain an initial guess for the BHMF and the ERDF.

For the XLF, the space density in luminosity bin $j$ is given by the sum over
all $N_{\rm bin}$ weighted objects within this bin:

$\Phi_{\rm L,j}\,d\log L_{\rm X}=\sum_{i}^{N_{\rm
bin}}\frac{1}{V_{\rm{max},i}}.$ (7)

We compute $\Phi_{\rm M}(\log M_{\rm BH})$ and $\xi(\log\lambda_{\rm E})$
similarly and use the bin sizes $d\log L_{\rm X}$, $d\log M_{\rm BH}$, and
$d\log\lambda_{\rm E}$ given in Table 1.

To compute $V_{\mathrm{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: $\Omega_{\rm sel}(\log L_{\rm X},N_{\rm H},z)$.
For each object $i$ at redshift $z_{i}$ with intrinsic luminosity $\log
L_{\mathrm{X},i}$ and column density $N_{\mathrm{H},i}$, the maximum comoving
volume $V_{\mathrm{max},i}$ is then given by (e.g., Hogg, 1999; Hiroi et al.,
2012):

$\displaystyle V_{\rm{max},i}(\log L_{\mathrm{X},i},N_{\mathrm{H},i},z_{i})=$
(8) $\displaystyle 4\pi\times\frac{c}{H_{0}}\times\int_{z_{\rm min,s}}^{z_{\rm
max,s}}\Omega_{\rm sel}(\log
L_{\mathrm{X},i},N_{\mathrm{H},i},z)\frac{D_{C}(z)^{2}}{E(z)}dz.$

In this expression, $D_{C}(z)$ corresponds to the comoving distance to
redshift $z$, $E(z)$ is given by $\sqrt{\Omega_{\rm
M}(1+z)^{3}+\Omega_{\rm\Lambda}}$. The integration limits $z_{\rm min,s}$ and
$z_{\rm max,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_{\mathrm{max},i}$. Instead, this information is
encoded in the flux-area curve, which ultimately provides the same outcome as
$\Omega_{\rm sel}=0$ when the source observed flux becomes too faint [i.e.,
$\log(F_{\rm 14{-}195,obs}/{\rm erg\,cm}^{-2}\,{\rm s}^{-1})\lesssim-11.3$].

Note that the 1/$V_{\mathrm{max}}$ method is very sensitive to the flux-area
curve, which is non-uniform (as described in §3.1). As a specific example, the
AGN in NGC 5283 (BAT ID 684), with $z=0.01036$ and $\log(F_{\rm
14{-}195,obs}/{\rm erg\,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 area flux curve
drops to zero. As a result, the $V_{\mathrm{max}}$ within which this object
can be detected above $z>0.1$ is very low, and this drives up the
corresponding 1/$V_{\mathrm{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 §3.4) is
not as dramatically affected by a single object as the direct
1/$V_{\mathrm{max}}$ approach. We show the effect of including this single
object on 1/$V_{\mathrm{max}}$ values in Appendix F, whereas the
1/$V_{\mathrm{max}}$ values shown in 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_{\mathrm{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_{\mathrm{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 to
estimate luminosity functions while correcting for sample truncation was
suggested by Kelly et al. (2008). However, most previous (X)LF works used the
1/$V_{\mathrm{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
uncertainties on key quantities (in all three distribution functions), we
chose to present the 1/$V_{\mathrm{max}}$ results, if only as a first order
(albeit somewhat biased) estimate of the XLF.

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

###  3.3 Functional forms for XLF, BHMF and ERDF

After having determined the XLF and having determined an initial guess for the
BHMF and the ERDF via the 1/$V_{\mathrm{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:

$\Phi_{\rm L}(\log L_{\rm X})=\frac{dN}{d\log L}=\Phi^{*}_{\rm
L}\times\left[\left(\frac{L_{\rm X}}{L^{*}_{\rm
X}}\right)^{\gamma_{1}}+\left(\frac{L_{\rm X}}{L^{*}_{\rm
X}}\right)^{\gamma_{2}}\right]^{-1}\,.$ (9)

For the fitting procedure we parametrize the second power-law slope as
$\gamma_{2}=\gamma_{1}+\epsilon_{\gamma}$ with $\epsilon>0$ to avoid
degeneracy between the two exponents. Given what is known about the rather
universal nature of the AGN LF (e.g., Shen et al., 2020, and references
therein), in practice this means that our $\gamma_{1}$ and $\gamma_{2}$
correspond to what is often referred to as the faint- and bright-end LF
slopes, respectively.

For the BHMF, $\Phi_{\rm M}(\log M_{\rm BH})\equiv\frac{dN}{d\log M_{\rm
BH}}$, we use a modified Schechter function, defined as:

$\displaystyle\Phi_{\rm M}(\log M_{\rm
BH})\propto\ln(10)\Phi^{*}\left(\frac{M_{\rm BH}}{\mathrm{M}^{*}_{\rm
BH}}\right)^{\alpha+1}\exp{\left[-\left(\frac{M_{\rm BH}}{\mathrm{M}^{*}_{\rm
BH}}\right)^{\beta}\right]}.$ (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:

$\xi(\log\lambda_{\rm E})=\frac{dN}{d\log\lambda_{\rm
E}}\propto\xi^{*}\times\left[\left(\frac{\lambda_{\rm E}}{\lambda_{\rm
E}^{*}}\right)^{\delta_{1}}+\left(\frac{\lambda_{\rm E}}{\lambda_{\rm
E}^{*}}\right)^{\delta_{2}}\right]^{-1}\,.$ (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-$\lambda_{\rm E}$ slope of
the ERDF (e.g., Caplar et al., 2015; Weigel et al., 2017). Similar to what was
done with the XLF, we parametrize $\delta_{2}$ as
$\delta_{1}+\epsilon_{\lambda}$ with $\epsilon_{\lambda}>0$, thus linking
$\delta_{1}$ and $\delta_{2}$ with the low- and high-$\lambda_{\rm E}$ ERDF
slopes, respectively.

In the initial part of our analysis, we fit the 1/$V_{\mathrm{max}}$
measurements of the BHMF and 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_{\mathrm{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_{\mathrm{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
1/$V_{\mathrm{max}}$ approach (discussed in § 3.2), where we fit the
aforementioned functions to the data directly.

###  3.4 Correcting the BHMF and the ERDF for sample truncation and for
uncertainties

####  3.4.1 The basic principle

To correct for sample truncation—i.e., the bias against low mass and low
Eddington-ratio AGN 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
should be considered as a parametric maximum likelihood approach.

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

$\displaystyle\Phi_{\rm M}(\log M_{\rm BH},z)=$
$\displaystyle\int_{\log\lambda_{\rm E,min,s}}^{\log\lambda_{\rm
E,max,s}}\Psi(\log M_{\rm BH},\log\lambda_{\rm E},z)d\log\lambda_{\rm E},$
(12) $\displaystyle\xi(\log\lambda_{\rm E},z)=$ $\displaystyle\int_{\log
M_{\rm BH,min,s}}^{\log M_{\rm BH,max,s}}\Psi(\log M_{\rm BH},\log\lambda_{\rm
E},z)d\log M_{\rm BH}.$

The integration limits $\log\lambda_{\rm E,min,s}$, $\log\lambda_{\rm
E,max,s}$, $\log M_{\rm BH,min,s}$, and $\log M_{\rm BH,max,s}$ correspond to
the minimum and maximum Eddington-ratios and black hole mass values that we
are considering in our analysis (see Table 1). The bivariate distribution
$\Psi(\log M_{\rm BH},\log\lambda_{\rm E})$ has the physical units of space
density, namely h3 Mpc-3 dex-2 (i.e., per dex in $M_{\rm BH}$ and per dex in
$\lambda_{\rm E}$). As the redshift range is very small, we consider only a
single redshift bin and do not model the redshift evolution of the XLF, the
BHMF, and/or the ERDF; this is why our bivariate distribution function
$\Psi(\log M_{\rm BH},\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 $\Psi(\log
M_{\rm BH},\log\lambda_{\rm E})$ (the expressions in Eq. 12) reduces to:

$\displaystyle\Phi_{\rm M}(\log M_{\rm BH})=$ $\displaystyle\tilde{\Phi}(\log
M_{\rm BH})\times\int_{\log\lambda_{\rm E,min,s}}^{\log\lambda_{\rm
E,max,s}}\tilde{\xi}(\log\lambda_{\rm E})d\log\lambda_{\rm E}$ (13)
$\displaystyle\xi(\log\lambda_{\rm E})=$
$\displaystyle\tilde{\xi}(\log\lambda_{\rm E})\times\int_{\log M_{\rm
BH,min,s}}^{\log M_{\rm BH,max,s}}\tilde{\Phi}(\log M_{\rm BH})d\log M_{\rm
BH}.$

The assumed functional forms for $\tilde{\Phi}(\log M_{\rm BH})$ and
$\tilde{\xi}(\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 $\Phi_{\rm M}(\log M_{\rm BH})$ and
$\xi(\log\lambda_{\rm E})$, we have to define the $M_{\rm BH}$ and
$\lambda_{\rm E}$ range that is being considered and marginalize over the
other distribution. As equation 13 shows, $\tilde{\Phi}(\log M_{\rm BH})$ and
$\tilde{\xi}(\log\lambda_{\rm E})$ retain their assumed shape, however their
normalizations are adjusted. In §3.4.4, we explain how these normalizations
are determined after constraining the functional forms.

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

$\ln\mathcal{L}=\sum_{i}^{N_{\rm obs}}\ln p_{i}(\log M_{\mathrm{BH},i},\log
N_{\mathrm{H},i},\log\lambda_{\mathrm{E},i},z_{i})\,,$ (14)

where $p_{i}$ is the probability of observing object $i$ with black hole mass
$\log M_{\mathrm{BH},i}$, Eddington-ratio $\log\lambda_{\mathrm{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_{\rm obs}$. For
an assumed $\Psi(\log M_{\rm BH},\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 which 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 $\Psi(\log M_{\rm BH},\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_{\rm obs}$ sources.

####  3.4.2 The probability of observing a given source

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

$\displaystyle p_{i}$ $\displaystyle(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log N_{\mathrm{H},i},z_{i})$
(15) $\displaystyle=\frac{N_{i}(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log
N_{\mathrm{H},i},z_{i})}{N_{\rm tot}}$ $\displaystyle=\frac{1}{N_{\rm
tot}}\Psi(\log M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i})$
$\displaystyle\leavevmode\nobreak\ \Omega_{\rm sel}(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log N_{\mathrm{H},i},z_{i})$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ p(\log
N_{\mathrm{H},i})\leavevmode\nobreak\ p(z_{i})\leavevmode\nobreak\
\frac{dV_{C}(z_{i})}{dz}$

Note that this expression by itself does not correct for measurement
uncertainties. We explain the mechanism that does ultimately account for these
uncertainties in §3.4.3 below. For simplicity, we assume that $p(z_{i})$ and
$p(\log N_{\mathrm{H},i})$ are independent of $\Psi$ 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}$:

  * •

$\Psi$, the intrinsic bivariate distribution function of BH mass and Eddington
ratio: For a given set of $\log M_{\mathrm{BH},i}$ and
$\log\lambda_{\mathrm{E},i}$ values, the bivariate distribution function
returns the space density of objects with those properties, i.e. the intrinsic
number of such AGN per unit comoving volume, per dex of BH mass and Eddington
ratio.

  * •

$dV_{C}/dz$, the comoving volume element: By multiplying the space density of
objects with $\log M_{\mathrm{BH},i}$ and $\log\lambda_{\mathrm{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):

$\frac{dV_{C}(z)}{dz}=4\pi\frac{c}{H_{0}}\frac{D_{C}(z)^{2}}{E(z)},$ (16)

where $D_{C}$ and $E(z)$ were defined in Section 3.2.

  * •

$N_{\rm 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 selection
effects have been taken into account. $N_{\rm tot}$ is given by an integral
that accounts for all the quantities that affect the observability of our AGN
within the survey:

$\displaystyle N_{\rm tot}=\iiiint$ $\displaystyle\Omega_{\rm sel}(\log
M_{\mathrm{BH}},\log\lambda_{\mathrm{E}},\log N_{\rm H},z)$ (17)
$\displaystyle\leavevmode\nobreak\ \Psi(\log M_{\rm BH},\log\lambda_{\rm
E})\leavevmode\nobreak\ p(\log N_{\rm H})\leavevmode\nobreak\ p(z)$
$\displaystyle\leavevmode\nobreak\ \frac{dV_{C}(z)}{dz}d\log M_{\rm
BH}\leavevmode\nobreak\ d\log N_{\rm H}\leavevmode\nobreak\ d\log\lambda_{\rm
E}\leavevmode\nobreak\ dz\,.$

Here, $d\log M_{\rm BH}$, $d\log\lambda_{\rm E}$, and $dz$ are computed over
the corresponding minimum and maximum values considered in the sample (see
Table 1).

  * •

$p(\log N_{\rm H})$: A bias-corrected intrinsic absorption function needs to
be assumed from which the $\log(N_{\rm H}/{\rm cm}^{-2})$ values of the
intrinsic AGN population are drawn. As explained above, while the intrinsic
luminosities of AGN are the product of their $M_{\rm BH}$ and $\lambda_{\rm
E}$, any selection or distribution function that depends on their observed
fluxes would also have to depend on (the distribution of) line-of-sight
obscuration, which we generally associate with circumnuclear (torodial) dusty
gas. Ricci et al. (2015) derived an intrinsic $\log N_{\rm H}$ distribution
specifically for Swift-BAT detected AGN considering X-ray reflection from a
torus. This model assumed an opening angle of 60∘ and considered two
luminosity bins - above and below $\log(L_{\rm
14-195\,{\text{keV}}}/{\text{erg\,s}}^{-1})=43.7$. This intrinsic $\log N_{\rm
H}$ distribution is shown in Figure 5.

To test the impact of different models, we also calculate 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 §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^{\circ}$) 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 $\log(N_{\rm H}/{\rm
cm}^{-2})=20-25$ in discrete bins of 1 dex width.

While calculating $N_{\rm tot}$ for all AGN, we assign $\log(N_{\rm H}/{\rm
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 chose to incorporate the luminosity dependent absorption
function as provided by Ricci et al. (2015), and do not impose any direct
$\lambda_{\rm E}$ dependence on the absorption function. While a $\lambda_{\rm
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 differences in the BHMF and 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 ERDF for Type 1 AGN, we assume
the intrinsic absorption function to be identical to the observed one, which
is justified given the negligible effect of absorption in such sources. For
Type 2 AGN, we subtract the distribution of Type 1 AGN from the bias-corrected
overall distribution in each luminosity bin. These distributions are shown
using 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(\log N_{\rm H})$
used in this work is in fact $p(\log N_{\rm H}|L_{\rm X}(M_{\rm
BH},\lambda_{\rm E}))$.

  * •

$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 numerator and
denominator of $p_{i}$).

  * •

$\Omega_{\rm sel}$: In equation 15, the term $\Omega_{\rm sel}(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log N_{\rm H,i},z_{i})$
corresponds to the selection function of the survey. In its simplest form,
$\Omega_{\rm sel}$ returns the value of the flux-area curve, $\Omega_{\rm
sel}(\log F_{\rm X})$, for each source (see Section 3.2). To predict the X-ray
flux $\log F_{\mathrm{X},i}$ for source $i$ with BH mass $\log
M_{\mathrm{BH},i}$, an Eddington ratio $\log\lambda_{\mathrm{E},i}$, a column
density $\log N_{\mathrm{H},i}$ and a redshift $z_{i}$, we perform the
following sequence of calculations.

First, we compute the corresponding bolometric luminosity, $L_{\rm bol}$,
using Eq. 5 (in Appendix E we experiment with a variable, luminosity-dependent
bolometric correction); from $L_{\rm bol}$, we then calculate the intrinsic
ultra-hard X-ray luminosity over the 14–195 keV range,
$L_{14-195\,{\text{keV}}}$, using Eqn. 2; from that we deduce the luminosity
and flux as measured by BAT, taking into account the column density (using the
curve shown in Fig. 4) and the luminosity distance to the source (using
$z_{i}$). The calculated $\log F_{\mathrm{X},i}$ is then used to obtain the
selection function based on the curve shown in Fig. 3.

Figure 5:  Top panel: Observed histograms of absorbing column density, $\log
N_{\rm H}$, for Type 1 (blue lines) and Type 2 (red lines) AGN, and the
overall distribution (black dotted lines). Bottom panel: Normalized
distributions observed for Type 1 (blue stacked bars) and Type 2 AGN (red
stacked bars), along with the bias-corrected column density distribution of
the BASS 70-month survey sample (Ricci et al., 2015), calculated for two
luminosity bins: $\log(L_{14-195\,{\text{keV}}}/{\text{erg\,s}}^{-1})\leq
43.7$ (black solid lines) and
$\log(L_{14-195\,{\text{keV}}}/{\text{erg\,s}}^{-1})>43.7$ (black dashed
lines). For Type 1 AGN (blue lines), we assume the intrinsic and the observed
column densities are the same, as Type 1 AGN tend to be relatively unobscured,
then we calculate the Type 2 distribution (red lines) by subtracting the Type
1 distribution from the overall distribution in each luminosity bin. We use
the intrinsic distributions when calculating the BHMF and ERDF for the sample.

####  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
uncertainty in each of the four underlying parameters. Specifically, the
selection effect corrected $p_{i}$ is as follows:

$\displaystyle p_{\rm i,conv}$ $\displaystyle(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log
N_{\mathrm{H},i},z_{i})=\frac{N_{\rm i,conv}}{N_{\rm tot}}$ (18)
$\displaystyle=\frac{1}{N_{\rm tot}}\iiiint\Psi(\log
M_{\mathrm{BH}},\log\lambda_{\mathrm{E}},\log N_{\mathrm{H}},z)$
$\displaystyle\ \times\Omega_{\rm sel}(\log
M_{\mathrm{BH}},\log\lambda_{\mathrm{E}},\log N_{\rm H},z)$ $\displaystyle\
\times p(\log N_{\rm H})\leavevmode\nobreak\ p(z)\leavevmode\nobreak\
\frac{dV_{C}(z)}{dz}$ $\displaystyle\ \times\omega(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log N_{\mathrm{H},i},z_{i}|\log
M_{\rm BH},\log\lambda_{\rm E},\log N_{\mathrm{H}},z)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ d\log M_{\rm
BH}\leavevmode\nobreak\ \leavevmode\nobreak\ d\log\lambda_{\rm
E}\leavevmode\nobreak\ \leavevmode\nobreak\ d{\log N_{\rm
H}}\leavevmode\nobreak\ \leavevmode\nobreak\ dz\,.$

Here the denominator $N_{\rm tot}$ is the normalization constant of the
original probability distribution $p_{i}$ and retains the value presented in
Eq. 17, while $\omega$ is the four-dimensional Gaussian distribution function
that reflects the uncertainties related to the $i$-th object, for all four
directly or indirectly measured quantities. That is:

$\displaystyle\omega$ $\displaystyle(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log N_{\mathrm{H},i},z_{i}|\log
M_{\rm BH},\log\lambda_{\rm E},\log N_{\mathrm{H}},z)$ (19)
$\displaystyle=\frac{1}{4\pi^{2}\sigma_{\log M_{\rm BH}}\leavevmode\nobreak\
\sigma_{\log\lambda_{\rm E}}\leavevmode\nobreak\ \sigma_{\log
N_{\mathrm{H}}}\leavevmode\nobreak\ \sigma_{z}}\times$
$\displaystyle\exp\left(-\frac{(\log M_{\mathrm{BH},i}-\log M_{\rm
BH})^{2}}{2\sigma_{\log M_{\rm
BH}}^{2}}-\frac{(\log\lambda_{\mathrm{E},i}-\log\lambda_{\rm
E})^{2}}{2\sigma_{\log\lambda_{\rm E}}^{2}}\right)$
$\displaystyle\times\exp\left(-\frac{(\log N_{\mathrm{H},i}-\log
N_{\mathrm{H}})^{2}}{2\sigma_{\log
N_{\mathrm{H}}}^{2}}-\frac{(z_{i}-z)^{2}}{2\sigma_{z}^{2}}\right)\,,$

where $\sigma_{\log M_{\rm BH}}$, $\sigma_{\log\lambda_{\rm E}}$,
$\sigma_{\log N_{\mathrm{H}}}$ and $\sigma_{z}$ correspond to the assumed
uncertainties on $\log M_{\rm BH}$, $\log\lambda_{\rm E}$, $\log N_{\rm H}$
and $z$ (respectively; see Table 1). Note that—at least for
$\log\leavevmode\nobreak\ M_{\rm BH}$ and $\log\leavevmode\nobreak\
\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 $\log N_{\rm H}$ are derived from extensive
spectral decomposition of multi-mission X-ray data (combining Swift-BAT and
ancillary spectra at $E<10$ keV), the convolution over
$\log\leavevmode\nobreak\ N_{\rm H}$ and $z$ (with
$\sigma_{\log\leavevmode\nobreak\ N_{\rm H}}=0.25$ and $\sigma_{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:

$\displaystyle\omega$ $\displaystyle(\log
M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i}|\log M_{\rm BH},\log\lambda_{\rm
E})$ (20) $\displaystyle=\frac{1}{2\pi\sigma_{\log M_{\rm
BH}}\sigma_{\log\lambda_{\rm E}}}\times$ $\displaystyle\exp\left(-\frac{(\log
M_{\mathrm{BH},i}-\log M_{\rm BH})^{2}}{2\sigma_{\log M_{\rm
BH}}^{2}}-\frac{(\log\lambda_{\mathrm{E},i}-\log\lambda_{\rm
E})^{2}}{2\sigma_{\log\lambda_{\rm E}}^{2}}\right)\,.$

We note that convolving over $z$ and/or $\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_{\log M_{\rm
BH}}=\sigma_{\log\lambda_{\rm E}}=0.3$ dex (see Table 1). The assumption of no
uncertainties (i.e., $\sigma_{\log M_{\rm BH}}=\sigma_{\log\lambda_{\rm
E}}=0$) is used only to test and demonstrate our analysis framework in
Appendix D. A higher value of $\sigma_{\log M_{\rm
BH}}=\sigma_{\log\lambda_{\rm E}}=0.5$ dex is required for Type 2 AGNs, to
take into account rotation and aperture effects in the optical spectroscopy
that is used to measure $\sigma_{\star}$ (see, e.g., Gültekin et al. 2009 and
Shankar et al. 2019). Additionally, we report results for a scenario where
$\sigma_{\log M_{\rm BH}}=0.3$ dex, along with a total uncertainty of 0.2 dex
in bolometric luminosity (${\sigma_{\log L,{\rm scatt}}}$), which reflects
both measurement and systematic uncertainties (dominated by the latter). As
$\log\lambda_{\rm E}$ is calculated using two observed quantities (luminosity
and mass), the uncertainty in luminosity measurement also contributes to the
total uncertainty in $\log\lambda_{\rm E}$. Therefore, the uncertainty in
$\log\lambda_{\rm E}$ is $\sigma_{\log\lambda_{\rm E}}\simeq 0.36$ (i.e. log
uncertainties in mass and luminosity added in quadrature). We note that this
approach to calculating $\sigma_{\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_{\rm BH}$, and
$\lambda_{\rm E}$ means that any systematic errors on $\log M_{\rm BH}$ and
$\log\lambda_{\rm E}$ would be anti-correlated, rather than independent. As
shown in § 4, the (conservative) various levels of uncertainty we assume do
not significantly affect our key results, which attests to the robustness of
our conclusions.

In Appendix E, we report the results of an even more complex scenario,
assuming luminosity dependent bolometric correction (from Duras et al. 2020),
and an even larger error on $\log\lambda_{\rm E}$. Our framework has also been
tested 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
either luminosity, $M_{\rm BH}$ or $\lambda_{\rm E}$ in the main analysis.
Thus, any difference between the XLF, BHMF, and/or ERDF derived for the
obscured and unobscured sub-samples (Type 1 and 2 AGNs) would occur
independent of our model assumptions. In Appendix E, where we present results
for a variable bolometric correction, some more complex assumptions are
introduced. For example, Type 1 and Type 2 AGN have different bolometric
corrections. This might artificially introduce difference in results between
the two populations, therefore we keep such assumptions to a minimum in our
main analysis.

####  3.4.4 Constraining $\Psi$

To determine the intrinsic bivariate distribution function, $\Psi(\log M_{\rm
BH},\log\lambda_{\rm E})$, we maximize the likelihood of observing our main
input sample. Expressing the $\log$-likelihood (Eq. 14) using Eq. 15:

$\displaystyle\ln\mathcal{L}=\sum_{i}^{N_{\rm obs}}\ln\left[N_{\rm
i,conv}(\log M_{\mathrm{BH},i},\log\lambda_{\mathrm{E},i},\log
N_{\mathrm{H},i},z_{i})\right]-\ln N_{\rm tot}$ (21)

Similar to the fitting procedure for the 1/$V_{\mathrm{max}}$ values, we use
the MCMC python package emcee to maximize $\ln\mathcal{L}$ (Foreman-Mackey et
al., 2013). The number of free parameters in these fits, which is six
($\gamma_{1}$, $\epsilon_{\lambda}$, $\lambda_{\rm E}^{*}$, $\alpha$, $\beta$
and $\mathrm{M}^{*}_{\rm 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_{\mathrm{max}}$ values, and 50 walkers are allowed to take 3,000 steps
(the chains usually converge within 2,500 steps).

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

$\Psi^{*}=a_{\rm rescale}\times\frac{N_{\rm obs}}{N_{\rm tot,best-
fit}(\Phi^{*}_{\rm init},\xi^{*}_{\rm init})}\,.$ (22)

Here $N_{\rm obs}$ corresponds to the total number of objects observed in our
sample, while $N_{\rm tot,best-fit}$ is determined by using Equation 17, the
best-fitting parameters for $\Psi$, and the initial normalizations
$\Phi^{*}_{\rm init}$ and $\xi^{*}_{\rm init}$. The parameter $a_{\rm
rescale}$ corresponds to an additional re-scaling factor which can be used to
correct for only partially available data (i.e., $z$ and/or $M_{\rm BH}$
measurements). As our sample is spectroscopically complete (excluding the non-
beamed, non-Galactic sources, along with others discussed in § 2.1.2), we
assume $a_{\rm rescale}=1$ for this work.

Finally, we use Eq. 13 to determine the bias-corrected intrinsic BHMF and
ERDF:

$\displaystyle\Phi_{\rm M}(\log M_{\rm BH})=\Psi^{*}\tilde{\Phi}(\log M_{\rm
BH},\Phi^{*}_{\rm init})$ (23) $\displaystyle\int_{\log\lambda_{\rm
E,min,s}}^{\log\lambda_{\rm E,max,s}}\tilde{\xi}(\log\lambda_{\rm
E},\xi^{*}_{\rm init})d\log\lambda_{\rm E}$ $\displaystyle\xi(\log\lambda_{\rm
E})=\Psi^{*}\tilde{\xi}(\log\lambda_{\rm E},\xi^{*}_{\rm init})$
$\displaystyle\int_{\log M_{\rm BH,min,s}}^{\log M_{\rm
BH,max,s}}\tilde{\Phi}(\log M_{\rm BH},\Phi^{*}_{\rm init})d\log M_{\rm BH}$

As $\Psi^{*}$ has the units of space density, namely h3 Mpc-3 dex-2, so do
$\Phi^{*}$ and $\xi^{*}$.

We 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 AGN and tested our method
using three different levels of uncertainty on the (mock) $\log M_{\rm BH}$
and $\log\lambda_{\rm E}$ values. These tests and mock catalogs are described
in detail in Appendix D.

##  4 Results - the Intrinsic Distributions Governing the BASS Sample

In this section we present our results for the XLF, the BHMF and the ERDF of
the BASS/DR2 sample, including the distributions corresponding to the subsets
of Type 1 and Type 2 AGN. We discuss some of the considerations made when
determining the XLF, present the XLF for AGN in different $\log N_{\rm H}$
bins (in addition to the Type 1 and Type 2 AGN XLFs), and compare to the
results of previous studies of AGN selected in the ultra-hard X-ray regime. We
then determine the bias-corrected BHMF and ERDF for BASS/DR2 AGN and again
compare our results to those of previous studies. We finally demonstrate that
the XLF that we derive from our fundamental bivariate, bias-corrected $M_{\rm
BH}$\- and $\lambda_{\rm E}$-dependent distribution can reproduce the XLF that
we measure directly from the observations.

We note that the 1/$V_{\mathrm{max}}$ based XLF reported in this work is a way
to understand the distribution of observed data with respect to AGN
luminosity, rather than an involved derivation of the intrinsic XLF (such as
the XLFs presented in Gilli et al. 2007, Ueda et al. 2014, Aird et al. 2015,
Buchner et al. 2015, and Ananna et al. 2019). Indeed, we report the
1/$V_{\mathrm{max}}$ based estimates (and fits) for all three distribution
functions (BHMF, ERDF and XLF) as a first-order approximation, and as a way to
make our results more directly and readily comparable with the observed
distributions reported in previous studies (e.g., Greene & Ho, 2009; Ajello et
al., 2012; Schulze & Wisotzki, 2010; Schulze et al., 2015).

###  4.1 The XLF of low-redshift AGN

To determine the AGN XLF of ultra-hard X-ray selected AGN, we use the
1/$V_{\mathrm{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 (Fig. 3, 4). Once we have determined
$\Phi_{\rm L}$, we fit the XLF with a double power-law (Eq. 9). The
1/$V_{\mathrm{max}}$ space densities are given in Appendix F (specifically
Table 11), and the best-fitting double power-law parameters of the various
sub-samples we consider are given in Table 2. The analysis steps are further
discussed in what follows.

In Figure 6 we present the XLF determined for BASS/DR2 AGN, 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\,{\text{keV}}}$), and column densities ($N_{\rm H}$) of all 672
non-beamed, high Galactic latitude BASS/DR2 AGN at $z\leq 0.3$; their XLF is
shown in the left panel with blue open squares and dotted line. Second, the
left panel of Fig. 6 shows the XLF of the 619 sources that meet our redshift
restrictions ($0.01\leq z\leq 0.3$), without any corrections due to
obscuration applied to their selection function. The only noticeable
difference between the XLF of all BASS/DR2 AGN and the redshift-restricted
sub-sample is in the low-$L$ end. The redshift-restricted sub-sample lacks the
15 lowest-$L$ AGN (i.e., in the $\log[L_{\rm X}/{\text{erg\,s}}^{-1}]<42.25$
bins), which are also some of the lowest redshift AGN in BASS/DR2 (i.e.,
$z<0.01$). At higher luminosities, the number of AGN in each luminosity bin
for the second sample is non-zero, and usually equal to the number of AGN in
the supersample. As shown in Figure 6 and Table 11, 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_{\mathrm{max}}$
(and therefore higher 1/$V_{\mathrm{max}}$) values in some cases.

We analyse the XLF of our samples further in the right panel of Figure 6. In
addition to the XLF of the redshift-restricted sub-sample, we show how this
XLF changes when obscuration corrections are applied to the selection
function. The obscuration correction results in slightly higher values of
1/$V_{\mathrm{max}}$ at all luminosities (as the $V_{\mathrm{max}}$ value of
the obscured object is smaller). The 1/$V_{\mathrm{max}}$ values and the
corresponding errors are given in Tables 8 and 11. The best-fitting parameters
for each XLF are shown in Table 2.

The correction due to obscuration is generally very small, which is expected
as the radiation in the 14$-$195 keV regime only starts to get attenuated
beyond $\log(N_{\rm H}/{\rm 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
$\log(N_{\rm H}/{\rm 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_{\mathrm{max}}$ calculations except
for the two cases shown in Fig. 6. In the Figure, we show the
1/$V_{\mathrm{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 AGN by Sazonov et al. (2007), Tueller et al. (2008), and Ajello et
al. (2012). We also show the binned 1/$V_{\mathrm{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 $\Gamma=1.8$. For the results by
Ajello et al. (2012), this implies $L_{14{-}195\,\rm{keV}}=2.31\times
L_{15-55\,\rm{keV}}$. To convert from $L_{20-40\,\rm{keV}}$ and
$L_{17{-}60\,\rm{keV}}$ to $L_{14{-}195\,\rm{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
AGN reaches down to hard X-ray luminosities that are $\approx 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\leq z\leq
0.3$) only goes down to $\log(L_{\rm X}/{\text{erg\,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_{\mathrm{max}}$ based $\Phi_{\rm L}$
values for our BASS/DR2 AGNs are in excellent agreement with previous results.

Table 2: Best-fitting XLF parametersa.

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

aafootnotetext: 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$.

Figure 6: Effect of the low redshift cut and obscuration on the BASS XLF. Left
panel: The BASS XLF determined without a low-$z$ cut (blue squares; $672$
sources), and the corresponding fit (dotted blue line). Green plus signs show
the XLF after excluding nearby sources with $z<0.01$ (leaving $659$ sources in
the sample), with the green dashed line fit to these $\Phi_{\rm L}$ values.
For comparison we show the XLFs by Sazonov et al. (2007), Tueller et al.
(2008), and Ajello et al. (2012). Right panel: The obscuring column density
$N_{\rm H}$ has a small but systematic effect on the XLF. Green plus signs and
green dashed line show the XLF computed by not taking obscuration into account
when computing the 1/$V_{\mathrm{max}}$ values, while the orange-filled
circles (shifted left by 0.05 dex for clarity) illustrate obscuration-
corrected $\Phi_{\rm L}$ values. The $N_{\rm H}$ correction consistently
increases the space densities in most luminosity bins. The turquoise crosses
(shifted right by 0.05 dex for clarity) and the turquoise solid line show the
XLF after applying all the cuts in Table 1 (586 sources). Figure 7: XLFs for
BASS AGN in different regimes of $N_{\rm H}$. Top panels: From left to right,
black points, black lines show the XLF for unobscured Compton-thin ($N_{\rm
H}<10^{22}\ \rm{cm}^{-2}$), obscured Compton-thin ($10^{22}\,\rm{cm}^{-2}\leq
N_{\rm H}<10^{24}\ \rm{cm}^{-2}$), and Compton-thick ($N_{\rm H}>10^{24}\
\rm{cm}^{-2}$) AGN. The XLF for the entire BASS sample is shown by blue filled
circles and blue solid lines. Also shown in the top right panel is the XLF of
BAT Compton-thick AGN from Akylas et al. (2016), which includes objects at
$z<0.01$, below the limit used in our analysis. Bottom panels: Ratios between
the $N_{\rm H}$ sub-sample XLFs and the XLF of all AGN in each $\log L_{\rm
X}$ bin (black data points), along with the associated Wilson score intervals
(black vertical lines). The numbers refer to the objects in each bin, with
upper and lower rows corresponding to the sub-sample and the full sample,
respectively). These XLFs and ratios should be treated as representing the
AGNs observed within BASS/DR2, in the corresponding absorption bins, rather
than intrinsic XLFs (which would also account for AGNs completely missed by
our survey; see main text).

###  4.2 The XLF of AGN with Various Levels of Absorption

The size of the BASS sample allows us to determine the XLF for various subsets
of AGN. In Figure 7, we split the AGN by $\log N_{\rm H}$, and calculate the
binned XLF separately for unabsorbed Compton-thin [$\log(N_{\rm H}/{\rm
cm}^{-2})<22$], absorbed Compton-thin [$22\leq\log(N_{\rm H}/{\rm
cm}^{-2})<24$] and Compton-thick [$\log(N_{\rm H}/{\rm cm}^{-2})\geq 24$]
sources. For comparison, the blue solid lines on 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 AGN reported by Akylas et al. (2016). They used 53 Compton-thick
AGN selected from the the same parent catalog as the one we use (Swift-BAT
70-month catalog), and determined a relatively low “break” luminosity in the
14$-$195 keV range, $\log(L^{*}_{\rm X}/{\text{erg\,s}}^{-1})\simeq 42.8$,
compared to our value of $43.59^{+0.38}_{-0.23}$. Note that we also have 53
Compton-thick objects in our parent sample, 11 of which are at $z<0.01$, one
falls above $z>0.3$ and another one which is a faint dual source. Our
individual 1/$V_{\mathrm{max}}$ $\Phi_{\rm L}$ values are generally consistent
with the fit by (Akylas et al., 2016, within errors) in all but the highest
luminosity bin, although at
$\log(L_{14-195\,{\text{keV}}}/{\text{erg\,s}}^{-1})>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 maximum likelihood estimator that considers individual
sources rather than a fit to 1/$V_{\mathrm{max}}$ (see also Loredo, 2004).

To further demonstrate the fractional densities of the $\log N_{\rm H}$ sub-
samples, the bottom panels of Figure 7 show the ratios between the space
densities of each sub-sample compared to the entire sample, as a function of
$L_{14-195\,{\text{keV}}}$. These ratios are computed using the
1/$V_{\mathrm{max}}$ $\Phi_{\rm L}$ values, and the errors on the ratio are
calculated using 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 Fig. 7, are
estimated using the 1/$V_{\mathrm{max}}$ method and thus 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 (which are accounted for in more sophisticated studies, e.g., Ueda et
al., 2014 and Ananna et al., 2019). Deriving such intrinsic $N_{\rm
H}$-dependent XLFs and ratios requires a much more involved approach than the
1/$V_{\mathrm{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 §5.

###  4.3 BHMF and ERDF of local AGN

Figure 8:  BH mass, Eddington ratio, and luminosity distribution functions of
BASS/DR2 Type 1 (left panels) and Type 2 (right panels) AGNs. For each type of
AGNs, combining the BHMF (top row) and ERDF (middle row) reproduces the
observed XLF (bottom row). Data points come from the 1/$V_{\mathrm{max}}$
analysis and dashed lines represent fits to these data points. The solid lines
and shaded areas show the final, bias-corrected intrinsic distribution
functions. Lines with different colors trace the intrinsic distributions
assuming various levels of uncertainty (see legend for details). Specifically,
the blue (left) and red (right) solid lines assume uncertainties of
$\sigma_{\log M_{\rm BH}}=\sigma_{\log\lambda_{\rm E}}=0.3$; the orange solid
lines assume $\sigma_{\log M_{\rm BH}}=0.3$ and ${\sigma_{\log L,{\rm
scatt}}}=0.2$ (or $\sigma_{\log\lambda_{\rm E}}=0.36$); and the green dotted
lines assume $\sigma_{\log M_{\rm BH}}=\sigma_{\log\lambda_{\rm E}}=0.5$ dex.
The intrinsic distributions of Type 1 AGNs (assuming $\sigma=0.3$) are shown
also in the right panels, to highlight differences between Type 1 and 2 AGNs.
Note that the 1/$V_{\mathrm{max}}$ points shown for each subset of AGNs were
calculated considering only the AGNs of the respective type (i.e., Type 1 or 2
AGNs). The 1/$V_{\mathrm{max}}$ values are reported in Tables 8, 9 and 10.
Figure 9:  Distribution functions for all BASS/DR2 AGNs that fall within the
criteria specified in Table 1 (586 sources). As in Fig. 8, the BHMFs, ERDFs
and reconstructed XLFs are shown in the top, middle and bottom panels,
respectively. We show the intrinsic distributions derived for all AGNs
assuming uncertainties of either $\sigma=0.3$ dex or $\sigma=0.5$ dex in $\log
M_{\rm BH}$ and $\log\lambda_{\rm E}$ (black solid and green dotted lines,
respectively). The orange solid lines show intrinsic functions, derived
assuming $\sigma_{\log M_{\rm BH}}=0.3$ and ${\sigma_{\log L,{\rm
scatt}}}=0.2$ (or $\sigma_{\log\lambda_{\rm E}}=0.36$). The black data points
in the bottom panel show the direct 1/$V_{\mathrm{max}}$ XLF estimates, and
the black dashed line is the fit to those points. Solid black lines show the
final, bias-corrected, intrinsic distribution functions for the complete AGN
sample. For comparison, we also show the distribution functions of Type 1 and
Type 2 AGNs (blue and red lines, respectively) from Figure 8.

To determine the BHMF and ERDF for local AGN we use the $M_{\rm BH}$ and
$\lambda_{\rm E}$ measurements of 586 ultra-hard X-ray selected BASS/DR2 AGN,
as described in §2.2. Since the $M_{\rm BH}$ (and thus $\lambda_{\rm E}$) of
Type 1 (broad line) and Type 2 (narrow line) BASS/DR2 AGN were 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 ERDF of the total BASS/DR2 AGN
sample.

To gain an initial guess for the two distribution functions we use the
1/$V_{\mathrm{max}}$ approach, assume functional forms, and fit the individual
$\Phi$ and $\xi$ values independently. For the BHMF and the ERDF we assume a
modified Schechter function and a double power-law, respectively (see Eqs. 10
and 11 in §3.3). To correct for sample truncation, i.e., the bias against low
mass and low Eddington-ratio AGN, we use the parametric maximum likelihood
approach outlined in §3.4. The 1/$V_{\mathrm{max}}$ values for both the BHMF
and ERDF are given in Appendix F.

Figure 8 shows the BHMF (top panels) and ERDF (center panels) for Type 1 and
Type 2 BASS AGN (left and right columns, respectively). The two bottom panels
in 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 §3.4.1. We are thus left with six free parameters: the break and
two slopes of the BHMF ($M^{*}_{\rm BH}$, $\alpha$, $\beta$), and the break
and two slopes of the ERDF ($\lambda_{\rm E}^{*}$, $\delta_{1}$,
$\delta_{2}$). We note that the high-$\lambda_{\rm E}$ slope, $\delta_{2}$, is
parametrized as $\delta_{1}+\epsilon_{\lambda}$ with $\epsilon_{\lambda}>0$.
We vary $\delta_{1}$ and $\epsilon_{\lambda}$ in the MCMC. We also recall that
for the bias correction we assume uncertainties of 0.3 dex and 0.5 dex on both
$\log M_{\rm BH}$ and $\log\lambda_{\rm E}$ (see Equations 18 and 20), as well
as an additional scenario with $\sigma_{\log M_{\rm BH}}=0.3$ and
${\sigma_{\log L,{\rm scatt}}}=0.3$ (i.e., $\sigma_{\log\lambda_{\rm
E}}=0.36$). The results from all these different $\sigma$ values reassuringly
converges on the same solution, as shown in the Figure. In Appendix E we
present BHMF/ERDF calculated 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.

The contour plots presenting the likelihoods resulting from our MCMC analysis
are shown in Figure 20 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 §3.4.2) in these tables. As shown in the Figure 21, our final
results for both torus geometries are consistent with each other.

Table 3: Sample truncation corrected BHMFa, and fit to /$V_{\mathrm{max}}$
values for the BHMF for Type 1, Type 2 AGN, and both samples togetherb.

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

aafootnotetext: We assume a modified Schechter function (see equation 10) for
the BHMF.

bbfootnotetext: We use a constant bolometric correction (Equation 2) to
compute these results.

Table 4: Sample truncation-corrected ERDFa, and fit to /$V_{\mathrm{max}}$
values for the ERDF for Type 1 AGN, Type 2 AGN and the full AGN sampleb.

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

aafootnotetext: We assume a double power-law shape for the ERDF (see equation
11).

bbfootnotetext: We use a constant bolometric correction (Equation 2) to
compute these results.

As discussed in §1 and shown in detail in Appendix D, the bolometric
luminosity function 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 BHMF and ERDF 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 Eqn. D1). 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_{\mathrm{max}}$ values would have led
us to underestimate the space density of low mass AGNs at the lowest mass bin
($\log[M_{\rm BH}/M_{\odot}]=6.65$) by $\gtrapprox 1$ dex for both Type 1 and
Type 2 sources. At low $M_{\rm BH}$ and/or $\lambda_{\rm E}$, the
1/$V_{\mathrm{max}}$ method underestimates the intrinsic space density of AGNs
due to the survey incompleteness. At high $M_{\rm BH}$ and/or $\lambda_{\rm
E}$, the 1/$V_{\mathrm{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 15 and 16), and Figure 17 of S15, as
measurement uncertainty ($\sigma_{\log M_{\rm BH}}$, $\sigma_{\log\lambda_{\rm
E}}$) increases, this overestimation increases at the high $M_{\rm BH}$ and/or
high $\lambda_{\rm E}$ end. This effect is a manifestation of the so-called
Eddington bias (Eddington, 1913): the uncertainty causes objects from lower
$M_{\rm BH}$ (or $\lambda_{\rm E}$) bins to scatter into higher $M_{\rm BH}$
(or $\lambda_{\rm E}$) bins, and vice versa. Intrinsically, there are always
fewer objects with higher $M_{\rm BH}$ (or $\lambda_{\rm E}$), therefore
scattering from the lower to higher bins causes significant overestimation of
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 AGN
is higher than that of Type 1 AGN at most masses. The bias-corrected ERDF of
Type 2 AGN has higher space densities at $\log\lambda_{\rm E}\ \leq-1.7$.
Beyond this point, the ERDF of Type 2 AGN drops off rapidly below the ERDF of
Type 1 AGN. This is because the break in ERDF of Type 2 AGN is at $\log$
$\lambda_{\rm E}$= $-1.657^{+0.087}_{-0.064}$, which is significantly below
the break of Type 1 AGN at $\log$ $\lambda_{\rm E}$=
$-1.152^{+0.089}_{-0.053}$. Figure 9 shows the Type 1 and Type 2 BHMF, ERDF
and XLF along with the overall sample results.

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

Figure 10:  Contour plot of the best-fit bivariate distribution functions (top
panels), BHMF marginalized over bins of $\log\lambda_{\rm E}$ (middle panels),
and ERDF marginalized over bins of $\log M_{\rm BH}$ (bottom panels). Type 1
AGN are shown on the left, Type 2 AGN on the right. In the top panels,
contours highlight levels of constant $\Psi$ values, with log-uniform spacing
(see color bar); black crosses mark the positions of $[\log(\lambda_{\rm
E}^{*})$, $\log(M^{*}_{\rm BH}/M_{\odot})]$; and black lines highlight the
flux limits at different redshifts. In the middle panels, we show the
integrals for $\log\lambda_{\rm E}=-2\pm 0.15$ (dotted pink lines),
$\log\lambda_{\rm E}=0\pm 0.15$ (dashed pink lines), and all $\log\lambda_{\rm
E}$, (solid pink lines). In the bottom panels, we show the integrals for
$\log(M_{\rm BH}/M_{\odot})=7\pm 0.15$ (dotted orange lines), $\log(M_{\rm
BH}/M_{\odot})=8.5\pm 0.15$ (dashed orange lines), and all $\log(M_{\rm
BH}/M_{\odot})$. Each $\log\lambda_{\rm E}$ and $\log(M_{\rm BH}/M_{\odot})$
bin over which a function was marginalized is shown in the top panels in the
corresponding color and line style (e.g., the purple dashed lines show the
$\log\lambda_{\rm E}$ range over which the BHMF is marginalized, and the
result of the marginalization is shown with dashed purple lines in the middle
panels). For the purpose of this plot we assume a constant flux limit of
$\log(F_{14-195\,{\text{keV}}}/{\text{erg\,s}}^{-1})=-11.1$ (which corresponds
to a sky coverage completeness of just 2–3%; see Fig. 3). The lines of
constant flux limit in top left panel are also lines of constant luminosity,
and integrating along these lines produces the local XLF. This figure
demonstrates that the bivariate distribution function fully captures the
statistical properties of AGN.  Figure 11: Comparison of the BASS BHMF and
ERDF for Type 1 AGN (left panels) and Type 2 AGN (right panels) to previous
studies. The blue/red solid lines show the intrinsic BASS BHMF and ERDF. Top
panels: The BHMF of Type 1 and Type 2 AGN. The green triangles shows the bias-
uncorrected 1/$V_{\mathrm{max}}$ results from Greene & Ho (2007, 2009),
calculated using a sample of local broad-line AGN from SDSS. The yellow
squares and dashdotted lines show the 1/$V_{\mathrm{max}}$ and incompleteness
corrected (but not uncertainty corrected) BHMF of Schulze & Wisotzki (2010),
calculated using local broad-line AGN from the Hamburg/ESO survey. The fully
bias-corrected Type 1 BHMF and predicted Type 2 BHMF from Schulze et al.
(2015), evaluated at $1<$ z $<2$, are shown with purple dashed lines. The
1/$V_{\mathrm{max}}$ points from Schulze et al. (2015), constrained using
VVDS, SDSS and zCOSMOS data sets, are also shown using open purple circles,
open purple triangles and open purple squares, respectively. For reference, we
also plot the total BHMF (including inactive black holes) from Shankar et al.
(2009) (turquoise shaded region). Bottom panels: The ERDF of Type 1 and Type 2
AGN. The data points and lines from Schulze & Wisotzki (2010) and Schulze et
al. (2015) follow the same scheme as in the top panels.

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 20). For all three AGN populations (Type 1, 2, and
overall), the parameters of the BHMF show significant correlation. The pairs
of parameters $(\alpha,\beta)$ and $(\alpha,\log M_{\rm BH}^{*})$ seem to be
anti-correlated, while the pair $(\beta,\log M_{\rm BH}^{*})$ is positively
correlated. These trends are rather expected: while fitting the same intrinsic
population, if we fix the break in BHMF at a higher mass, the slope at low
mass ($\alpha$) has to be shallower and the slope at high mass ($\beta$) has
to be steeper, to compensate for the higher mass break and produce a good fit.
For the ERDF, the slopes ($\delta_{1}$–$\epsilon_{\lambda}$) are negatively
correlated for all three samples. The break of the ERDF ($\log\lambda_{*}$) is
weakly positively correlated with $\delta_{1}$ and shows no significant
correlation with $\epsilon_{\lambda}$. The Figure also shows that even when
these functions converge to the same distribution for one parameter, it may
occupy distinctly different locations in the 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 2 AGN), 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 similarities (or lack thereof) in the shapes of
certain fitting functions shown in the Figure.

In Figure 11 we compare the BASS BHMF and ERDF for Type 1 AGN (left panels)
and Type 2 AGN (right panels) to previously published determinations of these
distributions. In the top panels we show the 1/$V_{\mathrm{max}}$ $\Phi$
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
$\Phi$ values determined by Greene & Ho (2007, 2009), based on a large sample
of SDSS AGNs with broad H$\alpha$ lines, which is not corrected for sample
truncation. The yellow data points and dash-dotted line illustrate the
1/$V_{\mathrm{max}}$ $\Phi$ values and incompleteness-corrected distribution,
respectively, determined by SW10 based on the quasar sample of the Hamburg/ESO
survey (HES; Wisotzki et al. 2000). Note that SW10 does not correct for
measurement uncertainties. The purple dashed lines in Fig. 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 AGN are predicted from the
corresponding distributions for Type 1 AGN, 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., 2018a; Shen et al., 2020). The
turquoise shaded region in the top panels of Fig. 11 also shows the total
local black hole mass function—that is, including both active and inactive
SMBHs determined by Shankar et al. (2009). The total BHMF naturally occupy a
much higher space density than our BASS AGN-only BHMF. The recent study by
Shankar et al. (2020) used a different set of $M_{\rm BH}-\sigma$ relations
(discussed in §2.2) to derive an updated total BHMF. However, as our analysis
relies on the more widely used $M_{\rm BH}-\sigma$ relation of Kormendy & Ho
(2013), we do not compare our results with the Shankar et al. (2020) BHMF.

###  4.4 Mass Independence of the ERDF Shape

Figure 12:  ERDF of Type 1 (top panel) and Type 2 (bottom panel) AGN divided
into two mass bins: $\log$ $M_{\rm BH}$$\leq 7.8$ (purple lines and points)
and $\log$ $M_{\rm BH}$$\geq 8.2$ (green lines and points). The data points
represent 1/$V_{\mathrm{max}}$ in each $\lambda_{\rm E}$ bin, the dashed lines
represent fit to the 1/$V_{\mathrm{max}}$ points, and the solid lines
represent the bias-corrected, intrinsic distribution functions. Dividing the
sample in this way shows that, for both Type 1 and Type 2 AGN, the shape of
the ERDF is independent of mass.

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 sub-samples into two
broad mass bins, $\log(M_{\rm BH}/M_{\odot})\leq 7.8$ and $\log(M_{\rm
BH}/M_{\odot})\geq 8.2$. The 0.4 dex wide gap in $\log M_{\rm BH}$ was imposed
to minimize “mixing” between the mass bins due to uncertainties on $M_{\rm
BH}$ estimation. We have then derived the 1/$V_{\mathrm{max}}$ measurements
(which are susceptible to selection biases, as discussed in detail in §3), the
associated functional fits of these measurements, and the bias-corrected,
intrinsic functional forms, for each of the four sub-samples (Type 1s and 2s,
low- 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- and high-mass
bins—that is, power-law exponents and the the location of the break—are in
excellent agreement (for each of the AGN sub-types), as can also be seen from
the best-fit parameters (Table 7). 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_{\rm BH}$. This analysis
indicates that the shape of the ERDF is indeed mass independent, for both
obscured and unobscured AGN.

##  5 Discussion

We have derived the XLF, BHMF and ERDF of a large, complete sample of ultra-
hard X-ray-selected, low-redshift AGNs from the BASS/DR2 dataset. Our analysis
included the direct 1/$V_{\mathrm{max}}$ approach as well as an elaborate
inference scheme that allowed us to recover the intrinsic distribution
functions, accounting for numerous generic, AGN-related, and survey-specific
potential biases. In both cases, we derived functional fits of the
distribution functions, considering the combined BASS/DR2 sample of (non-
beamed) AGNs, as well as sub-samples based on the AGN optical spectral classes
(Type 1s and 2s) and line-of-sight obscuration (parametrized by their $N_{\rm
H}$). We finally demonstrated that the intrinsic BHMF and ERDF can be combined
to reproduce the XLF, and that the assumption of mass-independence of the
shape of the ERDF is justified, at least for our sample.

Before moving to a higher-level discussion, we briefly list the main results
of our main analyses:

  * •

Our main result is the intrinsic BHMF and ERDF for Type 1 and Type 2 AGNs,
presented in Figure 8, and Tables 3 and 4. The method used to derive these
functions overcomes the limitations of the 1/$V_{\mathrm{max}}$ approach
(e.g., the Eddington bias; Eddington, 1913) and incompleteness at low
luminosities, masses, and/or accretion rates.

  * •

We show the source number counts for all the AGN samples drawn from the
Swift/BAT 70-month catalog in Figure 2, and provide the best-fit power-law
slopes and normalizations of the differential number counts for all samples
under study in §2.3.

  * •

We present the overall XLFs in Figure 6. The key observed quantities for these
XLFs (i.e., 1/$V_{\mathrm{max}}$ based space densities) are given in Tables 8
and 11, while their best-fit, broken-power-law parameters are given in Table
2.

  * •

We present the XLFs of unabsorbed objects, Compton-thin objects and Compton-
thick objects separately in Figure 7. The corresponding 1/$V_{\mathrm{max}}$
values and best-fits parameters are given in Tables 12 and 2, respectively.

  * •

The 1/$V_{\mathrm{max}}$ fits to the BHMF, ERDF, and XLF for Type 1, Type 2,
and combined AGN samples are given in Tables 9, 10, and 8, respectively. The
parameter fits to the observed, 1/$V_{\mathrm{max}}$ based measurements for
all three distributions, as well the intrinsic BHMF and ERDF are given in
Tables 3, 4, and 2 (again, respectively). The top (BHMF), middle (ERDF), and
bottom (XLF) panels in Figures 8 and 9 show the observed quantities and the
(intrinsic) best-fit distribution functions.

We next discuss several higher-level topics pertaining to the population of
accreting SMBHs in the local Universe, based on the results of our main
analysis. In particular, we (1) discuss the relation between accretion power
and (circumnuclear) obscuration; (2) compare the intrinsic BHMF and ERDF to
earlier studies and discuss possible ways to interpret their shapes; and (3)
discuss the AGN duty cycle.

###  5.1 AGN Obscuration and Demographics

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 ERDF). This is motivated
by the significant attenuation expected for highly obscured sources
[$\log(N_{\rm H}/{\rm cm}^{-2})>23$], even in the ultra-hard 14$-$195 keV band
(i.e., Fig. 4).

Second, the BASS sample allowed us to construct and explore the XLF for AGN of
several $N_{\rm H}$ regimes. The top panels of Fig. 7 show the XLF in each of
the three $N_{\rm H}$ bins, while the bottom panels show the fractions of
objects in each $N_{\rm H}$ bin relative to all AGN (in a luminosity-resolved
way; see details in §4.2). We find that (1) the fraction of unabsorbed AGN
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 suggested that unobscured AGN with higher
bolometric luminosities may have less (dusty) torus material, based on the
ratio of 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 AGN 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 AGN 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 often made the (pragmatic) assumption that Compton-thin
AGN 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 $\log(N_{\rm H}/{\rm
cm}^{-2})=24-26$ regime, thus predicting significant space densities of
$\log(N_{\rm H}/{\rm cm}^{-2})>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 $\log(N_{\rm H}/{\rm cm}^{-2})\ \geq 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 AGN would comprise $\approx$ 50% of all AGN 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 $\log(N_{\rm H}/{\rm cm}^{-2})<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, where the intrinsic $N_{\rm H}$ distribution shown