CONSTRAINTS ON BLACK HOLE GROWTH, QUASAR LIFETIMES, AND EDDINGTON RATIO DISTRIBUTIONS FROM THE SDSS BROAD-LINE QUASAR BLACK HOLE MASS FUNCTION

We model the BLQSO BHMF as a mixture of K log-normal distributions:

$$\begin{eqnarray} \phi (\log M_{\rm BH},\log z) &=& N \left(\frac{dV}{dz}\right)^{-1} \sum _{k=1}^{K} \frac{\pi _k}{2 \pi |\Sigma _k|^{1/2}} \nonumber\\ &&\times\exp \left[ -\frac{1}{2} ({\bf y} - \mu _k)^T \Sigma _k^{-1} ({\bf y} - \mu _k) \right],\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where ∑^K_{k = 1}π_k = 1. In this work, we choose K = 4, based on our previous experience in working with simulated data sets (KVF09), and because we did not notice a significant difference in the results obtained using larger values of K. Here, N is the total number of BLQSOs in the universe that could be seen by an observer on Earth at the time of the survey, y = (log M_BH, log z), μ_k is the two-element mean vector for the kth Gaussian functions, Σ_k is the 2 × 2 covariance matrix for the kth Gaussian function, and x^T denotes the transpose of x. In addition, we denote π = (π₁, ..., π_K), μ = (μ₁, ..., μ_K), and Σ = (Σ₁, ..., Σ_K). The variance in log M_BH for Gaussian function k is σ²_m,k = Σ_11,k, the variance in log z for Gaussian function k is σ²_z,k = Σ_22,k, and the covariance between log M_BH and log z for Gaussian function k is σ_mz,k = Σ_12,k. The parameters for the mass function are N, π, μ, and Σ.

The distribution of luminosity density at a given black hole mass and wavelength λ is assumed to also follow a mixture of J log-normal distributions:

$${\selectfont\fontsize{8}{9}{\begin{eqnarray} &&\lefteqn{p(\log L_{\lambda }|M_{\rm BH})}\nonumber\\ &&= \sum _{j=1}^{J} \frac{\gamma _j}{\sqrt{2 \pi \sigma _{l,j}^2}}\exp \left[ -\frac{1}{2} \left(\frac{\log \lambda L_{\lambda } - \alpha _{0,j} - \alpha _{m,j} (\log M_{\rm BH} - 9)}{\sigma _{l,j}} \right)^2 \right].\nonumber\\ \end{eqnarray}}} \tag{ 2 }$$

This represents an extension of the model of KVF09, which only used J = 1 log-normal distribution. We made this extension to incorporate additional flexibility in p(L_λ|M_BH), ensuring that the luminosity distribution, and therefore the black hole mass incompleteness correction, is robust to the particular parametric form. For each log-normal distribution, the parameters for the distribution of L_λ at a given M_BH are γ_j, α_0,j, α_m,j, and σ_l,j. We used J = 3 log-normal distributions, as we did not notice a significant change in the estimated values of p(L_λ|M_BH) when using J ⩾ 3. We further assess the robustness of our assumed form for p(L_λ|M_BH) in Section 3.3, and show that our results should not be significantly altered if the true form of p(L_λ|M_BH) is a power law, as might be expected from some physical models for BLQSO light curves.

In our analysis, we use the luminosity density at 1350 Å in order to minimize bias at the highest redshifts introduced from extrapolating the power-law continuum beyond the spectral window. At z ≳ 1.8, the rest frame λ = 1350 Å falls within the observed spectral window for the SDSS sources used in this work, while a smaller redshift window is available when using the luminosity densities calculated at other wavelengths. Furthermore, the z ∼ 1 distributions of bolometric luminosity did not exhibit any significant difference when using the luminosity density at 1350 Å, as compared to that calculated using the luminosity density at λ>1350 Å, suggesting that significant biases at lower z are not introduced by extrapolating the power-law continuum.

We can connect the parameters in Equation (2) to the distribution of the Eddington ratio and bolometric correction. The monochromatic luminosity is related to the Eddington luminosity ratio Γ_Edd ⁹ and bolometric correction C_λ as

$$\begin{equation} \lambda L_{\lambda } = 1.3 \times 10^{38} \frac{\Gamma _{{\rm Edd}}}{C_{\lambda }} \frac{M_{\rm BH}}{M_{\odot }} \;({\rm erg\;s^{-1}}). \end{equation} \tag{ 3 }$$

Therefore, Equation (2) implies that for the jth log-normal distribution we are assuming that on average log(Γ_Edd/C_λ) = α_0,j − 47.11 + (α_m,j − 1)log M_BH, with a Gaussian scatter about this mean value of standard deviation σ_l,j. We do not make any formal attempt to prohibit Equation (2) from allowing values of L/L_Edd > 1, as this would require us to make an assumption about the bolometric correction. However, as we will show in Section 4.5, our estimate for p(L_λ|M_BH) implies only a negligible fraction of BLQSOs with L/L_Edd > 1, assuming a constant bolometric correction of C₁₃₅₀ = 4.3 (Vestergaard & Osmer 2009).

The distribution of emission line widths at a given luminosity density and black hole mass is modeled as a log-normal distribution:

$${\selectfont\fontsize{9}{10}{\begin{eqnarray} &&\lefteqn{p(\dsty\log {\rm FWHM}_{\rm BL}|L_{\rm BL},M_{\rm BH})} \nonumber\\ &&= \dsty\frac{1}{\sqrt{2 \pi \big(\sigma _{\rm BL}^2 + \sigma ^2_{\rm FWHM}\big)}} \nonumber \\ &&\quad\!\!\times\!\exp\!\! \left\lbrace {-}\!\frac{1}{2} \!\frac{\big(\log {\rm FWHM}_{\rm BL} - \beta ^{\rm BL}_0 + 1/4 \log L_{\rm BL} - 1/2 \log M_{\rm BH}\big)^2}{\sigma ^2_{\rm BL} + \sigma ^2_{\rm FWHM}}\! \right\rbrace{.}\nonumber\\ \end{eqnarray}}} \tag{ 4 }$$

Here, FWHM_BL is the line width for a particular broad emission line, L_BL is the luminosity density used as an estimate for the broad-line region size for a particular broad emission line, σ_FWHM is the measured uncertainty in FWHM due to measurement error, and β^BL₀ and σ_BL are the parameters for Equation (4) for a particular broad emission line. We do not make any attempt to correct for radiation pressure on the broad emission line clouds (Marconi et al. 2008), as its importance and effects are currently poorly understood (e.g., for a discussion see Vestergaard & Osmer 2009).

Equation (4) implies that on average FWHM ∝ M^1/2_BH/L^1/4_BL, or equivalently M_BH ∝ L^1/2_BLFWHM²_BL. Therefore, we can use the results obtained for the broad emission line scaling estimates of M_BH to fix β₀. We use the mass scaling relationship for C iv that is presented by Vestergaard & Peterson (2006), and the relationship for Mg ii that is presented by Vestergaard & Osmer (2009). As noted in KVF09, these scaling relationships imply β^Mg ii₀ = 10.61 and $\beta ^{{\rm C\,\mathsc {iv}}}_0 = 11.33$, and we fix β₀ to these values.

In addition, the dispersion in FWHM_BL at a given L_BL and M_BH can be related to the scatter in the mass estimates based on the scaling relationships. Vestergaard & Peterson (2006) find the statistical uncertainty in the broad-line mass estimates to be 0.36 dex for C iv. Therefore, because FWHM ∝ M^1/2_BH, $\sigma _{\rm BL}({\rm C\,\mathsc {iv}}) \approx 0.18$ dex. Likewise, the intrinsic uncertainty in the broad-line mass estimate for Mg ii is found to be ∼0.4 dex (M. Vestergaard et al. 2010, in preparation), and therefore $\sigma _{\rm BL}({\rm Mg\,\mathsc {ii}}) \approx 0.2$ dex. However, there have been indications from previous analysis of flux-limited surveys, which probe higher redshifts and a narrower range in luminosity than that exhibited by the objects with reverberation mapping data, that the intrinsic scatter in the mass estimates may be smaller than ≈0.4 dex (Kollmeier et al. 2006; Shen et al. 2008; Fine et al. 2008; Steinhardt & Elvis 2010a, 2010b). This may be caused by correlated scatter in the mass estimates with luminosity (e.g., Shen & Kelly 2010) or redshift, or a dependence on σ_BL with luminosity or redshift. Both of these possibilities would decrease the dispersion in the scatter in the mass estimates when only probing a smaller range in L, or higher redshifts. In addition, Marconi et al. (2008) find a smaller scatter in the masses estimated using Hβ when one corrects the reverberation mapped masses for radiation pressure. Because of the possibility for a smaller scatter in the mass estimates, we perform our analysis with both σ_BL fixed to ≈0.2 dex, and with σ_BL as a free parameter.

The technique for estimating the BHMF developed by KVF09 takes a Bayesian approach for performing statistical inference, meaning that it calculates the probability distribution of the mass function, given the observed data. All the information regarding the BHMF and the parameters for the statistical model is contained within the posterior probability distribution, which is defined as the probability distribution of the model, given the observed data. KVF09 derived the posterior distribution for the statistical model described in the previous section. Denoting the model parameters for the shape of the BHMF as θ = (π, μ, Σ, γ, α₀, α_m, σ_l), the posterior distribution is

$$\begin{eqnarray} p(\theta |{\rm FWHM},L,z) &\propto & p(\theta) \left[p(I=1|\theta) \right]^{-n}\nonumber\\ && \times\prod _{i=1}^n p({\rm FWHM}_i,L_{\lambda,i},z_i|\theta), \end{eqnarray} \tag{ 5 }$$

