Constraining the Low-mass End of the Black Hole Mass Function and the Active Fraction of the Intermediate-mass Black Holes

We investigate the black hole mass function (BHMF) and the Eddington ratio distribution function (ERDF), focusing on the intermediate-mass black holes (IMBHs) with masses down to M • ∼ 104 M ⊙. Based on the active galactic nuclei (AGNs) with a detected broad Hα emission line, we construct a sample of 14,242 AGNs at redshift z < 0.35, including 243 IMBHs with M • < 106 M ⊙. By jointly modeling the BHMF and ERDF via the maximum posterior estimation, we find that the BHMF peaks at ∼106 M ⊙ and exhibits a relatively constant value of 10−4 Mpc−3 dex−1 at the low-mass end. By comparing the derived BHMF of type 1 AGNs with the galaxy mass function based on the updated black hole mass–host galaxy stellar mass relation, we derive the active fraction. We also determine the active fraction for all AGNs using the upper and lower limit of the type 1 fraction. The active fraction decreases from 15%–40% for massive galaxies (M ⋆ > 1010 M ⊙) to lower than ∼2% for dwarf galaxies with M ⋆ ∼ 108 M ⊙. These results suggest that the black hole occupation fraction is expected to be ∼50% for low-mass galaxies (M ⋆ ∼ 108.5–109 M ⊙) if the duty cycle is similar between IMBHs and supermassive black holes.


INTRODUCTION
The origin and evolution of supermassive black holes (SMBHs) are not well understood yet, despite their significance in galaxy evolution.The presence of SMBHs in almost all massive galaxies, as well as the correlation with the host galaxies, suggests the coevolution of black holes and galaxies (e.g., Rees 1984;Kormendy & Richstone 1995;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Kormendy & Ho 2013;Woo et al. 2013;Heckman & Best 2014).
At the low-mass end of the mass spectrum are the intermediate-mass black holes (IMBHs), which are typically defined as black holes with mass of 10 2 M ⊙ ≤ M • < 10 6 M ⊙ .IMBHs are considered to hold a key to understanding the origin of SMBHs, as various theoretical models predict different mass functions and occupation fractions of IMBHs in dwarf galaxies even in the local Universe (see Greene et al. 2020 for a review).For instance, the light-seed scenario argues that SMBHs originated from the remnants of Population III (Pop III) stars with a typical mass of ∼100M ⊙ (e.g., Haiman & Loeb 2001;Madau & Rees 2001;Heger et al. 2003;Volonteri et al. 2003).Since Pop III stars were prevalent in the early Universe, most of the galaxies in the current Universe should harbor SMBHs at their centers if SMBHs were formed by the light-seed channel.However, the presence of black holes heavier than 10 10 M ⊙ at z ∼ 6.3 (Wu et al. 2015) poses challenges for the light-seed scenario, as growth at super-Eddington rates would be necessary in a relatively short timescale.In contrast, the heavy-seed scenario suggests that the rapid collapse of halo gas in the early Universe formed IMBHs of M • ∼ 10 5 -10 6 M ⊙ (e.g., Loeb & Rasio 1994;Bromm & Loeb 2003;Begelman et al. 2006;Lodato & Natarajan 2006).In this case, SMBHs are not as prevalent in dwarf galaxies as the light-seed scenario predicts (e.g., Bellovary et al. 2019).However, this channel requires a mechanism that prevents the fragmentation of the gas (see Inayoshi et al. 2020b, and references therein).
Thus, it is important to investigate the black hole mass function (BHMF) and the occupation fraction to understand the origin of SMBHs.Previous constraints on the BHMF, however, have been insufficient for testing the black hole seed scenarios.One approach is to infer the BHMF from the mass function (or anything equivalent, such as the luminosity function (LF)) of galaxies.Shankar et al. (2004), for instance, derived local BHMF down to 10 6 M ⊙ by convolving the local galaxy LF with the M • -L sph relation (McLure & Dunlop 2002), implicitly assuming that all galaxies harbor a black hole at their center.Gallo & Sesana (2019) advanced this approach further by accounting for the occupation fraction in the galaxies observed in the local Universe (Miller et al. 2015).While this approach is useful, it relies on the choice of the scaling relation of the black hole mass against the host galaxy property, which is not well con-strained for the low-mass black holes (e.g., Kormendy & Ho 2013;Reines & Volonteri 2015;Woo et al. 2019).Furthermore, this approach cannot be employed to derive the occupation fraction, which is required for computing the BHMF.
The virial mass estimates of active galactic nuclei (AGNs) can provide constraints on the BHMF since the number of AGNs is the lower limit of the number of black holes.For instance, Kelly et al. (2009) demonstrated a Bayesian method to jointly model the BHMF and Eddington ratio distribution function (ERDF) as Gaussian mixtures.By applying the method to the 87 quasars in the Bright Quasar Survey (Schmidt & Green 1983), including 16 with reverberation mass estimates, they constrained the local (z < 0.5) BHMF down to M • ∼ 10 8 M ⊙ .Kelly & Shen (2013) applied the same method to the quasars in the Sloan Digital Sky Survey (SDSS; York et al. 2000) Data Release 7 (DR7), resulting in a similar BHMF for M • > 10 8 M ⊙ .Schulze & Wisotzki (2010) utilized AGNs in the Hamburg/ESO Survey (Wisotzki et al. 2000) and constrained the BHMF down to 10 6 M ⊙ based on a similar parametric approach assuming various functional forms of the BHMF and ERDF.Greene & Ho (2007) were one of the first to provide BHMF below 10 6 M ⊙ threshold by selecting AGNs in SDSS and measuring their M • based on the broad Hα emission.However, the BHMF and LF they derived with the 1/V max method (Schmidt 1968) were substantially smaller than those in other studies, including the aforementioned ones, which may suggest that the selection completeness they used was overestimated.Ultimately, these "BHMFs" only trace the mass function of active black holes; thus, they cannot be directly compared with the galaxy mass function to derive the occupation fraction without assumptions on the black hole activity, such as the duty cycle.
In this work, we present an improved estimate of the BHMF and ERDF based on the broad Hα-selected AGNs in SDSS DR7 by modeling the sample completeness more accurately than in previous studies.In Section 2, we define the sample and quantify the selection bias in our sample.Section 3 describes the mathematical details of our modeling.The results of our modeling are presented in Section 4, followed by a discussion in Section 5. Throughout this paper, we adopt a flat ΛCDM cosmology with H 0 = 72 km s −1 Mpc −1 and Ω m = 0.3.