where the number of data points is n, p(θ) is the prior on θ, and p(I = 1|θ) is the probability as a function of θ that a BLQSO makes it into the SDSS DR3 catalog. We note that if the dispersion in the mass estimates, σ_BL, is also treated as a free parameter, then θ also contains σ_BL. The joint distribution of FWHM, L, and z, p(log FWHM_i, log L_λ,i, log z_i|θ), is given by Equations (30)–(40) in KVF09, modified to use the mixture form for p(L_λ|M_BH). We use the Mg ii line at 1 < z < 1.6, both the Mg ii and the C iv line at 1.53 < z < 1.6, and the C iv line at z > 1.6. We do not use Hβ because we limit our analysis to z > 1. At z ∼ 0.8 Hβ is shifted into the water vapor bands, which tends to decrease the FWHM accuracy; Mg ii is similarly affected at higher redshifts. In addition, we see systematic changes in the Mg ii FWHM distribution above a redshift of 1.6, suggesting the presence of biases, which need further investigation (M. Vestergaard et al. 2010, in preparation). The posterior distribution for the BHMF normalization, given θ and n, is given by Equation (16) in KVF09.

The inclusion probability as a function of θ is calculated by averaging the SDSS DR3 selection function over the distribution of L_λ and z (see Equations (46)–(49) in KVF09). In order to simplify our analysis, we ignore the lower limit of FWHM>1000 km s⁻¹ on the emission line width for the SDSS DR3 sample. The distribution of FWHM for the SDSS DR3 falls off before reaching FWHM = 1000 km s⁻¹, and it does not appear that imposing the lower-limit results in a non-negligible fraction of the BLQSO population being missed. Therefore, any correction for the lower limit in FWHM is negligible, and we ignore it. In this case, the inclusion probability is

$$\begin{eqnarray} \hspace{-12pt}p(I=1|\theta) & = & \frac{\Omega }{4\pi } \int _{L_{\lambda }=0}^{\infty } \int _{z=1}^{4.5} \frac{s(L_{\lambda },z)}{L_{\lambda } z \ln 10} \nonumber \\ & & \times \sum _{k=1}^{K} \pi _k \left[ \sum _{j=1}^J \gamma _j N(\log L_{\lambda }|\bar{l}_{kj}(z),V_{l,kj})\right] \nonumber\\ &&\times N (\log z|\mu _{z,k},\sigma ^2_{z,k})\ dz\ dL_{\lambda }, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \bar{l}_{kj}(z) & = & \alpha _{o,j} + \alpha _{m,j} \mu _{m,k} + \frac{\alpha _{m,j} \sigma _{mz,k}}{\sigma _{z,k}^2} \left(\log z - \mu _{z,k} \right), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} V_{l,kj} & = & \alpha _{m,j}^2 \sigma ^2_{m,k} \big(1 - \rho _{mz,k}^2\big) + \sigma ^2_{l,j}. \end{eqnarray} \tag{ 8 }$$

Here, Ω = 1622 deg² is the effective sky area of the SDSS DR3 sample (Richards et al. 2006), s(L_λ, z) is the SDSS selection function, N(x|μ, σ²) denotes a normal distribution with mean μ and variance σ², as a function of x, μ_m,k and μ_z,k are the mean values of log M_BH and log z for the kth Gaussian function, respectively, and ρ_mz,k is the correlation between log M_BH and log z for the kth Gaussian function. Note that $\bar{l}_k(z)$ and V_l,k define the mean and variance in log L_λ at a given redshift for the kth Gaussian function. The integral in Equation (6) is over 1 < z < 4.5 because we have removed the sources at z < 1, and there are no useable broad emission lines at z ≳ 4.5.

Equation (6) is in terms of the selection function with respect to the luminosity density. As mentioned above, we use the luminosity density at 1350 Å in this work. However, Richards et al. (2006) report their selection function in terms of the i-band magnitude. We can convert the selection function of Richards et al. (2006) to be in terms of the BLQSO power-law continuum L₁₃₅₀ through the equation

$$\begin{equation} s(L_{1350},z) = \int _{-\infty }^{\infty } s(i,z) p(i|L_{1350},z)\ di, \end{equation} \tag{ 9 }$$

where p(i|L₁₃₅₀, z) is the distribution of i-band magnitude at a given L₁₃₅₀ and z, and s(i,z) is the SDSS DR3 selection function in terms of i and z. We model the distribution of i-magnitudes at a given L₁₃₅₀ in different redshift bins as a student's t distribution with mean that depends linearly on log L₁₃₅₀:

$$\begin{eqnarray} &&p(i|L_{1350},z) = \frac{\Gamma [(\nu (z) + 1)/2]}{\Gamma (\nu (z)/2) \sigma _i(z) \sqrt{\nu (z) \pi }} \nonumber\\ &&\quad \times\bigg(1 + \frac{1}{\nu (z)}\bigg(\frac{i - A_i(z) - B_i(z) \log L_{1350}}{\sigma _i(z)}\bigg)^2\bigg)^{-(\nu (z) + 1) / 2}.\nonumber\\ \end{eqnarray} \tag{ 10 }$$

The student's t distribution converges to the normal distribution for ν → ∞; for finite ν the t distribution is more heavy tailed than the normal distribution, and we have found it to better describe the distribution of i-magnitudes at a given L₁₃₅₀ and z. The parameters A_i(z), B_i(z), σ_i(z), and ν(z) are fit by maximizing their posterior probability distribution, given the observed set of i-magnitudes and L₁₃₅₀. The posterior distribution is given by inserting the assumed distributions into Equation (40) of Kelly (2007), and for simplicity we assume a simple flux limit of i = 19.1 for z < 2.7 and i = 20.2 for z > 2.7 (e.g., Richards et al. 2006). At z < 2.7, the redshift bins used in Equation (10) have width Δz = 0.1, while at z > 2.7 the redshift bins have width Δz = 0.3.

In this work, we constrain the parameters of our statistical model to be within certain limits, but in general assume a uniform prior on θ. Our prior is uniform on the parameters within the limits 44.5 < α₀ < 46.5, − 1 < α_m < 3, 0.1 < σ_l < 2, 7 < μ_m,k < 10, log 1 < μ_z,k < log 4.5, 0.1 < σ_m,k, σ_z,k < 2, and −0.98 < ρ_mz,k < 0.98. We also assume a uniform prior on π and γ, subject to the elements of π and γ summing to unity. The limits on α_m were chosen because we did not consider it realistic that L₁₃₅₀ would depend on M_BH outside of the range of dependences spanning L₁₃₅₀ ∝ 1/M_BH to L₁₃₅₀ ∝ M³_BH, and the limits on α₀ were chosen to restrict the average value of L/L_Edd to be within 0.01 < L/L_Edd < 1 for BLQSOs with M_BH = 10⁹ M_☉, assuming a bolometric correction of C₁₃₅₀ = 4.3. The limits on σ_l were chosen because Δlog L ≈ 0.1 is comparable to the grid spacing on which the selection function was computed by Richards et al. (2006), and we did not consider it likely that the dispersion in L₁₃₅₀ at a given M_BH would be greater than 2 dex. The limits on the BHMF parameters were chosen so as not to extrapolate the BHMF very far beyond the detectable range of M_BH. As such, in this work we limit our analysis of the BHMF to M_BH ≳ 10⁷ M_☉.

As mentioned before, there exists the possibility that, in the range of L and z we probe, the uncertainty in the mass estimates may be smaller than the commonly quoted ∼0.4 dex. If the error in the mass estimates is correlated with luminosity, then the variance in the mass estimates at a given luminosity and mass is reduced to

$$\begin{equation} {\rm Var}(\log \hat{M}_{\rm BL}|L) = {\rm Var}(\log \hat{M}_{\rm BL}) \big(1 - \rho _{\rm BL}^2\big). \end{equation} \tag{ 11 }$$

Here, $\hat{M}_{\rm BL}$ denotes the broad-line mass estimate, ${\rm Var}(\log \hat{M}_{\rm BL})$ is the variance in the mass estimates about the true mass, typically thought to be ∼0.4² dex², and ρ_BL is the correlation with luminosity in the scatter in the mass estimates about the true mass. Because we probe a somewhat narrow range in L_λ, when we estimate the parameter σ_BL in Equation (4), we are really estimating $2\sigma _{\rm BL} \approx \sqrt{{\rm Var}(\log \hat{M}_{\rm BL}|L)}$, and we therefore need to also construct a prior for σ_BL. We do this by first noting that Vestergaard & Peterson (2006) estimated the dispersion in the scatter in the broad-line mass estimates about the reverberation mapping estimates using 27 AGNs. Therefore, the appropriate prior distribution for Var$(\log \hat{M}_{\rm BL})$ is the posterior probability distribution of the variance in the mass estimates, given the reverberation mapping sample. Since we assume that the scatter in the mass estimates about the true mass is log-normal, the probability distribution of their variance follows using a standard result from Bayesian statistics, and is a scaled inverse χ² distribution with 26 degrees of freedom (e.g., Gelman et al. 2004). However, there are currently no constraints on the value of ρ_BL, and we use a uniform prior on its value. Our prior on σ_BL is then calculated by combining these two priors according to Equation (11). This results in a broad prior which peaks at Var$(\log \hat{M}_{\rm BL}|L) \approx 0.4^2\;{\rm dex}^2$, and falls off slowly to zero as Var$(\log \hat{M}_{\rm BL}|L) \rightarrow 0$ and Var$(\log \hat{M}_{\rm BL}|L) \rightarrow 0.6^2$.

In addition to the above constraints, we also impose the prior constraint that the number density of BLQSO SMBHs must never exceed the local number density of all SMBHs. In principle, this constraint could be violated if a large number of "wandering" black holes are present. These wandering black holes would be ejected from galactic nuclei in the late stages of a merger due to asymmetric emission of gravitational radiation (e.g., Volonteri 2007). However, this constraint would only be violated if the binary black hole system is ejected after or during the BLQSO phase. While this would certainly be a very intriguing result, we do not test it here; indeed, our results are unaffected by this prior constraint as the estimated BHMF is always below the local value for all random draws from the posterior probability distribution.

As described in KVF09, we do not work with Equation (5) directly, but instead use a Markov Chain Monte Carlo (MCMC) sampler algorithm to obtain a set of random draws of θ, distributed according to Equation (5). Because the posterior distribution described by Equation (5) is multimodal due to ambiguities in labeling the Gaussian functions, we label the BHMF Gaussian functions in order of increasing implied mean flux, and the p(L|M_BH) Gaussian functions in order of increasing mean luminosity at M_BH = 10⁹ M_☉. In addition, in order to make our MCMC sampling algorithm robust against additional possible multimodality, we include parallel tempering (e.g., Liu 2004) in our MCMC algorithm to facilitate sampling from the different modes. For each value of θ, we obtained from our MCMC random number generator, we also obtain a value of the BHMF normalization, N, by drawing from a negative binomial distribution with parameters n and p(I = 1|θ) (KVF09). The random realizations of θ and N then define a sample of random realizations of the BHMF obtained from the posterior probability distribution of the BHMF, given the observed set of mass estimates, luminosity densities, and redshifts. We ran our MCMC sampler for 2 × 10⁵ iterations each, keeping every 20th iteration.

Although we have assumed a flexible parametric form for the BLQSO BHMF and p(L|M_BH), this form is inconsistent with more physically motivated BLQSO p(L|M_BH), such as might be expected from a power-law decay in the BLQSO accretion rate (e.g., Yu et al. 2005; Hopkins & Hernquist 2006, also see discussion in Section 5.1). Instead, our assumed log-normal mixture form for p(L|M_BH) was motivated by mathematical convenience, as some of the integrals necessary for evaluation of the posterior distribution can be done analytically (KVF09). If we were to use a form which required numerical integration, we would have to perform over 10 billion numerical integrals in our MCMC sampler, which is computationally prohibitive. Moreover, we also used the mixture form to allow flexibility in the estimated p(L|M_BH), so that our assumed form should be able to approximate many different forms for p(L|M_BH).

In order to assess the impact of our chosen parametric form on the inferred BHMF, we simulated a data set where the distribution of Eddington ratios was assumed to have a power-law form p(Γ_Edd|M_BH) ∝ Γ^−(1+1/β)_Edd with β = 2 (see discussion in Section 5.1). The distribution of bolometric corrections was log-normal with geometric mean C₁₃₅₀ = 5 and dispersion of 0.2 dex, and the BHMF normalization was N = 2 × 10⁶; all other aspects of the simulation were done in the same manner as described in Section 6.1 of KVF09.

The results are illustrated in Figure 1, where we compare the true and estimated BHMF for the simulated sample at z = 2, and the true Eddington ratio distribution with that inferred assuming the statistical model described in Section 3.1. Here, and elsewhere in this paper, in addition to a single "best-fit" mass function, we also plot 100 random draws of the mass function from its probability distribution, generated by our MCMC random number generator. The spread and density of the random draws of the BHMF, and any quantities derived from it, give a visual representation of the uncertainty in these quantities, with the spread on the random draws constraining the BHMF. Because we use 100 random draws, the probability of a quantity falling within a certain area on a plot can be estimated by counting the number of random draws that intersect that area. For this example, the BLQSO BHMF inferred assuming a mixture of J = 3 log-normal distributions for p(L_λ|M_BH) is able to recover the true BHMF and Eddington ratio distributions, at least when the true distribution of L/L_Edd is p(Γ_Edd) ∝ Γ^−1.5_Edd. Although this test is far from exhaustive, we consider it reasonable to conclude that our flexible form for the BHMF and p(L_λ|M_BH) is able to accurately approximate the true forms, and therefore our results are robust against errors resulting from using an incorrect parametric form for these distributions.

Zoom In Zoom Out Reset image size

We used our Bayesian method to derive the BLQSO BHMF from the SDSS DR3 quasar sample, both holding the magnitude of the scatter in mass estimates fixed to 0.4 dex, and treating the amplitude of the scatter as a free parameter. However, before discussing the results, we first evaluate how well the statistical model described in Section 3.1 fits our sample. In order to do this, we compare the distributions of redshift, luminosity, and broad-line mass estimates of our sample to the distributions implied by our model, as described in KVF09. Figure 2 compares the observed distributions of z, λL_λ(1350 Å), and $\hat{M}_{\rm BL}$ to those implied by the statistical model described in Section 3.1. The implied distributions were calculated by first simulating a sample of M_BH and z from a BHMF randomly output from the MCMC sampler. Then, for each value of M_BH, we simulated a value of λL_λ(1350 Å) and mass estimate $\hat{M}_{\rm BL}$. Finally, we applied the selection function given by Equation (9) to the simulated data set. This was repeated for each realization of the BHMF obtained from the MCMC output, in order to account for the uncertainty in our estimated model parameters.

Zoom In Zoom Out Reset image size

We are not able to fit the mass estimate distributions if we keep the statistical scatter in the mass estimates fixed to 0.4 dex. In particular, this predicts a distribution of mass estimates that is too broad compared to the actual distribution. In fact, the observed dispersion in the mass estimates for our SDSS sample is ≈0.35 dex, smaller than that expected even if all objects had the same mass. However, if we allow the dispersion in the mass estimate error to be a free parameter, we are able to obtain a good match to the data. Our best-fit value for the scatter in the mass estimates at a given mass is 0.18 ± 0.01 dex for Mg ii and 0.13 ± 0.01 for C iv. Our result that the scatter in the mass estimates for high L and z must be smaller than ∼0.4 dex is consistent with what has been found in previous work (Kollmeier et al. 2006; Shen et al. 2008; Fine et al. 2008), and our best-fit values of the standard deviations in the mass estimate errors are consistent with the upper limits recently calculated by Steinhardt & Elvis (2010b).

The origin of this smaller scatter is unclear, and there may be several different possibilities. One possibility, as mentioned earlier and by Shen et al. (2008) and Shen & Kelly (2010), is that the error in the mass estimates may be correlated with luminosity. If this is true, then an error of ∼0.4 dex represents the error in the mass estimates when averaging over a broad range in L, as was done for the reverberation mapping sample, while a smaller scatter of ∼0.15 dex represents the error in the mass estimates when one is limited to a more narrow range in L, as the SDSS quasar sample is. It is unclear why the error in the mass estimates would be correlated with luminosity, but one possible source unaccounted for is radiation pressure. Marconi et al. (2008) argue that virial mass estimates should be corrected for radiation pressure. They find a correction that implies a steeper dependence on luminosity than $\hat{M}_{\rm BL} \propto L^{0.5}$, especially for sources with high L/L_Edd. Under their model, if one does not correct for radiation pressure then one will tend to underestimate the mass with increasing L, producing a correlation between the error in the mass estimates with luminosity, and possibly producing a smaller scatter in the mass estimate errors over a narrow range in luminosity. It is interesting to note that Marconi et al. (2008) find a smaller scatter in the mass estimates of ∼0.2 dex when correcting for radiation pressure for the reverberation mapped sample, which covers a larger range in luminosity. However, more work is need to understand the importance of radiation pressure, and if it can produce the observed smaller scatter in the mass estimates over the range in luminosity we probe.

Another possibility is that the error in the mass estimates may not be correlated with L, but the dispersion in the errors may decrease with increasing L or z. Unfortunately, the number of AGNs with M_BH estimated from reverberation mapping is too small to test this, and dominated by sources at lower L and z. The third possibility is that the virial mass estimates are biased at high L and z, at least for Mg ii and C iv, and are only marginally related to the actual masses. The R–L relationship is well established for Hβ, and does not require a large extrapolation to the luminosities probed in our sample, and thus we do not expect Hβ-based mass estimates to be significantly biased in this range (Vestergaard 2009). The situation is less clear for Mg ii and C iv. There is only one reliable time lag for the Mg ii emission line (Metzroth et al. 2006). An R–L relationship has been estimated for the C iv line over a broad range in luminosity and redshift (Kaspi et al. 2007), including those probed in our study, and the R–L relationship is similar for both Hβ and C iv. Unfortunately, the C iv R–L relationship is estimated from only eight data points, and further work is needed in order to understand the Mg ii- and C iv-based mass estimates and their errors.

Throughout the rest of this work we will focus on the results obtained from allowing the dispersion in the mass estimate error to be a free parameter, a value of ∼0.4 dex is clearly ruled out. However, we note that we have performed the same analysis for both cases, and while the quantitative details change, the scientific conclusions are unaffected by treating the amplitude of the scatter as a free parameter. The only exception is that we infer a much more narrow distribution of the Eddington ratio if we fix the amplitude of the scatter in the mass estimates to be ∼0.4 dex. In addition, the results reported in this section highlight the need for more reverberation mapping studies in order to better understand the nature of the errors in the mass estimates.