THE SAMPLE
The sample of AGNs we used for measuring the BHMF is based on SDSS (York et al. 2000), which is a comprehensive optical survey of photometry and spectroscopy.Its seventh data release covers an area of over 10,000 deg 2 , with photometry of 357 million objects and 1.6 million spectra.The details of the survey were described by Abazajian et al. (2009).SDSS used the ugriz filter system (Fukugita et al. 1996) for its photometry.Based on their photometric measurements, they compiled two sets of spectroscopic targets covering 8032 deg 2 .The main galaxy sample (Strauss et al. 2002) consists of objects with r-band Petrosian (1976) magnitudes brighter than 17.7 mag (r Petro ≤ 17.7) and not a point source (r PSF − r model ≥ 0.3).Similarly, the quasar sample (Richards et al. 2002) consists of objects with i-band PSF magnitudes brighter than 19.1 mag (i PSF < 19.1).Furthermore, objects that would saturate the spectrograph were rejected (i fiber < 14.5 or r fiber < 15).Liu et al. (2019) carried out an exhaustive search of AGNs within the galaxy and quasar samples that exhibit broad Hα emission lines in their spectra.After removing the stellar continuum and modeled broad Hα lines, they selected AGNs based on the following criteria: (1) p-value smaller than 0.05 from the F -test for the broad component fitting, (2) broad Hα flux larger than 10 −16 erg s −1 cm −2 , (3) broad Hα signal-to-noise ratio larger than 5, (4) flux density at the peak of broad Hα larger than twice the continuum-subtracted rms, and (5) an FWHM of the broad Hα line profile larger than that of narrow lines.The resulting sample consists of 14,584 broad-line AGNs within z < 0.35, which corresponds to the redshift limit for Hα in the SDSS spectrograph, with their broad Hα luminosity spanning 10 38.5 -10 44.3 erg s −1 .
Based on the sample and measurements by Liu et al. (2019), we constructed our sample as follows.First, we eliminated any object with broad Hα flux larger than 7.5×10 −13 erg s −1 cm −2 as such objects on average would be brighter than the saturation limit of the survey (see Section 2.1).Then, we calculated broad Hα luminosity based on our cosmology and rejected objects with L Hα < 10 39 erg s −1 , as the number of objects below it was too small to produce meaningful statistics.We also estimated the black hole mass and the Eddington ratio from the luminosity and FWHM of broad Hα based on the single-epoch mass estimator and Hα bolometric correction of L bol = (190 ± 10) × L Hα provided by Cho et al. (2023) and removed any object with Eddington ratio smaller than 10 −3 .The resulting sample consists of 14,242 AGNs, with their mass, luminosity, and Eddington ratio spanning 4.5 < log 10 M • /M ⊙ < 10.5, 39 < log 10 L Hα /erg s −1 < 44.5, and −3 < log 10 L bol /L Edd < 0.5, respectively.Note that our sample includes 243 active IMBHs with M • < 10 6 M ⊙ .We present the dis-tribution of the redshift, mass, and Eddington ratio of our sample in Figure 1.
To correctly estimate the underlying population of AGNs, especially for the low-mass and low-luminosity end of the spectrum, the detection probability of each object in a sample, or the selection function, needs to be determined as accurately as possible.This probability then can be used to estimate the number of AGNs that are not included in the sample because of the selection procedure.Our sample was constructed based on two distinct samples (i.e., SDSS quasar and main galaxy samples), and the selection process utilized the measurements of various optical bands, which could be contaminated by the stellar emissions of the host galaxy.While we cover a wide range of mass and luminosity by combining two distinct AGN samples, it is challenging to correctly model the selection function of our sample, which was selected based on the broad Hα emission lines.In contrast, host galaxy correction is not required for AGN samples selected from X-ray surveys, although X-ray binaries can be a challenge to properly constructing an AGN sample, particularly at low luminosity (e.g., Gallo et al. 2008;Baumgartner et al. 2013).
We derive the selection function of our sample as follows.In Section 2.1, we construct the spectral energy distribution (SED) of pure AGNs without host galaxy contamination.Based on this, we derive the saturation limit of the AGNs in terms of Hα flux.Then, we model the distribution of the host galaxy fraction in Section 2.2.Finally, we derive the selection function considering the r-band flux from both the AGN and the host galaxy.

Spectral Energy Distribution of Quasars
The selection process of SDSS involves magnitudes of broadband photometry, which cannot separate the AGN from its host.To construct the selection probability based on the AGN luminosity only, the flux of broad Hα emission lines is a more robust proxy, as its distinctive feature can easily be separated from the host galaxy spectra.In this subsection, we investigate the correlation between r-or i-band flux and the broad Hα flux of AGNs, which we simply refer to as the SED of AGNs hereafter.
First, we compiled a list of objects with little to no host contamination.We defined the quasar subset of our sample, consisting of 7826 objects, to be the AGNs that were classified as quasar candidates but not as galaxy candidates based on the photometric selection criteria SDSS used (Richards et al. 2002;Strauss et al. 2002).However, a significant fraction of objects in this subset, typically fainter than 17.7 mag in the r band, appear to have Petrosian magnitudes brighter than PSF magnitudes by more than 0.3 mag, indicating that their flux is dominated by the extended host galaxy.We note that the r magnitude of 17.7 coincides with the magnitude limit of the galaxy sample of SDSS.Hence, we removed all objects with r Petro > 17.7 and r PSF − r model > 0.1.The latter condition is stricter than the corresponding selection criterion of the galaxy sample by Strauss et al. (2002) in order to remove all potentially extended sources.The final subset of pure AGNs consists of 1462 AGN-dominated objects.
To simplify the problem, we assumed that the shape of the SED does not depend on the luminosity of the AGN, i.e., the ratio of the r-or i-band flux to the broad Hα flux is the same for all AGNs.This is motivated by the Shakura & Sunyaev (1973) disk model, where the SED of the accretion disk follows a power law F ν ∼ ν 1/3 for wavelengths λ ≫ 10 3 Å for a wide range of accretion rates.Observationally, Richards et al. (2001) demonstrated that the r − i color of SDSS quasars is consistent across different magnitudes except for the redshift effect.Furthermore, the luminosity of the broad Hα emission line is proportional to the optical continuum luminosity over a wide range of luminosities (e.g., Yee 1980;Shen & Liu 2012;Jun et al. 2015;Cho et al. 2023).Given the narrow range of redshifts of AGNs in our sample, we expect the ratio of r-or i-band flux to Hα flux to be similar among the AGNs.Note that more distant AGNs may have different r −i colors owing to the 3000 Å bump (Grandi 1982) or the Lyman-α forest affecting the magnitudes, which should be taken into account when modeling the selection function of the optically selected AGNs at high z.
Under the assumption of the shared SED, we expected the PSF magnitudes of r or i bands to follow a linear function of the broad Hα magnitude, −2.5 log 10 F Hα , with slopes of unity.Indeed, the slopes obtained from ordinary least-square fits were consistent with unity within 1σ uncertainty in both bands, as presented in Figure 2. Fixing the slope to unity, the best-fit relations read r PSF = −2.5 log 10 F Hα 10 −6.25±0.01erg s −1 cm −2 (1) i PSF = −2.5 log 10 F Hα 10 −6.30±0.01erg s −1 cm −2 (2) We note that the relation did not change significantly even when we included all 7826 objects from the quasar subset.We used the SED to calculate the saturation limit of broad Hα flux.First, we found that for the subset of    Estimating the typical SED for the AGNs in the quasar subset.Black points were used for the fitting as described in Section 2.1, while gray points show the entirety of the objects that were flagged by the SDSS photometric pipeline as quasars but not as galaxies.The best-fit relations, fixing the slope to unity, are given as solid lines.The linear fits without fixing slope are shown with dashed lines for comparison.
pure AGNs r fiber − r PSF = 0.31 ± 0.01 [mag] (3) Then, saturation limits in both bands can be expressed as i PSF < 14.2 or r PSF < 14.7.Converting both magnitudes using Eqs. 1 and 2, we derived the saturation limit of Hα to be the smaller of two values, 7.5 × 10 −13 erg s −1 cm −2 .We removed any AGN in our sample above this threshold as described in Section 3.