Figure 3 shows the BLQSO BHMF at several redshifts. Our estimated BHMF is compared with an estimate of the local mass function of all SMBHs, and the BHMF reported by Vestergaard et al. (2008), obtained from binning up the broad-line mass estimates. Following Merloni & Heinz (2008), the local BHMF was computed to be near the middle of the uncertainty range reported by Shankar et al. (2009) by convolving a Schechter function with a log-normal distribution with standard deviation 0.3 dex, chosen to be consistent with the scatter about the M_BH–σ relationship; the parameters for the Schechter function are those reported by Merloni & Heinz (2008). Also, we note that early-type galaxies dominate the local BHMF at M_BH ≳ 4 × 10⁷ M_☉ (Yu & Lu 2008).

Zoom In Zoom Out Reset image size

In order to focus on the region of the BHMF that is robust against uncertainties in the selection function, as well as against uncertainty on the Eddington ratio distribution, we estimate the black hole mass completeness for our sample as a function of z. Our best estimate of the SDSS completeness as a function of M_BH and z is shown in Figure 4, and the 10% completeness limit is marked by a vertical line in Figure 3. The black hole mass completeness depends on both the completeness in luminosity, and the assumed distribution of L at a given M_BH. As can be seen from Figure 4, at z ≳ 2 the SDSS quasar sample is only ≈10% complete at M_BH ∼ 10⁹ M_☉, becoming more incomplete at lower masses. At masses much lower than M_BH ∼ 10⁹ M_☉, the estimated mass function almost completely depends on extrapolation from the set of BHMFs and Eddington ratio distributions that fit the observed data well, constrained by our assumed parametric forms. Therefore, we stress that below M_BH ∼ 10⁹ M_☉ the BLQSO BHMF must be interpreted with caution, and in this work we will try to focus on what we can infer from the high-mass end of the mass function.

Zoom In Zoom Out Reset image size

We also estimate the completeness of our sample as a function of L/L_Edd, assuming a constant bolometric correction of C₁₃₅₀ = 4.3. The estimated completeness depends on the selection function and the estimated BHMF, as the selection function depends on luminosity, which is defined by the BHMF at a given L/L_Edd. The Eddington ratio completeness is also shown in Figure 4. The SDSS is only ∼50% complete for sources radiating at the Eddington limit, and ≲10% complete for sources radiating at ≲10% of Eddington. This heavy incompleteness is due to the fact that, for our estimated BLQSO BHMF, half of BLQSOs have black holes that are not massive enough to make the SDSS flux limit even if they radiate at Eddington. We further discuss the Eddington ratio distribution in Section 4.5.

The difference between the binned estimate of the BLQSO BHMF, calculated from the estimate of Vestergaard et al. (2008), and our estimate, is the result of the difference in statistical methodology employed by Vestergaard et al. (2008) and our work. First, the statistical error in the broad-line mass estimates results in a broader inferred BHMF when one simply bins up these estimates (e.g., Kelly & Bechtold 2007, KVF09). Second, the 1/V_a technique corrects for incompleteness in flux, but not black hole mass. As a result, the 1/V_a corrections only partially correct for incompleteness in M_BH, and the estimated binned BHMF will still suffer from incompleteness, especially at the low-mass end. The two effects combined result in a systematic shift in the estimated BHMF toward higher M_BH (Shen et al. 2008; KVF09). However, the Bayesian approach outlined in KVF09 is able to self-consistently correct for these two effects in a statistically rigorous manner, conditional on the survey selection function and assumptions implicit within the statistical model. In spite of these differences in methodology, our estimated BHMF agrees fairly well with that of Vestergaard et al. (2008) at the high-mass end, considering the differences in the methodology, and the two do not strongly diverge until values of M_BH where the SDSS is highly incomplete. The poorer agreement between the two estimates at z < 2 is due to the fact that the BHMF is estimated using the Mg ii emission line in this redshift range, which we find to have a higher statistical error than that of C iv. As a result, our correction to the BHMF due to the error in the mass estimates is greater at these redshifts.

In Figure 5, we show the evolution in the comoving number density of SMBHs in BLQSO at four different masses. As is evident, the number density of higher mass BLQSO SMBHs peaks at higher redshift, where the number density of BLQSO SMBHs of M_BH ∼ 5 × 10⁸ M_☉ peaks at z ∼ 1.5, and the number density for M_BH ∼ 5 × 10⁹ M_☉ peaks at z ∼ 2.5. This result is consistent with what has commonly been referred to in the literature as black hole "downsizing," where the most massive SMBHs are active, and grow, at earlier epochs than lower mass black holes, an effect also seen by Vestergaard & Osmer (2009) and Steinhardt & Elvis (2010a).

Zoom In Zoom Out Reset image size

In Figure 6, we show the comoving mass density of SMBHs that reside in BLQSOs, as a function of redshift, ρ_QSO(z). The peak in the cosmic mass density of BLQSO SMBHs occurs at z ∼ 2. Our constraint on the location of the peak in ρ_QSO(z) is consistent with the peak in the mass density derived from the Large Bright Quasar Survey and high-z sample of Fan et al. (2001a), as calculated by Vestergaard & Osmer (2009). Previous work has not found any evidence for evolution in the Eddington ratio distribution at z > 1 for the most massive BLQSOs (e.g., McLure & Dunlop 2004; Vestergaard 2004; Kollmeier et al. 2006; Vestergaard & Osmer 2009). If there is no significant evolution in the Eddington ratio distribution at z > 1 for these systems, then we would expect the peak in their luminosity density to coincide with the peak in their black hole mass density. The luminosity density of quasars peaks at z ∼ 2 (e.g., Wolf et al. 2003; Hopkins et al. 2007c), matching the peak in black hole mass density observed in our work.

Zoom In Zoom Out Reset image size

We can use the quasar BHMF to place constraints on the fraction of black holes that are seen in the BLQSO phase as a function of M_BH; this quantity is commonly referred to as the BLQSO "duty cycle." Denote the duty cycle as δ(M_BH, z). Then,

$$\begin{equation} \delta (M_{\rm BH},z) \equiv \frac{\phi _{{\rm QSO}}(M_{\rm BH},z)}{\phi _{\rm BH}(M_{\rm BH},z)}. \end{equation} \tag{ 12 }$$

Here, ϕ_BH(M_BH, z) is the BHMF of all SMBHs at a given redshift. Ignoring mergers of SMBHs, we can compute a lower limit to the duty cycle by comparing the BHMF at a certain redshift with the local SMBH number density, since ϕ_BH(M_BH, z) ⩽ ϕ_BH(M_BH, 0). In Figure 7, we compute the lower limit of the BLQSO duty cycle at z = 1 for SMBHs with M_BH > 5 × 10⁸ M_☉. The duty cycle at z = 1 is constrained to be δ ≳ 0.01 at M_BH ∼ 10⁹ M_☉, falling steeply to δ ≳ 10⁻⁵ at M_BH ∼ 10¹⁰ M_☉. The decrease in the duty cycle with increasing black hole mass may be another reflection of SMBH downsizing, with the most massive BLQSO SMBHs being active at earlier cosmic epochs. Alternatively, it may be due to the fact that the most massive SMBHs spend a shorter amount of time in the broad-line phase, as expected from some simulations of black hole feedback (e.g., Hopkins et al. 2006b). However, we note that the local number density of the most massive SMBHs is poorly constrained and subject to considerable systematic uncertainty (Lauer et al. 2007), and therefore the duty cycle of BLQSOs with M_BH ∼ 10¹⁰ M_☉ may be subject to considerable systematic error. Indeed, the observed falloff in the lower limit on the duty cycle for the most massive systems may simply be due to the fact that we are overestimating the number density of local SMBHs.

Zoom In Zoom Out Reset image size

The quasar lifetime is not well constrained empirically (e.g., Martini 2004), but we can use our estimated BLQSO BHMF to estimate the BLQSO lifetime. For simplicity, we assume that all BLQSOs of a given mass have a single lifetime, t_BL. Furthermore, if an SMBH undergoes multiple episodes of BLQSO activity, then we assume that t_BL is the same for each episode. These assumptions are unlikely to be true, but for the purposes of our work we may still think of t_BL as a "typical" BLQSO lifetime. Because the duty cycle is the probability of observing an SMBH as a BLQSO at a given M_BH and redshift we can relate t_BL to δ(M_BH, z). The probability of observing an SMBH as a BLQSO at redshift z is the product of the probability that an SMBH becomes a BLQSO along our line of sight at some point in its evolution with the probability that that SMBH is a BLQSO between cosmic ages t(z) − t_BL and t(z), normalized by the probability that that SMBH is a BLQSO at t ⩽ t(z). The normalization results from the additional assumption that if an SMBH of mass M_BH exists at t(z), and if it goes through at least one BLQSO episode as some point in its growth, then it had to have gone through at least one BLQSO episode at t < t(z) in order to grow to M_BH by t(z). This is a reasonable assumption, at least for the most massive system, since all SMBHs at t(z) had to grow from much lower mass seeds, and if BLQSO activity occurs for an SMBH at some point in its life, BLQSO activity would have occurred during this growth period. In other words, the assumption implicit in the normalization correction is that SMBH seeds of mass, say, M_BH ∼ 10⁹ M_☉, do not simply emerge and then undergo BLQSO activity at a cosmic epoch later than t(z). The practical result of this assumption is that the BLQSO duty cycle of, say, SMBHs with M_BH ∼ 10⁹ M_☉, can decrease from δ ∼ 1 at z ∼ 7 to δ ∼ 10⁻³ at z ∼ 1.

Denote the time that a BLQSO phase is initiated in an SMBH as τ. Then, the probability that an SMBH BLQSO is seen between t(z) − t_BL and t(z) is the same as the probability that the time that a BLQSO phase was initiated for that source occurred at t(z) − t_BL < τ < t(z). Note that this allows for multiple BLQSO episodes, as multiple phases of BLQSO episodes simply alter the probability of observing an SMBH as a BLQSO at a certain redshift. The duty cycle is thus related to t_BL as

$${\selectfont\fontsize{9}{10}{\begin{equation} \delta (M_{\rm BH},z) = {\rm Pr}({\rm BLQSO}|M_{\rm BH}) \left[ \frac{ \int _{t(z) - t_{\rm BL}}^{t(z)} p(\tau |M\rm _{\rm BH},{\rm BLQSO})\ {\it d}\tau }{ \int _{0}^{t(z)} p(\tau |M_{\rm BH},{\rm BLQSO})\ d\tau } \right], \end{equation}}} \tag{ 13 }$$

where, Pr(BLQSO|M_BH) is the probability that an SMBH with a mass of M_BH is a BLQSO at some point in its evolution, and p(τ|M_BH, BLQSO) is the probability distribution of BLQSO episode initiation time for a given black hole mass. Note that Pr(BLQSO|M_BH) is not the probability that an object is currently seen as a BLQSO, but is the probability that it appears as a BLQSO to an observer on Earth at some point in its life. As mentioned above, p(τ|M_BH, BLQSO) is sufficiently general to include both multiple and single episodes of BLQSO activity, although the actual form of p(τ|M_BH, BLQSO) will depend on the distribution of the number of BLQSO episodes an SMBH goes through. All values of M_BH in Equation (13) refer to the mass of the SMBH at redshift z.

If the distribution of τ does not evolve significantly over a time t_BL, then the integral in the numerator of Equation (13) can be approximated as p(τ|M_BH, BLQSO)t_BL. Likewise, if p(τ|M_BH, BLQSO) is approximately constant over t_BL, then the distribution of BLQSO initiation times will be similar to the distribution of BLQSOs at the time we observed them. This is a reasonable assumption so long as t_BL is short compared to the timescale for a significant change in the process that initiates BLQSO activity (e.g., mergers and other fueling events). Making this assumption, and rearranging Equation (13), it follows that

$$\begin{equation} t_{\rm BL} \approx \frac{\delta (M_{\rm BH},z) \int _{0}^{t(z)} p(t^{\prime }|M_{\rm BH},{\rm BLQSO})\ dt^{\prime }}{p(t(z)|M_{\rm BH},{\rm BLQSO}) {\rm Pr}({\rm BLQSO}|M_{\rm BH})}, \end{equation} \tag{ 14 }$$

where p(t(z)|M_BH, BLQSO) is the probability distribution of BLQSOs as a function of cosmic age, t(z). The term p(t|M_BH, BLQSO) is related to the BHMF for BLQSOs according to

$$\begin{equation} p(t(z)|M_{\rm BH},{\rm BLQSO}) = \left|\frac{dt}{dz}\right|^{-1} \left[\frac{\phi (M_{\rm BH},z)}{\int _{0}^{\infty } \phi (M_{\rm BH},z)\ dM_{\rm BH}}\right]. \end{equation} \tag{ 15 }$$

From Equation (14), it follows that for the special case where all SMBHs experience a BLQSO phase, and where BLQSOs initiation times are uniformly distributed over the age of the universe, then t_BL = δ(M_BH, z)t(z). This is the definition of a quasar "lifetime" commonly found in the literature; however, as the assumption of a uniform distribution of BLQSO initiation times is inconsistent with the BLQSO BHMF or luminosity function, t_BL will in general not equal δ(M_BH, z)t(z) (see also Hopkins & Hernquist 2009).

Because we have assumed that t_BL is the same for all BLQSOs of a given mass, and therefore must be the same at all redshifts, Equation (14) may be calculated at any redshift. However, in reality we can only estimate a lower limit to t_BL. This is because we do not know the fraction of SMBHs that will experience a BLQSO phase, although if the SMBHs growth is self-regulated we might expect Pr(BLQSO|M_BH) ≈ 1. However, if some SMBHs of mass M_BH can never be seen as BLQSOs, possibly due to orientation-dependent obscuration, then Pr(BLQSO|M_BH) < 1. In addition, as discussed above, we can only calculate a lower limit to δ(M_BH, z) by comparing the BLQSO BHMF at redshift z with the local BHMF of all SMBHs. Moreover, implicit in these calculations is the assumption that none of these quantities depend strongly enough on mass to vary significantly during the growth that occurs in the BLQSO phase. In Figure 8, we show the probability distribution of the BLQSO age of SMBHs with M_BH = 10⁹ M_☉, calculated from Equation (14) at z  = 1. We calculate t_BL at M_BH = 10⁹ M_☉ because our sample is reasonably complete at this mass, especially at this redshift. In addition, we chose z = 1 because the local BHMF better approximates the z  = 1 BHMF than the BHMF at z > 1, therefore giving the tightest lower bound on the duty cycle at z  = 1. We estimate t_BL = 150 ± 15 Myr as a lower bound on the age of the BLQSO phase for SMBHs of M_BH = 10⁹ M_☉. This estimate is roughly consistent with other estimates of quasar lifetimes (e.g., Yu & Tremaine 2002; Martini 2004; Gonçalves et al. 2008; Hopkins & Hernquist 2009).

Zoom In Zoom Out Reset image size

The probability distribution for the maximum mass of an SMBH in a BLQSO may place important constraints on models of SMBH growth, and may be calculated directly from the BHMF. In the context of our work, we do not consider the maximum mass of an SMBH to be caused by a hard upper limit, above which it is impossible to make a more massive black hole, but rather the result of a finite number of black holes drawn from a mass function. The probability that the maximum mass of an SMBH in a sample of N BLQSOs is less than ${\newmathcal M}$ is simply given by the probability that all N SMBHs have $M_{\rm BH} < {\newmathcal M}$:

$$\begin{equation} {\rm Pr}\big(M_{\rm BH}^{{\rm Max}} < {\newmathcal M}|N\big) = \left[ {\rm Pr}(M_{\rm BH} < {\newmathcal M}) \right]^{N}. \end{equation} \tag{ 16 }$$

Note that N is the total number of BLQSOs that could be observed in an all-sky survey with no flux limit; i.e., N is the normalization of the BLQSO BHMF. The term ${\rm Pr}(M_{\rm BH} < {\newmathcal M})$ is calculated from the BHMF as

$$\begin{equation} {\rm Pr}(M_{\rm BH} < {\newmathcal M}) = \frac{1}{N} \int _{0}^{\newmathcal M} \int _{0}^{\infty } \phi (M_{\rm BH},z) \left(\frac{dV}{dz}\right) \ dz\ dM_{\rm BH}. \end{equation} \tag{ 17 }$$

The probability distribution of M^Max_BH is then found by differentiating Equation (16) with respect to ${\newmathcal M}$ and evaluating the result at ${\newmathcal M} = M_{\rm BH}^{{\rm Max}}$:

$${\selectfont\fontsize{9}{10}{\begin{equation} p\big(M_{\rm BH}^{{\rm Max}} | N\big) = \big[ {\rm Pr}\big(M_{\rm BH} < M_{\rm BH}^{{\rm Max}}\big) \big]^{N-1} \int _{0}^{\infty } \phi (M_{\rm BH},z) \left(\frac{dV}{dz}\right) \ dz. \end{equation}}} \tag{ 18 }$$

The posterior probability distribution for M^Max_BH, given the observed data, can be calculated by averaging Equation (18) over the MCMC output.

In Figure 9, we show the posterior probability distribution of the maximum SMBH in a BLQSO at 1 < z < 4.5. Here, M^Max_BH should be interpreted as the maximum mass of an SMBH in a BLQSO that would be observed in an all-sky survey without a flux limit over the redshift interval 1 < z < 4.5. Thus, M^Max_BH is a lower bound on the mass of the most massive SMBH in the universe. In addition, in Figure 9 we also show the probability distribution of the redshift at which the quasar with M^Max_BH would be found. As can be seen from these figures, the maximum mass of an SMBH in a BLQSO at 1 < z < 4.5 is M^Max_BH ∼ 3 × 10¹⁰ M_☉. The probability distribution for the redshift of M^Max_BH is rather broad, but we constrain the redshift for the BLQSO hosting this SMBH to be z ≳ 2.

Zoom In Zoom Out Reset image size

Our results are in agreement with what others have found using samples of broad-line mass estimates (e.g., Vestergaard 2004; Vestergaard et al. 2008; Netzer et al. 2007), although Labita et al. (2009a, 2009b) find a smaller value of M^Max_BH ∼ 5 × 10⁹ M_☉. However, in the model of Labita et al. (2009a) M^Max_BH is a parameter determining the shape of the mass function, and is not the actual maximum black hole mass of a sample of objects drawn from the distribution of black hole mass; i.e., there is nothing in the statistical model of Labita et al. (2009a) that prevents objects with M_BH > M^Max_BH. As such, the actual realized maximum mass in a large sample of objects will be greater than the value of M^Max_BH estimated using the model of Labita et al. (2009a), as we have found. Considering this, our results are consistent with those of Labita et al. (2009a, 2009b).

These highly massive black holes represent the extremes of black hole growth, and thus are important for constraining models of black hole growth. Furthermore, these massive black holes offer the best chance of probing the faint end of the BLQSO Eddington ratio distribution. Therefore, it is of use to investigate how large and deep a survey must be in order to find an useful number of these objects. In Figure 10, we show the expected number of high-mass SMBHs in BLQSOs at 1 < z < 4.5 in a 1 deg² survey as a function of limiting i-magnitude, assuming our best-fit statistical model. We stress that these are the expected number counts for the true black hole masses, and not the mass estimates. Because the errors in the mass estimates will scatter more sources into higher mass bins than into lower mass bins, the number of sources with estimated mass in a mass bin will be larger than the actual number of sources.