The Flux Selection Function
While a majority of our sample is selected from the SDSS quasar sample, a subsample of 4308 objects are not in the quasar sample because the host galaxy flux dominates over AGN continuum in the broadband photometry.However, they are distinctively AGNs as demonstrated by a prominent broad Hα emission line.Since the quasar sample excluded any extended object whose color was similar to that of galaxies (Richards et al. 2002), these objects were included in the galaxy sample because they appeared as galaxies in r-band images (Strauss et al. 2002).Therefore, to assign a detection probability for these AGNs based on the AGN property, it is necessary to adopt a correlation between AGN and host galaxy properties.Specifically, for a given flux of broad Hα, we need to constrain the probability of being included in the galaxy sample by having r Petro ≤ 17.7, which we express as P (r Petro ≤ 17.7|F bHα ).
We first compared the r Petro against the r PSF for our entire sample.We found that the ratio of "extended flux" in r band to the PSF flux, R ≡ 10 −0.4(rPetro−rPSF) − 1, follows an exponential distribution, as shown in Figure 3.The best fit reads k = 0.44 ± 0.01.This, in turn, can be used to compute the probability of observing r Petro for an AGN given the value of r PSF : Based on this, we calculated our fiducial selection function at a given broad Hα flux, s(F bHα ), as follows: 3. Calculate the i PSF using Eq. 2.