Zoom In Zoom Out Reset image size

Any survey with a flux limit of i < 19 should be able to detect BLQSOs with M_BH ≳ 5 × 10⁹ M_☉; however, the survey must have an area of Ω ≳ 10 deg² to expect to detect at least one of these objects. Similarly, in order to get a large number of BLQSOs with M_BH ≳ 5 × 10⁸ M_☉, one needs a deep, large area spectroscopic survey. For example, the BLQSO spectroscopic samples from the COSMOS survey (Scoville et al. 2007) cover ≈2 deg² at i < 24 (Trump et al. 2007) and i < 22.5 for zCOSMOS (Merloni et al. 2010), and, ignoring redshift incompleteness, are therefore expected to contain ∼60 BLQSOs at 1 < z < 4.5 with masses 5 × 10⁸ < M_BH/M_☉ < 5 × 10⁹, but none with masses M_BH > 5 × 10⁹ M_☉.

As discussed in Section 3.1, we also estimate the distribution of the ratio of the Eddington ratio to the bolometric correction at 1350 Å, Γ_Edd/C₁₃₅₀. Under the assumption that Γ_Edd/C₁₃₅₀ follows a mixture of log-normal distributions, with mean values that linearly depend on log M_BH, and that C₁₃₅₀ = 4.3 for all sources (Vestergaard & Osmer 2009), our model implies that the geometric mean Eddington ratio increases with M_BH according to

$$\begin{equation} \left\langle \frac{L}{L_{{\rm Edd}}} \right\rangle _{\rm geo} = 0.12 \pm 0.01 \left(\frac{M_{\rm BH}}{M_{\odot }}\right)^{0.48 \pm 0.04}. \end{equation} \tag{ 19 }$$

The dispersion in L/L_Edd decreases with increasing M_BH, having a value of ∼0.4 dex at M_BH ∼ 10⁸ M_☉, and decreasing to ∼0.3 dex at M_BH ∼ 10⁹ M_☉. Both the dispersion and slope are smaller than what was found by KVF09, who analyzed a sample of BLQSOs from the Bright Quasar Survey at z < 0.5, suggesting possible evolution in the Eddington ratio distribution. However, the statistical significance of a difference in the slopes is marginal at best, as the differences are only significant at ∼2σ.

In Figure 11, we show the implied distribution of BLQSO Eddington ratios assuming a constant bolometric correction of C₁₃₅₀ = 4.3 (Vestergaard & Osmer 2009). As can be seen, the estimated distribution of quasar Eddington ratios peaks at L/L_Edd ∼ 0.05. However, we note that the location of the Eddington ratio peak falls below the 10% completeness limit of the SDSS, and thus the exact location of the peak is highly uncertain. Our inferred Eddington ratio distribution is broader and shifted toward smaller values of L/L_Edd than what has been found in previous works, which were not able to fully correct for incompleteness in black hole mass and Eddington ratio (e.g., McLure & Dunlop 2004; Vestergaard 2004; Kollmeier et al. 2006). However, it is consistent with what Shen et al. (2008) found using a similar approach to ours. The major differences between our work and that of Shen et al. (2008) are that they assume a power-law distribution of M_BH, and they constrain the BHMF and Eddington ratio distribution by visually matching the model distributions to the observed distributions. Our results, as well as the results of Shen et al. (2008), show that the narrow distribution in quasar Eddington ratios seen in previous work, and peaking at L/L_Edd ∼ 0.25, is due to uncorrected incompleteness. Indeed, deeper samples of BLQSOs have observed Eddington ratio distributions that are broader and peak at lower values (Gavignaud et al. 2008; Trump et al. 2009). Therefore, we conclude that most BLQSOs are not radiating at or near the Eddington limit and that there is a large dispersion in the Eddington ratio for BLQSOs.

Zoom In Zoom Out Reset image size

Recently, Steinhardt & Elvis (2010a) have analyzed the SDSS DR5 sample of mass estimates calculated by Shen et al. (2008), and concluded that there is a dearth of objects at a high Eddington ratio and that there is a redshift-dependent systematic decrease in the Eddington ratio with increasing black hole mass for BLQSOs in the high L/L_Edd tail of the distribution. They have called this effect a sub-Eddington boundary. While we also find evidence that BLQSOs radiating near the Eddington limit are rare, as did KVF09, the dependence on black hole mass is seemingly in contrast to what we find here, in that we find that the average Eddington ratio increases with M_BH, assuming a constant bolometric correction. The most important difference between their work and ours is the fact that Steinhardt & Elvis (2010a) analyze separate redshift bins, while we do not model a redshift dependence in the Eddington ratio distribution. Steinhardt & Elvis (2010a) did not perform any calculations for the whole sample, averaged over all redshifts, making a more direct comparison difficult, and thus it is unclear how discrepant our conclusions are.

Another potential contributor to this discrepancy is the handling of the statistical error in the mass estimates. Statistical errors (e.g., measurement error) have the effect of flattening slopes and correlations when one analyzes the quantity contaminated by the error (e.g., Akritas & Bershady 1996; Kelly 2007; Kelly & Bechtold 2007), due to the fact that the distribution of the estimated quantity is a biased estimate of the distribution of the quantity of interest. Because the mass estimates are contaminated by statistical error, the dependence of luminosity on the mass estimates will be flatter than the intrinsic dependence of luminosity on M_BH. In particular, the statistical error will cause some masses to be underestimated, and some to be overestimated. The objects which appear to have the highest masses at a given luminosity will have the highest overestimates, and thus will have their Eddington ratios underestimated. As a result, there will be an apparent dearth of high Eddington ratio objects at the highest masses. Similarly, there will be an apparent excess of high Eddington ratio objects in the low-mass bins. Our method estimates the intrinsic dependence of luminosity on M_BH by, in a sense, "deconvolving" the distribution of the mass estimates and luminosity with the distribution of the error in the mass estimates; this is also why we find a somewhat steeper decline in the mass function at the highest masses, as compared to the estimate reported by Vestergaard et al. (2008).

Although the statistical errors in the mass estimates can have the effect described above, possibly contributing to the different conclusions, it is unlikely that the statistical errors alone are sufficient to account for the discrepancy, and other possibilities include differences in the methods and scaling relationship employed for estimating the masses, and other differences in the methods of data analysis employed. However, further investigation is beyond the scope of our work. In addition, we note that our samples only overlap at 1 < z < 2,¹⁰ and we cannot compare with the low-redshift results of Steinhardt & Elvis (2010a), where they find the greatest evidence for these effects.

Although we have used a flexible form for p(L|M_BH), our estimated distribution for L/L_Edd may be affected by significant systematics, as it relies on the assumption of a constant bolometric correction. There is evidence that the bolometric correction depends on both Eddington ratio (Vasudevan & Fabian 2007; Young et al. 2010) and black hole mass (Kelly et al. 2008), and the increase of L/L_Edd with M_BH we find may instead be a reflection of a decrease in the bolometric correction with M_BH. Indeed, Kelly et al. (2008) find that the ratio of optical/UV flux to X-ray flux increases with M_BH, implying that the bolometric correction to the UV flux decreases with M_BH, and therefore we would infer an increase in L/L_Edd with M_BH assuming a constant bolometric correction, as we do. In addition to these issues, we also note that the SDSS is at least partially incomplete at all Eddington ratios at z > 1,¹¹ as shown in Figure 4, and further work is needed using deeper surveys to confirm these results.

A significant amount of recent work has suggested that quasar feedback regulates the growth of SMBHs and affects the large-scale evolution of its host galaxy. Understanding the quasar light curve is thus of fundamental importance for understanding black hole growth and feedback, as well as placing constraints on the physics of quasar accretion flows. The distribution of luminosity at a given black hole mass can be related to the quasar light curve, and therefore one can use the estimated distribution of luminosity at a given black hole mass to place some empirical constraints on models for quasar light curves. We address this in the following section, and argue that our estimated Eddington ratio distribution is consistent with models where the BLQSO phase represents the final stages of an SMBH's growth, in which the QSO light curve is expected to decay until the object no longer appears as a BLQSO.

For a given fueling episode, denote the quasar light curve as L(t), and the mass fueling rate to the SMBH as a function of time as $\dot{M}(t)$. The fueling rate, $\dot{M}(t)$, gives the rate at which matter is externally supplied to the BLQSO accretion disk. The probability distribution of quasar luminosity at a given $M_{\rm BH}, \dot{M},$ and t is $p(L|M_{\rm BH},\dot{M},t)$, and the probability distribution of the fueling rate at a given M_BH and t is $p(\dot{M}|M_{\rm BH},t)$. The distribution $p(L|M_{\rm BH},\dot{M},t)$ is determined by the physics of the accretion disk, since the quasar luminosity is produced by viscous stresses in the accretion disk. For example, Siemiginowska & Elvis (1997) calculate $p(L|M_{\rm BH},\dot{M},t)$ implied by the model of Siemiginowska et al. (1996) for a thermal-viscous accretion disk stability. The distribution $p(\dot{M}|M_{\rm BH},t)$, on the other hand, depends on the physics and stochastic nature of the fueling mechanism and quasar feedback.

Because the quasar luminosity is determined by the physics of the accretion disk, it is reasonable to assume that given M_BH and $\dot{M}$, L is independent of time. This does not imply that the quasar light curve does not vary with time, as it certainly does, but rather that knowing the value of t does not give us any additional information on the value of L when we already know M_BH and $\dot{M}$ at a given t. We therefore drop the explicit dependence of $p(L|M_{\rm BH},\dot{M},t)$ on time.

The distribution of luminosity at a given M_BH is calculated as

$$\begin{eqnarray} p(L|M_{\rm BH}) &=& \int _{0}^{\infty } p(L|\dot{M},M_{\rm BH}) \nonumber\\ &&\times\left[\int _{t_0}^{t_o+t_{\rm BL}} p(\dot{M}|t,M_{\rm BH}) p(t|M_{\rm BH})\ dt \right] \ d\dot{M},\nonumber\\ \end{eqnarray} \tag{ 20 }$$

where p(t|M_BH) is the probability of seeing a BLQSO at time t, given M_BH, and t_BL is the length of time that a quasar would be classified as a BLQSO. Here, we have defined the broad-line phase of quasar activity to start at t = t₀. The term p(t|M_BH) is related to the model for black hole growth, since p(t|M_BH) ∝ p(M_BH|t)p(t). Because we observe quasars randomly during their BLQSO phase, p(t) is uniform and p(t|M_BH) ∝ p(M_BH|t). In general, BLQSOs with larger values of M_BH are more likely to be seen later in their lifetime, i.e., at larger values of t.

As an alternative to Equation (20), p(L|M_BH) may be directly calculated from the amount of time that a BLQSO spends above a given luminosity, as a function of M_BH. In this case, the term p(L|M_BH) is directly related to the "luminosity-dependent" lifetime interpretation advocated by Hopkins et al. (2005b, 2005c), where the quasar "lifetime" is the amount of time that a quasar spends above a given luminosity. Assume that p(t|M_BH) ∝ 1 for BLQSOs, and denote T(L|M_BH) to be the amount of time a BLQSO spends above a luminosity L, as a function of M_BH. Then,

$$\begin{equation} p(L|M_{\rm BH}) \propto \left| \frac{dT(L|M_{\rm BH})}{dL} \right|. \end{equation} \tag{ 21 }$$

For example, Hopkins & Hernquist (2009) suggest that p(L|M_BH) can be well approximated as a Schechter function, which generalizes the distributions implied by several common models of quasar light curves often found in the literature.

5.1.1. Effects of Evolution in the Quasar Fueling Rate

Recent work on quasar fueling and feedback suggests that the BLQSO phase occurs at the end of the fueling event, during which feedback energy from the AGN is able to unbind the surrounding gas, removing the obscuring material and halting the black hole's growth. Within the context of this model, the SMBH does not experience a large fractional increase in mass beyond the value of M_BH that is observed during its time as a BLQSO. Therefore, we approximate p(t|M_BH) as being uniform and independent of M_BH, i.e., p(t|M_BH) = 1/t_BL. In addition, this model predicts that during this so-called blow-out phase, the fuel supply is set by the evolution of the blast wave caused by the injection of energy from the SMBH, and is expected to be of the form $\dot{M}(t) \propto t^{-\beta }$ (Hopkins & Hernquist 2006; Hopkins et al. 2006b). Alternatively, if the fuel supply is suddenly removed, then the accretion rate during the broad-line phase is due to evolution of the viscous accretion disk. The fueling rate also has a power-law form under viscous evolution of the disk, but with a different value for β (Cannizzo et al. 1990; Pringle 1991). Under these models, it is reasonable to approximate the fuel supply as a deterministic process, $\dot{M}(t) \propto t^{-\beta }$. If the evolution of the fuel supply is a deterministic and one-to-one process, and assuming that p(t|M_BH) ∝ 1/t_BL, then $p(\dot{M}|M_{\rm BH})$ is

$$\begin{equation} p(\dot{M}|M_{\rm BH}) = \frac{1}{t_{\rm BL}} \frac{dt(\dot{M},M_{\rm BH})}{d\dot{M}}, \end{equation} \tag{ 22 }$$

where $t(\dot{M},M_{\rm BH})$ is the time implied by $\dot{M}(t,M_{\rm BH})$, i.e., $t(\dot{M},M_{\rm BH})$ is the inverse function of $\dot{M}(t,M_{\rm BH})$. If the quasar light curve is a many-to-one function, then the right side of Equation (22) is replaced by a sum of derivatives, as in Equation (15) of Yu & Lu (2004). Inserting Equation (22) into Equation (20) and dropping the time integral gives¹²

$$\begin{equation} p(L|M_{\rm BH}) = \frac{1}{t_{\rm BL}} \int _{0}^{\infty } p(L|\dot{M},M_{\rm BH}) \frac{dt(\dot{M},M_{\rm BH})}{d\dot{M}}\ d\dot{M}. \end{equation} \tag{ 23 }$$

For models where $\dot{M}(t) \propto t^{-\beta }$, a power-law distribution of fueling rates follows from Equation (22), $p(\dot{M}|M_{\rm BH}) \propto \dot{M}^{-(1 + 1/\beta)}$ (see also Equation (43) in Yu & Lu 2008). Therefore, assuming $\dot{M}(t) \propto t^{-\beta }$ gives

$$\begin{equation} p(L|M_{\rm BH}) \propto \int _{0}^{\infty } p(L|\dot{M},M_{\rm BH}) \dot{M}^{-(1 + 1/\beta)}\ d\dot{M}. \end{equation} \tag{ 24 }$$

The feedback-driven model of Hopkins & Hernquist (2006) predicts a value of β ∼ 2, while the disk evolution model predicts β ∼ 1.2 (Cannizzo et al. 1990; Yu et al. 2005; King & Pringle 2007), implying that $p(\dot{M}|M_{\rm BH}) \propto \dot{M}^{-1.5}$ or $p(\dot{M}|M_{\rm BH}) \propto \dot{M}^{-1.67}$, respectively.

5.1.2. Effects of Time-dependent Accretion Disks

Our results in this work are broadly in agreement with recent observational and theoretical work suggesting that SMBH growth and spheroid formation have a common origin. In particular, models where mergers dominate black hole growth predict that major mergers of gas-rich galaxies initiate a burst of star formation (Mihos & Hernquist 1994a, 1996), during which the black hole undergoes Eddington-limited obscured growth. Eventually, the black hole becomes massive enough for radiation-driven feedback to unbind the surrounding gas, halting the accretion flow and revealing the object as a BLQSO (e.g., Hopkins et al. 2006a). As the activity further declines, the remnant will redden and become quiescent, satisfying the black hole–host galaxy correlation and leaving a dense stellar remnant from the starburst (Mihos & Hernquist 1994b). This dense stellar remnant is identifiable as a second component in the light profiles of elliptical galaxies (Hopkins et al. 2008a, 2009a, 2009c). Techniques based on the argument of Soltan (1982) have concluded that the SMBH accretion rate density of the universe peaks at z ∼ 2 (e.g., Marconi et al. 2004; Merloni & Heinz 2008; Shankar et al. 2009), in agreement with the observed peak in the SMBH luminosity density of the universe at z ∼ 2 (e.g., Wolf et al. 2003; Hopkins et al. 2007c). We have found that the peak in the BLQSO cosmic mass density also occurs at z ∼ 2, in broad agreement with models where the SMBH undergoes Eddington-limited growth up to the BLQSO phase, at least on a cosmological level.

In addition, our results suggest that a lower bound on the typical length of time for which a massive (M_BH ∼ 10⁹ M_☉) SMBH could be seen as a BLQSO during a single BLQSO episode is t_BL ≳ 120 Myr, i.e., a few Salpeter times. If all M_BH ∼ 10⁹ M_☉ SMBHs go through a BLQSO phase, and if the number density of SMBHs with M_BH ∼ 10⁹ M_☉ does not significantly increase from z  = 1 to z  = 0, then Equation (14) is no longer a lower bound and t_BL ∼ 150 Myr represents an estimate of the lifetime of a single BLQSO phase. Indeed, if the SMBH's growth is self-regulated then there is good reason to assume that most SMBHs go through a BLQSO phase at the end of their growth (Hopkins & Hernquist 2006), as the obscured/unobscured dichotomy represents an evolutionary sequence, at least at z > 1. Furthermore, recent work employing the argument of Soltan (1982) suggests that the number density of M_BH ∼ 10⁹ M_☉ SMBHs does not significantly increase from z  = 1 to z  = 0 (Merloni & Heinz 2008). These two considerations imply that t_BL ∼ 150 Myr should represent a reasonable estimate of the length of the BLQSO phase.

Our estimated BLQSO lifetime of t_BL ∼ 150 Myr is longer than the value of t_BL ∼ 10–20 Myr predicted by Hopkins et al. (2006a). Part of this discrepancy may be due to the fact that they defined a BLQSO as having a B-band luminosity brighter than some fraction of the host galaxy's, while we define a BLQSO as simply having broad emission lines along our line of sight. Therefore, if a BLQSO has a decaying light curve, as assumed in the model of Hopkins et al. (2006a), then by our definition an SMBH may still be considered a BLQSO even after the B-band luminosity has fallen below its host's. This would imply that the predicted BLQSO lifetime under the definition of Hopkins et al. (2006a) would be shorter than our estimated value, as is the case. However, it is unclear if the lifetime of the BLQSO should be a factor of ∼10 shorter under the definition of Hopkins et al. (2006a). Unfortunately, we cannot invert our estimated Eddington ratio distribution to a quasar light curve to test this, as the inversion is not unique, and we would have to assume a distribution of host galaxy B-band luminosity during the BLQSO phase at z ∼ 1, which is poorly constrained; as Hopkins et al. (2006a) point out, uncertainty in the ratio of host galaxy to quasar B-band luminosity also introduces considerable systematic uncertainty in their predicted BLQSO lifetime. Moreover, our estimated lifetime is subject to the assumptions outlined in Section 4.3, which may introduce additional systematic uncertainties of a factor of a few. Considering this, a more quantitative comparison between our estimate and the prediction from Hopkins et al. (2006a) is difficult.