s(F
To test the susceptibility of our fitting result to the selection function, we obtained an alternative fit with a mixture of two exponential distributions, which can replace p (R|k) in Eq. 6.We determined the best-fit parameters in Eq. 7 as k 1 = 0.44 ± 0.01, k 2 = 0.10 ± 0.02, and f = 0.01 ± 0.01, where the uncertainties were calculated based on 10000 bootstrap repetitions (see red lines in Figure 3).Using this best fit, we also calculated the alternative selection function as shown in Figure 4 (magenta line).We discuss the effect of the choice of the selection function on the BHMF in Section 4.1.

METHODS
In this section, we present the method for determining the BHMF and ERDF.We describe maximum posterior modeling in Section 3.1 and then we elaborate the 1/V max method in Section 3.2.

Maximum Posterior Modeling
We model the distribution function of the number density of black holes at a given mass and Eddington ratio, ϕ(M • , Λ), using a joint distribution of the black hole mass (M • ) and the Eddington ratio (Λ ≡ L bol /L Edd ).By marginalizing the other distribution function, the BHMF (ϕ(M • )) and the ERDF (ϕ(Λ)) can be respectively expressed as While ϕ is a function of (M • , Λ), it can be described with any other pair of variables if the pair can fully describe (M • , Λ) (i.e., the Jacobian determinant between the pairs is not zero).For example, the LF is expressed as For mathematical convenience, we use log-variables Similarly, the luminosity is expressed as where the bolometric correction L bol /L bHα = 190 is adopted from Cho et al. (2023).
. Probability distribution of the flux ratio between Petrosian and PSF magnitudes (left) and the magnitude differences ∆r ≡ rPSF − rPetro (right).The normalized histograms are denoted with black bars.The best fit for an exponential distribution (Eq.5) is shown with a blue solid line, while the alternative fit (Eq.7) is expressed with a red solid line, along with its two exponential components (dashed and dotted lines).
10 −15 10 −14 10 −13 10 −12 . Selection function of AGNs with a broad Hα line (see Section 2.2).The selection function used for modeling the BHMF is denoted with a solid black line by combining quasar and galaxy samples, while the selection function of the quasar sample is presented with a green line.The probability that an AGN is included in the galaxy sample (Eq.6) is expressed with a solid blue line, while the selection function constructed based on the alternative fit (Eq.7) is shown by a red solid line.
We assume that ϕ(m, λ) can be fully described by a set of parameters θ.Then, the probability of any black hole having m can be described as where is the comoving number density of the black holes.The probability for λ or l can be defined in a similar manner.
If we assume no redshift evolution of the number density between z min = 0 and z max = 0.35, the probability of a black hole having a redshift z should be proportional to the comoving volume at a given redshift: Given a set of observations q i = (l i , m i , z i ), the posterior p (θ|q i ) is written as (e.g., Kelly et al. 2009 where p(θ) is the prior distribution for θ, which is discussed in Section 3.1.2.The detection probability, p (I = 1|θ), is a probability of any black hole being observed, under the model given by θ, thus, where s(F bHα ) is the flux selection function of our sample described in Section 2.2.The luminosity integral can be expanded in terms of the integrals of m and λ, We note that p (I = 1|θ), which is adopted from Kelly et al. (2009), is analogous to the normalization described by Schulze et al. (2015, Eq. 5).
The last term in the posterior, n i=1 p (q i |θ), is the product of likelihoods for individually observed AGNs, and its calculation is described in Section 3.1.1.
Lastly, we constrain the comoving number density (ϕ • ) by matching the number of AGNs in our sample (N obs ) with the number of AGNs expected to be observed given the joint distribution function and the selection function, where Ω s is the survey area in units of sr.Therefore, we find where V survey = Ω s 4π zmax zmin dV dz dz is the survey volume.We note that Kelly et al. (2009) modeled the total number of AGNs in the survey volume using the negative binomial distribution by assuming a logarithmic prior for the total number of AGNs (p(log N ) ∼ const.).For a fixed set of parameters (θ), one can show that the expected value of this approach is identical to Eq. 20.

Uncertainties to Likelihood
As mentioned in Section 2, we used the single-epoch mass estimator based on the luminosity and the FWHM of the Hα line derived by Cho et al. (2023), which provides accurate mass estimates for the low-mass black holes.This is written as The scatter for this relation σ is given by σ = σ 2 SL + σ 2 f = √ 0.28 2 + 0.12 2 = 0.30 [dex], where σ SL = 0.28 is the intrinsic scatter of the size-luminosity relation and σ f = 0.12 is the uncertainty of the f -factor used in the mass estimator (Woo et al. 2015).We assume that the scatter follows a normal distribution in logarithmic units.We further consider the observational uncertainty of the luminosity as follows.Let l be the true luminosity of an AGN.We assume that the observed luminosity follows a lognormal distribution, l i ∼ N l, ε i 2 .Then, the true mass and Eddington ratio of the black hole, m and λ, can be written as where B = 0.61 ± 0.04 is the slope of the size-luminosity relation and t ∼ N 0, σ 2 is the intrinsic scatter.If we write the Gaussian kernel of width w as 2 , the likelihood can be written as In practice, the Gaussian convolutions of the likelihood were performed using 11th-order Gauss-Hermite quadratures (Abramowitz & Stegun 1972).The detection probability was calculated with fast Fourier transform with 3 × 10 −5 dex resolution, using the pyFFTW wrapper for the FFTW3 library (Frigo & Johnson 2005;Gomersall 2016).

Modeling the Distribution Functions
Given that the maximum posterior method is a parametric approach, we must express the model for ϕ(m, λ) with a set of parameters θ.For instance, Kelly et al. (2009) used a mixture of bivariate Gaussians to model the BHMF and ERDF.However, this approach requires the number of Gaussians to be large enough to accurately model the shape of the functions, especially when the dynamic range is wide.Instead, we chose our model to follow specific functions commonly used to describe luminosity and mass functions, to reduce the number of parameters to describe them.
We assume that the BHMF is independent of the Eddington ratio and the ERDF is independent of the black hole mass, similar to many studies in the literature (e.g., Schulze & Wisotzki 2010;Ananna et al. 2022).Then, we can write the joint distribution function (ϕ(m, λ)) as We use the following three functions to model the shapes of the BHMF and ERDF in this study.For consistency with the literature, we present them in terms of X, which can be either M • or Λ.For probabilities for log-scale variable x = log 10 X, the probability shape changes as A Schechter (1976) function (abbreviated as 'S') is motivated by the galaxy LF: This function is described by two parameters (apart from the comoving number density), X c and α.
Lastly, motivated by the double power law used in the literature (e.g., Schulze & Wisotzki 2010), we define a broken power law (abbreviated as 'bpl') as This function was chosen over the double power-law function because the broken power-law function does not suffer from degeneracy between the power indices γ 1 and γ 2 ; ϕ ∼ X γ1 always holds when X ≪ X c while ϕ ∼ X γ2 when X ≫ X c .Integrating these functions over X ∈ [0, ∞), or equivalently x ∈ (−∞, ∞), usually diverges, unless specific combinations of parameters are used.Thus, we normalize these functions by restricting our integration within the following intervals: m ∈ [4, 11] and λ ∈ [−3, 0.5].As a result, our derived BHMF only describes the distribution of (type 1) "active" black holes.
We imposed flat priors for the parameters describing the models, bounded by intervals to be considered physical.For the parameters describing characteristic scales (M c or Λ c ), the priors were chosen so that they are flat in the logarithmic scale.When the maximum posterior estimate was found near the edge of the interval, we relaxed the constraining interval by several dex and performed our analysis again.If it still failed to be constrained, we simply presented the last result as is and described the convergence.
In this study, we consider three pairs of two fitting functions for the BHMF and ERDF: (1) a combination of a broken power law for the BHMF and a modified Schechter function for the ERDF (hereafter bpl-mS), which is our best-fit model; (2) a combination of a modified Schechter function for the BHMF and a Schechter function for the ERDF (hereafter mS-S), which is motivated by the work of Schulze & Wisotzki (2010); and (3) a combination of broken power laws for both the BHMF and the ERDF (hereafter bpl-bpl), which is adopted from Ananna et al. (2022).We will discuss the differences in results using these models in Section 4.1.
3.2.1/V max and 1/V survey Methods While we mainly use maximum posterior modeling for constraining the BHMF and ERDF, we also use the traditional 1/V max method (Schmidt 1968) for comparison.Nonparametric estimates of the BHMF and ERDF can be obtained by simply counting the effective number of sources within a bin of mass, Eddington ratio, or luminosity and then dividing the counts by the survey volume V survey .The effective number of an object is inflated by a factor of V survey /V max,i , where V max,i represents the volume within which the object in question can be observed given the selection function.However, the 1/V max method fails to correct for the observational uncertainties, leading to Eddington bias at high luminosities (Eddington 1913).Furthermore, because the selection is based on the flux, not M • or λ, the derived BHMF or ERDF is biased against the low-mass or low Eddington ratio objects, which is referred to as "sample censorship" or "sample truncation" in the literature (e.g., Schulze & Wisotzki 2010;Ananna et al. 2022).Despite these limitations, this method is popularly used in the literature because of its simplicity.
The average function, ⟨ϕ x ⟩, within an interval x ∈ [x l , x u ] for an arbitrary variable x is given as and its uncertainty is We calculate V max,i for each object as where p (I = 1|x i , z) is the selection probability for an object having x i .For instance, if x i = l i , the selection probability becomes the flux selection function p (I = 1|l i , z) = s (F bHα (l i , z)), which we used for calculating the "flux-corrected" estimate of the LF.On the other hand, p (I = 1|x i , z) = 1 within comoving distances d c ∈ [d c,min,i , d c,max,i ] reduces V max,i to the volume of a truncated cone (e.g., Weigel et al. 2016) as To correct for the sample censorship, we follow the recipe described by Schulze & Wisotzki (2010).The selection probability for an AGN at a certain black hole mass or Eddington ratio is computed as Replacing the flux selection function in Eq. 32 with the mass selection probability or the Eddington ratio selection probability yields "mass-corrected" or "Eddington-ratio-corrected" estimates.Note that the mass-corrected estimate requires a priori knowledge of the ERDF, while the Eddington-ratio-corrected estimate requires a priori knowledge of the BHMF, in contrast to flux-corrected estimates (see discussion by Schulze & Wisotzki 2010).In addition, the masscorrected estimates still suffer from the Eddington bias.Therefore, we use 1/V max estimates for a consistency check.
Alternatively, we define another metric, namely, 1/V survey , which can be obtained by replacing V max,i in Eqs. 30 and 31 with V survey .Effectively, this yields a histogram of the sample normalized by the survey volume and the interval size.Note that this histogram is not corrected for any of the biases we discussed previously.Since 1/V survey estimates do not depend on any assumptions, we compare them with the prediction from our models to assess the goodness of the fit.For this, we define the "observed functions" as where The Gaussian kernel K σ in Eqs.36 and 37 represents the scatter (uncertainty) of the single-epoch mass estimator.We adopt σ = 0.30 [dex] unless a different scatter is assumed for the particular model.A convolution with this kernel simulates the Eddington bias induced by the uncertainty of the single-epoch mass, as discussed in Section 3.1.1.In the case of luminosities, the observational uncertainty, ε i , is present in 1/V survey estimates, and the average luminosity uncertainty of our sample is ⟨ε i 2 ⟩ ∼ 0.03 [dex] (< 0.09 [dex]), which is far smaller than the scatter of the single-epoch mass estimator.Thus, the observed functions are expected to closely replicate the 1/V survey estimates.We discuss the goodness of our models by comparing 1/V max and 1/V survey estimates in Section 4.1.

RESULTS
In this section, we present the BHMF and ERDF estimates by applying the methods presented in Section 3. We performed Markov Chain Monte Carlo (MCMC) simulations using the emcee library (Foreman-Mackey et al. 2013) to obtain the maximum posterior solution.We carried out 10,000 iterations of the MCMC simulation with an ensemble size of twice the number of parameters.The set of these parameters did not include the comoving number density (ϕ • ), which is computed after the posterior distribution is obtained.By inspecting the trace plots of each model, we concluded that the parameters always converged before the first 5000 iterations.From the latter half of iterations, we chose the maximum posterior parameters and their central 68% interval (equivalent to 1σ for normal distribution) as our best-fit estimates and uncertainties.The comoving number density (ϕ • ) was computed for each set of parameters in the posterior sample using Eq.20, and its best-fit values and uncertainties were selected.The results are summarized in Table 1.
We also calculated 1/V max estimates for fluxcorrected, mass-corrected, and Eddington-ratiocorrected cases, as well as 1/V survey estimates, as described in Section 3.2.We use the interval sizes of 0.5 dex for estimating the BHMF and LF and 0.25 dex for the ERDF.The minimum nonzero value for a 1/V survey estimate, corresponding to a singular object in a given interval, is 1.0 × 10 −9 Mpc −3 dex −1 for the BHMF and 2.1 × 10 −9 Mpc −3 dex −1 for the ERDF.

Comparison of BHMF and ERDF Estimates
Depending on the Functional Forms We compare the maximum posterior estimates, 1/V max estimates, and 1/V survey estimates in Figure 5.The maximum posterior estimates for the intrinsic BHMF and ERDF based on the bpl-mS and bpl-bpl models are generally similar, while the mS-S modelbased estimates show a different shape.This is because the exponential nature of the mS-S model cannot reproduce the shape of the observed BHMF at the high-mass end.Note that the uncertainties of ϕ(m) and ϕ(λ) are relatively small except for the high Eddington ratio end because of a small number in the bin.For example, there is only one object at λ > 0.25.
The 1/V max estimates are larger than the maximum posterior estimates, particularly for the ERDF.For example, the 1/V max estimate of ϕ(λ) is larger than the maximum posterior ϕ(λ) by more than 1 dex at λ > 0. This is because the correction factors for the sample censorship (i.e., Eqs.34 and 35) do not correct for the Eddington bias, while the maximum posterior estimates suffer no Eddington bias.Thus, the 1/V max estimates are suppressed near the maxima of 1/V survey estimates and enhanced elsewhere.This demonstrates the bias inherent to the 1/V max and 1/V survey estimators.Hereafter we only present 1/V survey estimates when comparing the observed functions without further discussing 1/V max estimates.
We also present the simulated observed functions for the BHMF and ERDF based on the maximum posterior estimates, after taking into account the selection function and the uncertainties of mass and Eddington ratio (dotted lines in Figure 5).Then, we compare the simulated observed functions with the 1/V survey estimates.
As described in Section 3.2, the simulated observed functions are expected to properly reproduce 1/V survey estimates.We find that the simulated observed function of the BHMF based on bpl-mS and bpl-bpl models is virtually consistent with the 1/V survey estimates.In contrast, the simulated observed function based on the mS-S model predicts a smaller number of high-mass AGNs and a larger number of low-mass AGNs than the other two model-based simulations, as well as the 1/V survey estimates.Given that we aim to constrain the BHMF at the low-mass end, we conclude that the mS-S model does not properly fit the sample.In the case of the ERDF, the simulated observed functions based on all three models are similar and generally consistent with the 1/V survey estimates.At the lowest λ, the discrepancy is pronounced, with the bpl-mS model predicting the lowest number of AGNs.Considering that the 1/V survey estimate is even smaller than the simulation of the other two model-based estimates, we conclude that the bpl-mS model reproduces the ERDF the best, although the difference among three models is not significant.For the rest of the paper, we present the bpl-mS modelbased maximum posterior estimates as the representative model.

Dependency of the Shape of the Selection Function and Mass Uncertainty
We investigate the effect of the different assumptions that we applied to model the selection function and the uncertainty of black hole mass on the estimates of the BHMF and ERDF.We assume an uncertainty of 0.3 dex the 0.3 dex mass uncertainty as "Standard bpl-mS".In addition, we consider two alternative cases.First, we adopt the alternative selection function described in Section 2.2, which effectively increases the selection probability at low flux (see the red line in Figure 4), and calculate the BHMF and ERDF estimates based on maximum posterior modeling.Second, we artificially increased the mass uncertainty, σ, from 0.3 to 0.6 dex (see Section 3.1.1)to consider unaccounted-for uncertainties in the mass estimate, and we carried out maximum posterior modeling, which is expressed as "2× Int.Scatter." In other words, we used σ = 0.6 dex for Eqs.36 and 37.
We compare the results based on the assumptions in Figure 6.Note that the choice of the selection function does not change the BHMF and ERDF significantly since the increase of the selection probability of low-flux AGNs is insignificant for changing the number density even at low-mass and low Eddington ratio bins.In contrast, the increased intrinsic scatter substantially modifies the shape of the ERDF, while its effect on the BHMF is not clearly detected.Thus, we expect that the mass uncertainty is not likely to be as large as 0.6 dex.
Compared with the observed function constrained by 1/V survey estimates, the simulated observed function of the BHMF and ERDF from maximum posterior modeling with an assumption of 0.6 dex mass uncertainty shows a very different shape for both the BHMF and the ERDF, implying that the mass uncertainty is much smaller than 0.6 dex.Nevertheless, the change in the BHMF ϕ(m) caused by different assumptions is smaller than 0.5 dex over all mass scales.Thus, we conclude that our best fit for the BHMF is stable against the choice of different assumptions.

Comparison with the Literature
To compare our results with the literature, we compiled the reported density functions from Greene & Ho (2007, BHMF, LF(bHα) (1/V max )), Shankar et al. (2009, BHMF(all), LF(bol)), Schulze & Wisotzki (2010, BHMF, ERDF (mS, S)), Kelly & Shen (2013, BHMF, ERDF), and Ananna et al. (2022, BHMF, ERDF, LF(Xray) (type 1)).In the case of the LF, we collected not only the broad Hα LF but also LFs for bolometric luminosity (Shankar et al. 2009) and X-ray (Ananna et al. 2022).We adopted an X-ray bolometric correction of L bol /L 14-195 keV = 7.4 (Ananna et al. 2022) and an Hα bolometric correction of L bol /L bHα = 190 (Cho et al. 2023) to convert the given LF to an Hα LF.We also note that the BHMF by Shankar et al. (2009) is based on the galaxy LF and the black hole mass scaling relations.Thus, this is a total BHMF rather than a type 1 active BHMF, which we present along with the rest of the literature.Uncertainties for the functions were collected whenever available.When they were provided as parameter uncertainties, we determined the 1σ uncertainties based on the 1000 Monte Carlo realizations.Lastly, we multiplied the appropriate factor to compensate for the differences in the Hubble constant from the one we used (h = 0.72).
In Figure 7, we present our best-fit BHMF, ERDF, and broad Hα LF along with the results collected from the literature.Regardless of the method they used, we plotted the literature functions as solid lines.We showed the functions in the intervals in which they were originally presented.If the functions were presented as values, we simply showed them all.Conversely, if they were provided as functions and parameters, we considered those functions valid only within the x-axis intervals of the figures plotting the functions in their respective references.
The BHMF we found follows a broken power law with a peak near the mass M • ∼ 10 6 M ⊙ .For the intermediate-mass regime, we find the BHMF to be relatively constant in logarithmic scale near ∼ 10 −4 Mpc −3 dex −1 .Compared with the literature, we found a striking consistency with the functions reported by Schulze & Wisotzki (2010) and Ananna et al. (2022) for masses 10 6.5 M ⊙ < M • < 10 9 M ⊙ ; for higher mass, they predicted much smaller values of the BHMF, which may be because the Schechter-like models they used in their fits are asymptotic to exponential decay at the high-mass end.Kelly & Shen (2013) predicted a BHMF with a similar value to our fit at higher masses (> 10 8 M ⊙ ) with a steeper slope, possibly due to their choice of priors on the Gaussian mixture modeling.On the contrary, we predict a significantly larger number of black holes than what Greene & Ho (2007) found using the broad Hα from SDSS DR4 spectra.While there could be many explanations involving the spectroscopic completeness of DR4 itself, the choice of the decomposition method used by the authors, etc., the most likely explanation would be that the completeness of the sample was overestimated.This is further supported by the observation that the other studies using similar sample and selection criteria, notably Kelly & Shen (2013), who used SDSS DR7, predicted a much larger BHMF than what Greene & Ho (2007) found, similar to ours.Lastly, since we constructed our BHMF based on the AGNs, our BHMF is ∼0.1%-10% of the total BHMF estimated by Shankar et al. (2009) or Gallo & Sesana (2019).We discuss the active fraction in detail in Section 5.1.
Our ERDF is less consistent with the literature, although this is to be expected.Since our BHMF spans a larger dynamic range than that in the literature, the total number density involved in our sample should be much higher.Consequently, the ERDF should also be much higher in proportion to the total number density.Focusing on the shape, our ERDF is more concentrated toward an Eddington ratio of a few percent, while the literature functions predict a monotonically decreasing ERDF.As demonstrated in Section 4.1, assumptions on the mass and/or luminosity uncertainty affect the resulting ERDF significantly.Nevertheless, should the ERDF be smaller for Λ <0.1-1% as with our best-fit model, the physical mechanism behind it could be the change in the accretion state (e.g., Inayoshi et al. 2020a).Another possible explanation could be that there is a minimum required accretion rate for a broad-line region (BLR) to exist (e.g., Nicastro 2000).Then, the discrepancy in the shape between our ERDF and the X-ray-based ERDF by Ananna et al. (2022) can be explained as the BLR being absent for the low-λ AGNs.Lastly, our broad Hα LF is consistent with the X-ray and bolometric LF by Shankar et al. (2009) and Ananna et al. (2022).We observe the difference between our LF and the Hα LF by Greene & Ho (2007), and this is most likely to be of the same origin as the differences in the BHMF.

Active Fraction
The black hole occupation fraction of low-mass galaxies provides valuable information on the origin of SMBHs (e.g., Greene et al. 2020).However, it is challenging to observationally constrain the occupation fraction because of the lack of a complete sample of IMBHs in dwarf galaxies.The active fraction is defined as a fraction of galaxies hosting actively mass-accreting black holes, providing a lower limit of the occupation fraction.
Compared to the occupation fraction, the active fraction is easily obtained through AGN surveys.We measure the BHMF of AGNs across a wide dynamic range, including the intermediate-mass regime (M • < 10 6 M ⊙ ), with an overall precision of much less than 1 dex.Thus, our sample is ideal for constraining the active fraction at the low-mass end.In this section, we present the active fractions of galaxies in the local volume.
First, we investigate the correlation between the black hole mass and the total host galaxy mass.A number of previous studies presented the black hole mass correlation with bulge properties, focusing on relatively massive galaxies (e.g., Kormendy & Ho 2013;Woo et al. 2013).However, the bulge properties cannot trace the total stellar mass reliably, due to the diverse range of the bulge-to-total stellar mass ratio (e.g., Khochfar et al. 2011).
For low-mass galaxies, Reines & Volonteri (2015, see also Reines & Volonteri 2019) reported a broad correlation between the black hole mass and the total stellar mass based on 271 local AGNs, including a number of IMBHs.In their study, the total stellar mass was estimated based on the color-dependent mass-to-light ratio from Zibetti et al. (2009).We use this sample for our analysis after updating single-epoch black hole masses based on the new estimator Eq. 21).Note that the previous single-epoch mass of IMBHs could be overestimated by 0.2-0.5 dex (see the discussion by Cho et al. 2023).Thus, we redetermine the black hole mass of 255 AGNs in their sample.In addition, we adopt the reverberation mass of NGC 4395 from Cho et al. (2021).with an rms scatter of 0.61 dex (see Figure 8).Second, we estimate the active fraction as follows.We translate the BHMF to the corresponding mass function as a function of stellar mass using Eq.39 as Then, the translated MF is convolved with a Gaussian kernel with a width of 0.61 dex, which is the rms scatter of the M • -M ⋆,total relation, to obtain the mass function of active galaxies.Dividing this by a galaxy mass function yields the active fraction.We adopt the galaxy mass function from Weigel et al. (2016).
In addition, we consider two factors that may affect the active fraction.First, we assume a type 1 fraction (f type1 ), which is the fraction of type 1 AGNs among all AGNs.Since we use type 1 AGNs to derive the BHMF, the active fraction based on this BHMF only counts type 1 AGNs.Thus, we overcome this limitation by dividing the active fraction of type 1 AGNs by the type 1 fraction, obtaining the total active fraction.Previous studies suggest that the type 1 fraction increases with increasing black hole mass.For instance, Lu et al. (2010) reported that the type 1 fraction is ∼20% for M • < 10 8 M ⊙ and ∼30% for M • > 10 8 M ⊙ .Similarly, Oh et al. (2015) found that the type 1 fraction increases from ∼20% at M • ∼ 10 6 M ⊙ to ∼60% at M • ∼ 10 9 M ⊙ .In contrast, Moran et al. (2014) found only 2 type 1 AGNs out of 28 AGNs in dwarf galaxies, corresponding to f type1 ∼7%.Since we expect the type 1 fraction to be a function of the host galaxy mass or black hole mass, we assume that the upper and lower limits of the type 1 fraction are 20% and 7%, respectively.We find that the active fraction of type 1 AGNs is larger than ∼3% for massive galaxies with M ⋆ > 10 10 M ⊙ as shown in Figure 9.We also present the active fraction of all AGNs (blue shaded region), using the upper and lower limits of the type 1 fraction (i.e., 20% and 7%, respectively), constraining the active fraction of all AGNs to be between 40% and 15% for massive galaxies.In contrast, the active fraction decreases for lower-mass galaxies, particularly at M ⋆ < 10 9.5 M ⊙ .At the low end of stellar mass (M ⋆ = 10 8 M ⊙ ), the active fraction is ∼2%, even with a lower limit of the type 1 fraction of 7%.
Our results are broadly consistent with the previous study based on deep X-ray observations by Gallo et al. (2010), who reported the active fraction to be 24%-34% on average, varying from 0.7% to 87% depending on the host galaxy mass.Note that they used very low Eddington ratio AGNs (λ < −3) for deriving the active fraction, while the lower limit of the Eddington ratio is λ ≥ −3 in our sample.Gallo et al. (2010) described that the decreasing trend of the active fraction with decreasing host galaxy mass was caused by the fact that the fixed X-ray luminosity limit of their survey corresponds to a higher Eddington ratio for lower-mass black holes in lower-mass galaxies.Nevertheless, their average active fraction of 24%-34% is a rough approximation for the active fraction.If we apply the lower and upper limits of the type 1 fraction (7% and 20%) to their active fraction, then we obtain the active fraction of type 1 AGNs as 1.68%-6.8%.A more recent study by Ananna et al. (2022) based on the X-ray luminosity reported the active fraction of 10-16% for AGNs with λ > −3 AGNs with 10 6.5 ≤ M • /M ⊙ ≤ 10 10.5 .In the case of type 1 AGNs with λ > −3, they found the active fraction to be ∼4% (see Figure 13 in Ananna et al. 2022).Thus, our result of ∼3% active fraction based on type 1 AGNs is consistent with both studies.Pacucci et al. (2021) constrained the active fraction of galaxies with M ⋆ < 10 10 M ⊙ , using AGNs with λ > −4, by adopting simple assumptions on the accretion efficiency.They predicted 2%-20% of active fraction, with an increasing trend with host galaxy mass, which is broadly consistent with our result.However, the increasing slope of their active fraction is shallower than ours but consistent with the range of active fraction that we constrained with the upper and lower limits of the type 1 fraction.Their shallower slope can be explained with a varying type 1 fraction between dwarf galaxies (Moran et al. 2014) and the more massive galaxies (Lu et al. 2010;Oh et al. 2015).If the BLR is less common in active IMBHs, then the type 1 fraction can be smaller for IMBHs.

Black Hole Occupation Fraction
We derive the black hole occupation fraction from the active fraction by assuming a duty cycle, which is the fraction of active black holes among all black holes.We adopt a constant duty cycle of 13%, which is an average of 10% and 16% constrained based on a sample of AGNs selected with X-ray by (Ananna et al. 2022).As done for the active fraction, we also use two different fractions,  7% (Moran et al. 2014) and 20% for type 1 AGNs (Lu et al. 2010;Oh et al. 2015).
We present the occupation fraction in comparison with the previous results in the literature in Figure 10.The occupation fraction is close to 1 for with M ⋆ ≳ 10 9.5 M ⊙ .This is consistent with the expectation that all massive galaxies host SMBHs (e.g., Rees 1984;Kormendy & Richstone 1995).The occupation fraction departs from unity for galaxies with M ⋆ ≲ 10 9.5 M ⊙ , becomes 0.5 for galaxies with M ⋆ ∼ 10 8.5 -10 9 M ⊙ , and reduces to ∼0.1 for dwarf galaxies with M ⋆ ∼ 10 8 M ⊙ .Our result is broadly consistent with the occupation fraction reported by Miller et al. (2015) based on the X-ray detection, supporting the validity of the assumptions involved in the calculation of the occupation fraction.

Comparison with the Seed Scenarios
We compare the occupation fraction with the theoretical predictions based on heavy-and light-seed models in the literature in Figure 10.For the heavy-seed scenario, we adopt the occupation fraction by Bellovary et al. (2019), which was calculated based on a cosmological zoom-in simulation.For the light-seed scenario, we use a couple of occupation fractions by Ricarte & Natarajan (2018, see Figure 5 of Greene et al. 2020), which were calculated with two different assumptions on the accretion rate.
The occupation fractions of dwarf galaxies (M ⋆ < 10 9 M ⊙ ), predicted from the heavy-seed model (Bellovary et al. 2019) or from the optimistic lightseed model (Ricarte & Natarajan 2018, MS), are larger than the occupation fraction we found.On the contrary, the pessimistic light-seed model (Ricarte & Natarajan 2018, PL) is far smaller than our occupation fraction.
Given the large uncertainty of the occupation fraction, it is difficult to discriminate between different seed scenarios.The uncertainty arises largely from the assumptions on the activity of black holes, i.e., the duty cycle and the type 1 fraction, that we used to construct the occupation fraction.Therefore, it is necessary to quantify the activity of IMBHs in detail, to better constrain the occupation fraction.

Implications for Future Surveys
The BHMF derived in Section 4 can be used for estimating the number of AGNs expected in a specific survey, based on Eq. 19, along with the sensitivity function of the survey.In this section, we discuss the number of IMBHs (10 4 M ⊙ < M • < 10 6 M ⊙ ) expected from future surveys.
For instance, the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) covers an area of 18,000 deg 2 in the sky.Its limiting magnitude in r band is 24.7 in a single visit and 27.5 in the final co-added survey, while objects with r < 16 would saturate.Converting these magnitudes with Eq. 1 and constructing a simple top-hat-like selection function, we expect that ∼500,000 active IMBHs within z < 0.35 would be detected in a single visit and ∼2 million active IMBHs in the final survey.
One of the promising future all-sky surveys is SPHEREx (Doré et al. 2014), which will cover all sky with a depth of 18.5 AB mag and a 100 deg 2 area with a depth of ∼22 AB-mag over the wavelength range of 0.75-5µm (Doré et al. 2018).While Hα emissions of AGNs at z ≥ 0.143 can be detected in the survey, AGNs at lower redshift can be detected via alternative signatures, such as Paα.Considering that the flux of the Paα line is dimmer than that of Hα by a factor of 10 in type 1 AGNs (Netzer 1990), we can calculate the Paα flux for given λ.On the other hand, AB magnitude can be derived by averaging the Paα flux over the spectral resolving power R = 41 of SPHEREx at low wavelengths.Comparing these, we derive the 5σ limiting flux of Paα to be 5.6 × 10 −14 erg s −1 cm −2 for all sky and 2.2 × 10 −15 erg s −1 cm −2 for a 100 deg 2 area.From this, we expect the number of active IMBHs detected with Paα within a z < 0.143 volume to be ∼25 for a 100 deg 2 area and ∼100 for all sky, without considering Galactic obstruction.

Summary
In this paper, we derived the mass function and ERDF of type 1 AGNs with masses down to 10 4 M ⊙ , using the local AGNs in SDSS DR7.By constructing the selection function based on the detection probability of the broad Hα line fluxes and applying maximum posterior modeling, we obtained the intrinsic and simulated density functions.We summarize the main results as follows: 1.The BHMF peaks at M • ≃ 10 6 M ⊙ while it is flat or slightly decreasing for IMBHs, as similarly reported in the literature.
2. The ERDF peaks near 1%-10% and follows a Schechter-like or a broken power-law function, exponentially decreasing as the Eddington ratio exceeds 10%.This trend is in contrast to the monotonically decreasing function reported in the literature, many of which used X-ray data to constrain the ERDF.
3. The type 1 active fraction is ∼3% for galaxies with M ⋆ > 10 10 M ⊙ .Adopting that the type 1 AGNs constitute 7%-20% of all AGNs with M • < 10 8 M ⊙ , we constrain the active fraction as 40%-15%.In the case of dwarf galaxies with M ⋆ ∼ 10 8 M ⊙ , the active fraction is ∼2% even with a lower limit of the type 1 fraction of 7%.
4. We constrain the black hole occupation fraction using the derived active fraction.For dwarf galaxies with M ⋆ ∼ 10 8.5 -10 9 M ⊙ , we find that the occupation fraction is 50%, which is consistent with the constraints based on the X-ray detections in the literature.
We thank the anonymous referee for constructive comments, which helped improve the quality of the manuscript substantially.H.C. would like to thank Jubee Sohn, Ena Choi, Shu Wang, Dongkok Kim, and Mankeun Jeong for the helpful discussions.This work has been supported by the Basic Science Research Program through the National Research Foundation of Korean Government (2021R1A2C3008486).J.-H.W would like to thank his collaborators for various inputs, despite the serious budget cut suddenly imposed by the government.