Our estimated BLQSO BHMF suggests a typical value of M_BH ∼ 10⁸ M_☉ for BLQSOs. Theoretical estimates of the mass distribution of SMBH seeds suggest that typical values for the seed mass should be M_BH ∼ 10⁵–10⁶ M_☉ at z ∼ 10 (Lodato & Natarajan 2007; Pelupessy et al. 2007; Volonteri et al. 2008), and grow to M_BH ∼ 10⁸ M_☉ by z ∼ 3 (Volonteri & Natarajan 2009), consistent with our results. If the black hole's growth is Eddington limited, then at least several Salpeter times are required to grow a seed black hole to M_BH ∼ 10⁹ M_☉ from an M_BH ∼ 10⁶ M_☉ seed, which is inconsistent with our estimated BLQSO lifetime. Furthermore, we find that most SMBHs in BLQSOs are not radiating near the Eddington limit, with a typical value being L/L_Edd ∼ 0.1 at M_BH ∼ 10⁹ M_☉, as we might expect if BLQSOs are seen during a phase with a decaying fueling rate. Because most BLQSOs are radiating at considerably less than the Eddington limit, this therefore suggests that a factor of ∼10 longer is needed to grow these SMBHs from seeds of M_BH ∼ 10⁶ M_☉, i.e., a growth timescale of ∼70 Salpeter times. This corresponds to the age of the universe at z ∼ 2, implying that sources observed at z > 2 had to have been accreting at higher rates earlier in their growth.

An inferred growth timescale of ∼70 Salpeter times is significantly longer than our estimated BLQSO lifetime at z ∼ 1, implying one of a few explanations. First, it implies that if we assume that BLQSO SMBHs spend all of their growth as a BLQSO, then our assumption that all SMBHs go through a BLQSO phase along our line of sight is incorrect. In other words, this implies that some of the SMBH population that is growing at any redshift could never be observed by us to be a BLQSO at any time. Because the lifetime is related to what we can calculate from the BHMF by Equation (14), and if the SMBH spends all of its growth as a BLQSO, then an inferred SMBH growth timescale of t_BL ∼ 70t_salp implies that only ∼4% of SMBHs ever go through a phase where we could observe them as BLQSOs at any time. This is an unrealistically low number, and thus we conclude that t_BL is shorter than the growth time. Furthermore, it also illustrates that even if we can never observe a large fraction of active SMBHs, say ∼50%–75%, possibly due to orientation-dependent obscuration, then our estimated t_BL is only underestimated by a factor of ∼2–4.

Another possibility is that SMBHs in BLQSOs at z ∼ 1 built up their mass via numerous fueling events of length t_BL ≳ 150 Myr. In this case, one could grow an SMBH from a seed mass of M_BH ∼ 10⁶ M_☉ at z ∼ 10, say, to a mass of M_BH ∼ 10⁹ M_☉ by z ∼ 1, without the need for an earlier phase of obscured and accelerated growth. Similarly, we may have incorrectly assumed that the lifetime of the BLQSO phase is short compared to the timescale needed for any significant change in the BLQSO triggering time distribution, as this assumption was needed in order to use the observed distribution of BLQSOs as an estimate of their triggering distribution. This assumption may be incorrect if BLQSOs undergo a single long growth phase from z ∼ 10 to z ∼ 1. However, this possibility only exists for z ≲ 1.7, since at higher redshifts the growth time for objects radiating at L/L_Edd ∼ 0.1 is longer than the lookback time to z ∼ 10, and thus BLQSOs that are observed at z ≳ 1.7 would have had to experience an earlier phase of obscured growth at an enhanced accretion rate. Moreover, if SMBHs that are seen as BLQSOs at z ∼ 1 with M_BH ∼ 10⁹ M_☉ are grown in a single long BLQSO episode, or numerous repeated ones of t_BL ∼ 150 Myr, this would also imply that BLQSOs that are seen at z ≲ 1.7 would have had a different fueling mechanism than BLQSOs that are seen at z ≳ 1.7, and it is unclear why this should be true. Therefore, we conclude that while part of the mass of SMBHs in BLQSOs may have been accumulated via multiple BLQSO episodes, it is unlikely that most of these SMBHs did not also experience a phase of obscured growth at an enhanced accretion rate.

If most SMBHs go through a BLQSO phase along our line of sight at some point in their growth, the fact that our estimated BLQSO lifetime is short compared to the SMBH growth time implies that SMBHs spend a significant amount of their time growing in a non-BLQSO phase, such as an obscured phase. Note that this argument is unaffected by the fact that at any z we miss a considerable fraction of the SMBH population that is growing (e.g., due to obscuration), so long as these missed SMBHs undergo a BLQSO episode at some point in their life. The conclusion that a significant amount of growth occurs in an earlier obscured phase was also reached by Treister et al. (2010), and is expected from self-regulation models where mergers or other triggering mechanisms fuel and initiate Eddington-limited accretion and obscured SMBH growth, until the SMBH becomes massive enough to unbind the ambient gas, revealing it as a BLQSO for t_BL ∼ 150 Myr. Within this interpretation, the BLQSO black hole mass density shown in Figure 6 is proportional to the rate at which the relic SMBH mass increases as a function of redshift. Furthermore, considering that the growth timescale for an SMBH radiating at L/L_Edd ∼ 0.1 to grow from 10⁶ M_☉ to 10⁹ M_☉ corresponds to the age of the universe at z ∼ 2, we also conclude that SMBHs accrete at a significantly higher rate during the earlier obscured phase, as compared to the BLQSO phase.

In this work, we have found that the most massive SMBH that could be seen as a BLQSO is M^Max_BH ∼ 3 × 10¹⁰ M_☉ and would most likely be observed at z ≳ 2. These constraints on M^Max_BH are consistent with recent cosmological simulations of SMBH growth, as well as expectations from black hole feedback models. Cosmological simulations that follow the growth of SMBHs in bright z ∼ 6 quasars have been able to grow SMBHs to M_BH ∼ 10⁹ M_☉ by z ∼ 4 (Li et al. 2007; Di Matteo et al. 2008) and M_BH ∼ 2 × 10¹⁰ M_☉ by z ∼ 2 (Sijacki et al. 2009). Similarly, considerations based on quasar feedback and self-regulated SMBH growth, combined with the local distribution of bulge velocity dispersion, also suggest a value of M^Max_BH ∼ 10¹⁰ M_☉ (Natarajan & Treister 2009).

If the SMBHs growth is self-regulated, then the final mass of the SMBH after a fueling event is set by the binding energy of the bulge regardless of the fueling mechanism (Younger et al. 2008). In order to build up a mass as large as M^Max_BH, multiple fueling episodes would likely be necessary, as is seen in cosmological simulations (e.g., Li et al. 2007; Di Matteo et al. 2008; Sijacki et al. 2009). However, it is unlikely that M^Max_BH could be increased significantly beyond the observed value as additional fueling mechanisms would likely result in only a small increase in the binding energy of the bulge. Therefore, additional fueling episodes would immediately result in the SMBH unbinding the accreting gas, preventing significant additional growth. Indeed, if our estimated BLQSO lifetime of t_BL ∼ 150 Myr is correct for SMBHs with M_BH ∼ 10⁹ M_☉, and if the SMBH is radiating at a typical Eddington ratio of L/L_Edd ∼ 0.1 during the BLQSO phase, then it would only accrete an additional ∼10⁸ M_☉, assuming a radiative efficiency of _r = 0.1. If the distribution of Eddington ratios is a power law with p(Γ_Edd) ∝ Γ^−1.5_Edd, as suggested by the discussion in Section 5.1.1, then the typical Eddington ratio would be L/L_Edd ≲ 0.1, suggesting even less growth during the BLQSO phase. Moreover, additional growth via black hole mergers, as might be expected from "dry" mergers of galaxies, is also unlikely to lead to significant additional growth, as the mass of the additional SMBH is likely to be negligible compared to M^Max_BH. And finally, additional growth through radiatively inefficient low accretion rate modes does not contribute significantly to the black hole's final mass (Hopkins et al. 2006c; Cao 2007). Therefore, our estimated value of M^Max_BH ∼ 3 × 10¹⁰ M_☉ should be representative of the most massive SMBH for both active and inactive SMBHs.

Our results are consistent with previous data-based work that has attempted to map the distribution and growth of SMBHs. Most previous observational work has mapped SMBH growth by using the M_BH–σ relationship to infer the local BHMF, and then stepped backward in time using the quasar luminosity function to infer the contribution to the BHMF from accretion (e.g., Soltan 1982; Yu & Tremaine 2002; Shankar et al. 2004; Marconi et al. 2004; Merloni & Heinz 2008). In addition, there have been attempts to predict the BHMF and its evolution for all SMBHs, or all active SMBHs (e.g., Volonteri et al. 2003; Sijacki et al. 2007, 2009; Hopkins et al. 2008b; Di Matteo et al. 2008; Shen 2009). While not directly comparable to these studies, as we focus on BLQSOs, our results and conclusions are qualitatively consistent with previous observational and theoretical work in that we find evidence for self-regulated SMBH growth, black hole downsizing, and BLQSO lifetimes of t_BL ∼ 150 Myr.