Figure 1 .
Figure 1.Distributions of the mass (left) and the Eddington ratio (right) as a function of redshift, along with histograms of the redshift, mass, and Eddington ratio.

Figure 2 .
Figure2.Estimating the typical SED for the AGNs in the quasar subset.Black points were used for the fitting as described in Section 2.1, while gray points show the entirety of the objects that were flagged by the SDSS photometric pipeline as quasars but not as galaxies.The best-fit relations, fixing the slope to unity, are given as solid lines.The linear fits without fixing slope are shown with dashed lines for comparison.

Figure 7 .
Figure7.The best fit of the BHMF (left), ERDF (middle), and broad Hα LF (right), compared with the literature.The maximum posterior prediction of the bpl-mS model its 1σ uncertainties are denoted by the black solid line and the shaded region.Filled circles denote the mass-corrected estimates obtained using the Vmax method, as described in Section 3.2.Details of the individual functions from the literature are described in Section 4.3.

Figure 8 .
Figure8.The M•-M⋆ ,total relation.The solid black line is our best fit, while the solid yellow line represents the relation presented originally byReines & Volonteri (2015).The shaded region shows the 1σ uncertainty of the best-fit relation, and the dotted lines mark the rms scatter around the relation (0.61 dex).

Figure 9 .
Figure9.The active fraction of galaxies based on our BHMF.The x-axis is the total stellar mass M⋆ or corresponding black hole mass M• based on Eq. 39.The solid blue line is based on the type 1 BHMF, with its 1σ uncertainty denoted with dashed lines.The total active fraction, assuming type 1 fractions between 0.07 and 0.2, is shown as a hatched region.For comparison, the active fractions fromGallo et al. (2010, beige shading)  andPacucci et al. (2021,  magenta line)  are presented.

Table 1 .
The Best-fit Parameters of the BHMF Model Figure10.The black hole occupation fraction of galaxies, obtained by assuming type 1 fractions between 7% and 20% and a duty cycle of 13%, is shown as a hatched region.The x-axis is identical to that in Figure9.The occupation fractions from the literature, as discussed in Section 5.2 and Section 5.3, are presented by colored lines and a shaded region.