Abstract
The black hole mass function of supermassive black holes describes the evolution of the distribution of black hole mass. It is one of the primary empirical tools available for mapping the growth of supermassive black holes and for constraining theoretical models of their evolution. In this paper, we discuss methods for estimating the black hole mass function, including their advantages and disadvantages. We also review the results of using these methods for estimating the mass function of both active and inactive black holes. In addition, we review current theoretical models for the growth of supermassive black holes that predict the black hole mass function. We conclude with a discussion of directions for future research which will lead to improvement in both empirical and theoretical determinations of the mass function of supermassive black holes.
1. Introduction
Understanding how and when supermassive black holes (SMBHs) grow is currently of central importance in extragalactic astronomy. A significant amount of empirical work has established correlations between SMBH mass and host galaxy spheroidal properties, such as luminosity [1–3], stellar velocity dispersion (the - relationship, e.g., [4–6]), concentration or Sersic index [7, 8], bulge mass [9–11], and binding energy [12, 13]. These scaling relationships imply that the evolution of spheroidal galaxies and the growth of SMBHs are intricately tied together. The currently favored mechanism for linking the growth of SMBHs and their hosts is black hole feedback, whereby black holes grow by accreting gas in the so-called “active’’ phases, possibly fueled by a major merger of two gas-rich galaxies, until feedback energy from the SMBH expels gas and shuts off the accretion process [14–17]. Alternatively, it has been suggested that the origin of the scaling relationships does not necessarily require SMBH feedback but emerges from the stochastic nature of the hierarchical assembly of black hole and stellar mass through galaxy mergers [18, 19].
Feedback-driven “self-regulated” growth of black holes has been able to reproduce the local - relationship in smoothed particle hydrodynamics simulations [20–22]. Moreover, AGN feedback has also been invoked as a means of quenching the growth of the most massive galaxies [23, 24]. There have been numerous models linking SMBH growth, the quasar phase, and galaxy evolution [25–33]. While feedback is likely important for regulating the growth of SMBHs and galaxies, the fueling mechanisms that contribute to growing the SMBH are likely diverse. Major mergers of gas-rich galaxies may fuel quasars at high redshift and grow the most massive SMBHs. However, major mergers alone do not appear to be sufficient to reproduce the number of X-ray faint AGN [34], and accretion of ambient gas via internal galactic processes [35, 36] may fuel these fainter, lower AGN at lower . This is supported by the fact that many AGN are observed to live in late-type galaxies out to [37, 38], and the X-ray luminosity function of AGN hosted by late-type galaxies suggests that fueling by minor interactions or internal instabilities represents a nonnegligible contribution to the accretion history of the universe [39].
The black hole mass function (BHMF) provides a complete census of the mass of SMBHs and their evolution. Because of this, the BHMF is one of the primary empirical tools available for investigating the growth of SMBHs, and for constraining theoretical models for the growth of the SMBH population. Because SMBHs and galaxies are thought to be linked in their evolution, the BHMF provides insight into the fueling mechanisms that dominate black hole growth and therefore into the role of feedback in the evolution of the host galaxy. The BHMF is also an important tool in planning future surveys, as it provides an estimate of the distribution of SMBH mass expected for the survey. This in turn is important because mass is a fundamental quantity of the black hole and therefore is an important observational quantity for empirical studies of black hole accretion physics [40–45]. Of course, further improvement to our understanding of black hole accretion physics will further improve our modeling and understanding of black hole accretion and feedback, which in turn will improve our understanding of black-hole-galaxy coevolution. Therefore, the BHMF is an important empirical quantity for SMBH studies.
In this paper, we discuss the current status of BHMF estimation and theoretical modeling. In Section 2, we discuss the nontrivial task of estimating the BHMF. In Section 3, we discuss current estimates of the local SMBH BHMF. In Section 4, we discuss BHMF estimates derived by combining the local BHMF with the AGN luminosity function via a continuity equation. In Section 5, we discuss BHMFs estimated for AGN only. In Section 6, we review theoretical models for SMBH growth that predict the SMBH BHMF. Finally, in Section 7, we discuss directions for future improvements to the empirical and theoretical studies of the BHMF. We note that unlike, say, the luminosity function, the division between “observational” and “theoretical” studies is not as clear for the BHMF, as some amount of modeling is necessary in order to estimate the BHMF from strictly observational quantities. We have attempted to divide the studies according to whether the BHMF is constrained empirically, as in, say, a formal statistical fitting procedure, or if it is predicted from a theoretical model for SMBH growth. In reality, the line between theoretical and empirical studies is blurry and some procedures which we have considered to be empirical may be thought of as theoretical.
2. Estimating the Black Hole Mass Function
The black hole mass function, denoted as , is the number of sources per comoving volume with black hole masses in the range , . The black hole mass function is related to the joint probability distribution of and , , as The normalization of the BHMF is , the total number of SMBHs in the observable universe, and is given by the integral of over and .
2.1. Complications with Estimating Black Hole Mass Functions
Similar to luminosity function estimation, the BHMF may be estimated from astronomical surveys. However, while there are many well-established methods for estimating luminosity functions, there are two complications that make BHMF estimation a more difficult problem [46]. The first is the issue of incompleteness. Surveys are typically constructed by finding the set of objects of interest containing a SMBH that satisfy a flux criteria, for example, all objects brighter than some flux limit. Surveys are not constructed by selecting on mass. Because there is a distribution of luminosities at a given SMBH mass, whether it is the luminosity of the host galaxy or of the AGN, some SMBHs will scatter above the flux limit and some below. This creates a selection function which is less sensitive to , and it is possible that a survey may be incomplete in all mass bins.
The second complication is the large uncertainty in SMBH mass among mass estimators. Currently, it is not possible to obtain reliable mass estimators for large numbers of SMBHs through dynamical and modeling of the stellar or gaseous components, and thus scaling relationships are employed. Masses may be estimated using scaling relationships between and the properties of the host galaxy bulge or the luminosity and the width of the broad emission lines for AGN [47, 48]. It has also recently been suggested that the X-ray variability properties of AGN may also provide another scaling relationship for estimating [45, 49, 50], but further work is needed for developing this. While these scaling relationships enable one to estimate for large numbers of SMBHs, they also contain a significant intrinsic statistical scatter. Gültekin et al. [51] find that for early-type galaxies there is an intrinsic scatter in of dex and dex at fixed host galaxy bulge dispersion and luminosity, respectively; the amplitude of the scatter is larger for late-type galaxies. For AGN with broad emission lines, Vestergaard and Peterson [48] estimate the scatter in at fixed luminosity and line width to be ~0.4 dex, depending on which emission line is used.
The statistical uncertainty in the mass estimates can have a significant effect on the inferred BHMF. The distribution of the mass estimates is the convolution of the intrinsic BHMF with the error distribution in the mass estimates. In general, it is typically assumed that the error in the mass estimates is independent of the actual value of . This is not the case for estimated through dynamical modeling; however, independence between and its error is likely to be a good approximation for estimated using scaling relationships. Because scaling relationships are the only feasible manner to estimate for a large sample of SMBHs, which is necessary for any estimate of the BHMF, we will assume that and its error are independent. Under the assumption of independence between the estimated and its error, the BHMF that would be inferred directly from the distribution of the mass estimates is broader than the intrinsic BHMF and is thus biased. Figure 1 illustrates this effect, where an intrinsic mass function is compared with the distribution of an unbiased mass estimator having a statistical uncertainty of 0.3, 0.4, and 0.5 dex, respectively. As can be seen, the distribution of mass estimates is significantly different from the intrinsic mass function. In particular, the distribution of the mass estimates falls off more slowly with increasing , and overpredicts the number of SMBHs at the high end of the mass function. The bias is worse when the dispersion in the scatter in the mass estimates becomes larger.
2.2. Methodology for Estimating the Black Hole Mass Function
In order to estimate the SMBH mass function in an unbiased manner, it is necessary to match the mass function with the observed distribution of the mass estimates and any additional observational quantities that the selection function (The selection function is the probability of including a source in one’s sample as a function of its measured quantities) depends on. The basic idea is to start with an assumed mass function. Then, calculate the distribution of mass estimates implied by this mass function. In addition, calculate the distribution of observational quantities that one’s sample is selected on, say, flux, that is implied by the assumed mass function. This step allows one to correct for incompleteness but requires an additional assumption about how to relate the mass function to the quantity that one’s sample is selected on. Finally, impose the selection function for the sample and compare the predicted observed distributions of mass estimates and any other observables (e.g., flux) with the actual distributions. If they are not consistent, then the data rule out the assumed mass function and relationship between and the observable quantities.
We can make the above procedure more quantitative by deriving the likelihood function for the SMBH mass function. Kelly et al. [46] derived the likelihood function for the mass function when using masses estimated from AGN broad emission lines. They used this likelihood function for developing a Bayesian approach to estimating the SMBH mass function. Although their method was limited to broad-line mass estimates, it is straightforward to generalize their formalism for any generic mass estimator. Denote the black hole mass estimate as . In addition, denote as the set of observables that one uses to select one’s sample. In the majority of cases, this will be flux at one or more wavelengths. Then, the likelihood function for the BHMF based on a sample of SMBHs is where the BHMF is related to and the probability distribution of and via (1). Here, is the binomial coefficient, denotes the parameters for the BHMF, denotes the parameters for the distribution in at fixed and , and is the probability of including a SMBH in one’s sample as a function of and . Here, we have assumed that the distribution in the mass estimates at fixed , , and , is known, although one could include additional free parameters for this as well. The binomial coefficient arises from the fact that the number of objects included in one’s survey follows a binomial distribution (Often in the luminosity function literature the likelihood is assumed to be a Poisson distribution. A poisson distribution is an approximation to the binomial distribution when and , so (2) converges to the Poisson distribution when . See [52] for further details.) with “trials” and probability of success . The probability of including a SMBH in one’s survey as a function of the BHMF, , is calculated from the survey selection function as It is up to the researcher to choose the particular parametric form for the SMBH mass function, the distribution in the mass estimates at fixed , , and , and the distribution of the observable that the sample is selected on (e.g., flux) at fixed and . Typical choices are log-normal distributions, Schechter functions, and mixtures of log-normal distributions. Once one has done this, one can use (2) to compute a maximum-likelihood estimate for the BHMF or perform Bayesian inference.
An alternative form of estimating the BHMF can be used when the mass estimates are derived from an observational quantity, , and the intrinsic distribution of is known. This is commonly used to estimate the local mass function using host-galaxy scaling relationships [53]. In this case, the mass function is where is the comoving number density of SMBHs as a function of the quantity . When both and are known, then the BHMF follows directly from (4). As an example, if the mass function is derived from the scaling between and host galaxy spheroidal luminosity, , then , is the - relationship and is the luminosity function of stellar bulges hosting SMBHs. As with BHMFs determined from a mass estimator, improper treatment of the intrinsic scatter in at fixed will lead to a biased estimate of the BHMF. However, when calculating the BHMF from (4), ignoring the intrinsic scatter results in an estimated BHMF that is too narrow, underpredicting the number of SMBHs at the high-mass end of . This is opposite to the case when one estimates the BHMF directly from mass estimates.
3. Black Hole Mass Functions Derived from Host Galaxy Scaling Relationships
The observed scaling between and the properties of the SMBH host galaxy bulge have motivated several groups to estimate the local BHMF [53–65], with decreasing statistical uncertainties. These estimates of the local BHMF have formed the basis for many studies which have attempted to map black hole growth by comparing with the AGN luminosity function, this is further discussed in Section 4. Typically, the local BHMF is estimated using the local - relationship or the local - relationship, combined with the local number density of galaxies as a function of stellar velocity dispersion or bulge luminosity.
The scaling relationships between and host galaxy properties are only determined for the local universe, and thus most authors have limited their determination of the BHMF based on them to the local BHMF. There are, however, a couple of exceptions. Tamura et al. [66], estimated the BHMF out to assuming that evolution in the - relation is driven only by passive evolution in . Shankar et al. [67] (see also [68, 69]) used the local velocity dispersion function of spheroids in combination with their inferred age distributions to estimate the BHMF at . In order to do this, they assumed that most of the stars in nearby spheroids formed in a single event and that did not change once the spheroid was formed. In addition, Shankar et al. [67] allowed the normalization of the - relationship to evolve, with the degree of evolution being a free parameter. They found evidence for mild evolution in the normalization of the - relationship.
Evolution in the scaling relationships is currently an area of intense study, with most groups finding evidence that the normalization of the scaling relationships increases towards higher [70–75], at least for active SMBHs. However, there are still concerns regarding potential biases due to selection effects [76], but see Treu et al. [70] and Bennert et al. [75] for procedures aimed at modeling and correcting for selection. There may also be biases due to extrapolating the AGN mass estimates derived from the broad emission lines to luminous quasars at high [77]. As such, the uncertainties on the quantitative form of the evolution in the scaling relationships and their scatter are currently large, limiting their use for determining the BHMF outside of the local universe.
When the - relationship is used to estimate the local BHMF, it is common to use the velocity dispersion distribution derived from the SDSS by Sheth et al. [78], with an additional component representing the brightest cluster galaxies [58]. Sheth et al. [78] estimate the velocity dispersion distribution for late-type galaxies by using the Tully-Fisher relation to convert the luminosity function of late-type galaxies to a circular velocity distribution and then set . When the relation is used, it is typical to estimate the distribution of separately for early- and late-type galaxies by converting their respective luminosity functions to spheroidal luminosity functions using an assumed ratio of bulge luminosity to total luminosity. From this, it has been inferred that the local BHMF is dominated by early-type galaxies at [62]. Shankar et al. [64] present a compilation of recently determined local BHMFs based on a variety of methods, scaling relations used, and data sets used. In Figure 2, we show the range of local BHMFs estimated from the -, -, and - relationships, as presented in Shankar et al. [64]. In general, estimates of the local mass density of SMBHs span the range for [64].
While the procedure for estimating the local BHMF is, in theory, straightforward, a number of significant systematics remain. First, there is the observational difficulty that most BHMFs derived from the - relationship are based on SDSS spectra. Unfortunately, the SDSS velocity dispersions are based on a fixed aperture, and thus the size of the aperture relative to the bulge varies with the apparent size of the galaxy and its inclination. In addition, the spectral resolution of SDSS spectra is ~100 kms−1, making it difficult to reliably measure for SMBHs with . Another concern is that the local BHMF is derived by assuming that the - or - relations are single power laws with a constant scatter in at fixed or . However, recent work has shown these assumptions to be incorrect. For one, the - and - relations diverge at the high- end, which Lauer et al. [58] suggest implies that the - relation is not a single power law. This divergence creates an inconsistency in the BHMFs derived from these two scaling relationships [58, 61]. Similarly, the - relationships for the SDSS and dynamical SMBH samples are inconsistent, suggesting a possible selection bias in the estimated BHMFs [55, 79]. The scatter in the - relation is larger for spirals [51, 80], and appears to increase at low such that most SMBHs lie below the - relation, see (e.g., [81]). Several authors have found differences in the slope and scatter of the scaling relations for pseudobulges [80, 82–84]; however, it is unclear that this result is due to differences in the perceived bulge velocity dispersions for bulges as compared to pseudobulges or due to different scaling relationships. Recently, Kormendy et al. [85] argue that does not correlate with galaxy disks and only correlates weakly, if at all, with pseudobulges. On other hand, Graham et al. [86] analyzed a larger sample of barred galaxies and concluded that does correlate with , even though the - relationship for barred galaxies is offset from that of non barred galaxies. Although there is still much that we do not understand about the and host galaxy scaling relationship, these recent results suggest that the scaling relationships are not a single power law with constant intrinsic dispersion in , representing a significant source of systematic uncertainty in the estimated local BHMF, especially at the low-mass end.
4. Black Hole Mass Functions Derived from the Local Mass Function and the AGN Luminosity Function
By employing the argument of Soltan [87], numerous studies have attempted to estimate the BHMF at a variety of redshifts by comparing the accreted mass distribution implied by the quasar luminosity function with the local BHMF [53, 55, 63, 88–92]. These methods employ a continuity equation describing the evolution of the number density of SMBHs [93, 94]: Here, is the average growth rate of SMBHs as a function of and cosmic age, , and is the rate at which the number density of SMBHs changes due to mergers of black holes, or ejections of black holes from their host galaxies due to gravitational recoil. Technically, can also include a contribution from SMBHs which are created, but this has not been thought to occur over the redshift range in which (5) is typically applied, that is, . Because the merger rate of black holes is currently unknown, many studies that have employed (5) set .
Under the assumption that SMBHs grow during phases of AGN activity, AGN demographics in combination with the local BHMF may be used to compute . This is because the AGN luminosity function maps the accretion history onto SMBHs, and the local BHMF acts as a boundary condition on (5); it is also possible in principle to include the BHMF for AGN, which provides more information. Studies that have used (5) to estimate the BHMF generally fall into two categories: those that assume an AGN lightcurve and those that employ the BHMF of AGN. We discuss each of these separately.
4.1. Methods That Assume an AGN Lightcurve
Most authors employing (5) have assumed a parametric form for . The accretion rate is related to the bolometric luminosity output of the accretion flow onto the SMBH as , where is the radiative efficiency of the accretion flow, is the accretion rate of matter onto the SMBH, and is the speed of light. The growth rate of the SMBH is , due to the fact that a fraction of accreted mass is radiated away as energy. Making this substitution, the continuity equation becomes where we have ignored mergers of SMBHs. Equation (6) shows that it is possible to calculated the BHMF at a time given the local BHMF, an assumed average accretion flow lightcurve as a function of , , and an assumed radiative efficiency. Because and imply a luminosity function, the local BHMF and AGN luminosity function can be used to place constraints on and . This means that, in practice, one also has to assume a bolometric correction, which itself likely depends on both black hole mass [95] and [96, 97]. In addition, an estimate of also enables one to estimate the lifetime and duty cycle of AGN activity, modulo some luminosity-dependent definition of an AGN; note that the AGN duty cycle defines the fraction of SMBHs that are “active” at a given and .
A variety of lightcurve models have been used when employing (6) to reconstruct the evolution of the BHMF. The simplest model is that where SMBHs spend a fraction of their time radiating at a constant Eddington ratio and spend the remainder of their time in quiescence. The free parameters in this model are the Eddington ratio, AGN lifetime or duty cycle, and radiative efficiency. This model has been used by [53–55, 57, 64] (technically, [54] assumed that the Eddington ratio was a weakly increasing function of luminosity) to study the build-up of the local black hole mass function, although [64] also considered models where the average accretion rate relative to Eddington falls off toward lower and higher . Raimundo and Fabian [98] employed a variation on the constant models, assuming three different populations of AGN with their own Eddington ratio: a population of obscured low AGN, a population of obscured AGN with higher , and a population of unobscured AGN. Yu and Lu [62] modeled the quasar lightcurve as radiating at the Eddington limit for a period of time, and then transitioning into a power-law decay. Cao [92] also modeled the quasar lightcurve as undergoing a power-law decay. Lightcurves undergoing a power-law decay arise from self-regulation models and describe the evolution of the lightcurve after black hole feedback unbinds the accreting gas, therefore quenching its fuel supply. The power-law decay occurs either as a result of evolution of a blast wave [36, 99] or from viscous evolution of the accretion disk [100, 101].
In Figure 3, we compare the BHMF calculated by Shankar et al.[64] with that calculated by Cao [92]. For Shankar et al. [64], we show their reference model, which assumes a radiative efficiency of , an accretion rate relative to Eddington of , and that half of all SMBHs are active at . We show the model from Cao [92] which assumes a radiative efficiency of and a quasar lifetime of yr, as it better matches the Shankar et al. [64] estimates. The two estimates of the BHMF agree fairly well, despite the different quasar lightcurve models.
(a)
(b)
(c)
(d)
In general, most of the studies that have used (6) in combination with an assumed quasar lightcurve have concluded the following.(i)Most SMBH growth occurs in periods when the quasar is radiating near the Eddington limit.(ii)Most, if not all, of the local black hole mass function can be explained as the relic of previous AGN activity, implying that mergers of SMBHs are not important for building up the local mass function.(iii)SMBH growth is antihierarchical, with the most massive black holes growing first. This has also been termed “downsizing” of active SMBHs.(iv)The lifetime of AGN activity is ~ a few ×108 yr.(v)Most SMBHs have nonzero spin, as implied by inferred radiative efficiencies of .
However, while (6) has proven to be an important tool for studying SMBH growth and estimating the black hole mass function, it must be kept in mind that the use of (6) often entails some strong assumptions. These methods rely on the assumed form of the quasar lightcurve, distribution of radiative efficiencies, and bolometric corrections, all of which are subject to considerable uncertainty. Moreover, in general, these methods also rely on an estimate of the local black hole mass function, which, as discussed in Section 3, is itself subject to considerable uncertainty. Indeed, there is a strong degeneracy between the estimated radiative efficiency of accretion and the normalization of the local BHMF, and therefore the uncertainty in is linearly proportional to that in the normalization (or integral) of the local BHMF. All of these issues have the potential to introduce systematic error into methods based on (6), and further work is needed in reducing these systematics.
4.2. Methods That Include the Distribution of Active Supermassive Black Holes
An alternative to the methods described in Section 4.1 is to estimate the average value of the accretion rate onto SMBHs directly from the observational data. This avoids the issue of assuming a form for the quasar lightcurve, as instead is derived directly from the estimated distribution of . Techniques based on this approach require a means of linking the mass function of active SMBHs to observational quantities, which is done via scaling relationships. This was the approach of Merloni [89] and Merloni and Heinz [63], who employed the black hole “fundamental plane” (BHFP) [40, 41].
The BHFP is a scaling relationship between , radio luminosity, and X-ray luminosity, that exists for low-accretion rate black holes (i.e., ), extending from galactic black holes to supermassive ones. The BHFP likely reflects the connection between and the conversion of the accretion flow into radiative energy and jet power. It, in principle, enables one to connect the radio and X-ray luminosity functions to the mass function of active SMBHs. Having obtained a distribution of and X-ray luminosity at a given redshift for the active SMBH population, Merloni [89] and Merloni and Heinz [63] then convert this to a joint distribution of and assuming a conversion from X-ray luminosity to which depends on the Eddington ratio. The joint distribution of and at a given redshift for active SMBHs therefore enables calculation of the average growth rate , which can then be combined with the continuity equation to calculate the black hole mass function at the next redshift. Their estimated BHMF is shown in Figure 4, which is a recreation of their Figure 5. Similar to methods based on assuming a quasar lightcurve, Merloni and Heinz [63] concluded that SMBHs grow antihierarchically; however, in contrast to the lightcurve methods, Merloni and Heinz [63] concluded that most SMBHs have low spin as inferred from their derived radiative efficiency. In addition, Merloni and Heinz [63] concluded that the distribution of SMBH accretion rates is broad and that most SMBH growth occurs during a radiatively efficient accretion mode.
The method of estimating the BHMF from the BHFP developed by Merloni [89] and Merloni and Heinz [63] has the advantage that it derives the distribution of accretion rates empirically. However, there are also disadvantages to this approach. The uncertainties regarding the bolometeric correction, estimation of the local BHMF, and radiative efficiency also apply to the BHFP method as well. Moreover, as discussed in Merloni and Heinz [63], the BHFP is only defined for low-accretion rate objects, that is, objects with . Merloni and Heinz [63] extrapolate the BHFP to higher accretion rates, after rescaling the normalization to ensure that the radio luminosity is weak for AGN in the radiatively efficient mode (those objects with and lacking a jet). Unfortunately, the AGN in the radiatively efficient mode make a significant contribution to the X-ray luminosity function, from which is derived. Moreover, most studies, including those based on the BHFP, have concluded that most SMBH growth occurs at , which corresponds to the radiatively efficient mode. Because the radiatively efficient mode also corresponds to the regime of largest systematic uncertainty for the BHFP, there is the potential for significant systematic error in estimating the BHMF based on the BHFP, as well as in estimating the primary mode of SMBH growth. There is thus a need for further improvement to our understanding of the scaling relationships involving and the AGN SED.
5. Black Hole Mass Functions of AGN
Thus far, we have focused on methods for estimating the mass function of all SMBHs. In this section, we will describe methods for estimating the BHMF for those SMBHs in AGN and the results that have come from the application of these methods.
5.1. Methods Based on Scaling Relationships Involving the Broad Emission Lines
The steady improvement in reverberation mapping of AGN [102, 103] has revealed a correlation between the luminosity of AGN and the broad-line region radius [104, 105]. It is therefore possible, in principle, to obtain an estimate of for broad-line AGN (BLAGN) by combining a luminosity-based estimate of the broad-line region size with an estimate of the velocity dispersion of the broad-line region gas obtained from the width of the broad emission lines [47]. These virial mass estimates are then calibrated to the estimates of obtained from reverberation mapping, which themselves are calibrated to be consistent with the local - relationship [106, 107]. Currently, calibrations exist for [108], (e.g., [48]), Mg-II [109–111], and C IV [48]. The statistical scatter in the virial mass estimates is currently estimated to be ~0.4 dex [48], although there are indications that the scatter may be smaller, at least for the most luminous quasars [112–116]. Moreover, it should be noted that the calibration for Mg II is obtained by enforcing consistency in the mean values of the Mg II mass estimator and the Hβ and C IV ones, and therefore there is currently no direct estimate of the statistical scatter in Mg II-based virial mass estimates. In contrast, the amplitudes of the statistical scatter for Hβ and C IV are estimated by comparing mass estimates derived from these lines with the masses derived from reverberation mapping [48]. Although there is currently very little reverberation mapping data for C IV, the estimate of the dispersion in the C IV-based mass estimates should not be biased so long as the masses based on reverberation mapping are reliable estimates of the true , regardless of which emission line was used in the reverberation mapping campaign.
Early estimates of the mass function of SMBHs in BLAGN were obtained by binning up the virial mass estimates and applying a correction [110, 117–119], a technique borrowed from luminosity function estimation. Greene and Ho [118] estimated the local BHMF for BLAGN from the SDSS DR4, while Vestergaard et al. [119] estimated the BHMF for BLAGN over using the uniformly selected quasar sample from the SDSS DR3 [120]. Vestergaard and Osmer [110] estimated the BHMF for the brightest BLAGN using objects from a variety of surveys, as their sample was designed to complement the uniformly selected SDSS DR3 sample. Unfortunately, as discussed in Section 2.1, this method of binning up the mass estimates suffers from biases due to the large statistical scatter in the virial mass estimates, and due to the inability of a luminosity-based correction to correct for incompleteness in . Subsequent attempts have further improved in their methodology, providing more accurate BHMFs.
Shen et al. [113] employed a forward-modeling approach where the mass function and Eddington ratio distribution were estimated by matching the observed distribution of mass estimates and luminosity to that implied by the model BHMF and Eddington ratio distribution. Their method accounts for incompleteness and the statistical scatter in the mass estimates but lacked statistical rigor in that the matching was done visually. Schulze and Wisotzki [121] employed a maximum-likelihood technique for estimating the local BHMF for BLAGN. Their method corrects for incompleteness in but does not correct the BHMF for the broadening caused by the statistical scatter in the virial mass estimates. Kelly et al. [46] developed a Bayesian method that corrects for both the statistical scatter in the mass estimates and incompleteness and used their method to estimate the local BHMF of BLAGN from the Bright Quasar Survey [122]. Kelly et al. [115] used the method of [46] to estimate the BHMF of BLAGN at from the mass estimates in the SDSS DR3 quasar sample [119]. The BLAGN BHMFs from a variety of studies are compiled in Figure 5, showing the evolution of the BHMF from the local universe out to . More recently, Shen and Kelly [116] extended the Bayesian method of [46] to include a possible luminosity-dependent bias in virial mass estimates derived from the emission line FWHM, the existence of which was suggested by Shen and Kelly [77]. Shen and Kelly [116] applied their method to the SDSS DR7 uniformly-selected quasar sample, independently estimating the BHMF and Eddington ratio distribution in different redshifts bins.
Similar to the methods based on the continuity equation, investigations of the BHMF for BLAGN have found evidence for the anti-hierarchical growth of SMBHs, that is, cosmic “down-sizing” of BLAGN activity. The inferred Eddington ratio distributions are wide, and the density of SMBHs continues to increase toward Eddington ratios which are below the survey completeness limit. In addition, Kelly et al. [115] used the BLAGN BHMF to estimate the lifetime of broad-line quasar activity to be Myr among SMBHs with , which is similar to quasar lifetimes inferred from the continuity equation. Kelly et al. [115] also used their estimated BHMF to estimate the maximum mass of a SMBH to be , which is in agreement with theoretical expectations [123, 124].
Mass functions estimated from scaling relationships for BLAGN have the advantage that they are derived from estimates of that are obtained for individual sources, providing a more “direct” estimate of the mass function than those based on the continuity equation. However, they have the disadvantage that they are only available for a subset of the AGN population, which itself is only a subset of the SMBH population. This complicates comparison with other SMBH mass functions, as the fraction of AGN with broad emission lines is poorly constrained, especially as a function of mass. This being said, BHMFs of BLAGN represent a subset of SMBHs that are actively growing at the time that they are observed, and, as the aforementioned studies have demonstrated, their mass function still contains important information on SMBH growth.
As with all methods of BHMF estimation, the virial mass estimates and the mass functions derived from them still suffer from systematics. First, there is the usual problem of calculating a bolometeric correction, although this only affects the estimated Eddington ratio distribution and not the BHMF. Second, there are a few concerns with the virial mass estimates which could introduce systematic error; some of these have been discussed by Greene and Ho [125]. For one, most of the reverberation mapping data is only available for the H line. Because of the limited Mg II data, the Mg II scaling relationship is in general not calibrated using objects with black hole mass estimates from reverberation mapping. There may be systematic effects with luminosity or Eddington ratio when using the FWHM-based scaling relationships [116, 126], possibly due to a dependence of the broad-line region structure on these quantities. Systematic effects on broad-line region geometry, which can effect the inferred velocity dispersion, are a particular concern for C IV, which is thought to arise in an accretion disk wind [127]. Along these lines, unaccounting for radiation pressure on broad-line clouds may also bias the virial masses, especially among those AGN radiating near the Eddington limit [128]; however, its importance is still debated [129–131]. In addition, the reliability of line width measurements can rapidly deteriorate for low data [132]. And, finally, the BLAGN virial mass estimates are calibrated to the reverberation mapping derived masses, which themselves are calibrated to lie on the local - relationship. Most of the AGN that are used to calibrate the reverberation mapping masses to the - relationship have lower masses and are hosted by late-type galaxies, for which there is evidence that the - relationship begins to break down [81]. Greene et al. [81] argue that the normalization of the scaling relationships inferred when limiting the calibration to low-mass SMBHs hosted in late-type galaxies may be about a factor of ~1.5 lower than that used for the current broad-line mass estimates [81]. However, dynamical mass estimates exist for two reverberation mapped AGN: NGC 3227 [133, 134] and NGC 4151 [134, 135]. In both cases, the masses derived from dynamical modeling and reverberation mapping agree, so it is unclear if a smaller scaling factor is needed for late-type galaxies. These issues show that there are still many remaining questions regarding virial masses, highlighting the need for further study using high-quality reverberation mapping data.
5.2. Other Methods for Estimating the Black Hole Mass Function of AGN
Before broad-line mass estimates, there were two earlier attempts at estimating the BHMF for AGN, which we briefly mention here. Siemiginowska and Elvis [136] and Hatziminaoglou et al. [137] used a model for the AGN lightcurve arising due to thermal-viscous accretion disk instabilities [138] to calculate the expected distribution of luminosity at a given black hole mass. Based on this calculated distribution, they used the quasar luminosity function to constrain the quasar black hole mass function. Siemiginowska and Elvis [136] found evidence for SMBH downsizing in AGN, consistent with later work.
Franceschini et al. [139] found a tight correlation between and the total radio power observed in a sample of local galaxies. They then used their empirical relationship to estimate the local BHMF derived from the local radio luminosity function of galaxies. While many of the objects in their sample are not considered AGN in the traditional sense, Fraceschini et al. [139] argue that this correlation is a signature of an advection-dominated accretion flow, thought to dominate at low accretion rates relative to Eddington. Therefore, while these SMBHs may not be “active” in the quasar sense, the determination of their mass function relies on radio emission from the SMBH accretion flow, so this method may still be considered a method for estimating the BHMF for active SMBHs. Franceschini et al. [139] compared their BHMF to models of AGN activity and found that it was inconsistent with AGN activity being continuous and long lived, but consistent with AGN activity being transient and possibly recurrent.
6. Theoretical Models for Black Hole Mass Functions across Cosmic Time
There have been numerous theoretical models for the formation and growth of supermassive black holes, and coevolution with their host galaxies. Understanding this formation, growth, and coevolution is one of the current most important outstanding issues in extragalactic astrophysics. Because the black hole mass function provides a census of the SMBH population and its evolution, it is one of the most fundamental observational quantities available for constraining models of SMBH formation and growth. As such, many theoretical investigations have predicted a BMHF for comparison with the empirical BHMF. In this section, we review some of the models for SMBH formation and growth. There have been numerous theoretical models for SMBH growth and formation, and it is beyond the scope of this primarily empirically-focused review to review all of them; instead, we focus on those theoretical models that predict a BHMF.
6.1. Modeling the Coevolution of SMBHs and Galaxies: Predicted BHMFs
Early models for the coevolution of SMHBs and galaxies linked the growth of black holes to the properties of host dark matter halos, with periods of SMBH growth occurring in quasar phases initiated by mergers. In general, early studies that predicted a BHMF used various prescriptions to relate to the mass of the host halo [26, 140–143]. More recent models for the coevolution of SMBHs and galaxies has incorporated AGN feedback from the SMBH. In addition, the availability of empirical BHMFs have enabled modelers to compare their more recent models with observational data. In general, the models are qualitatively in agreement with the empirical results, in that they are able to match the local BHMF fairly well and predict downsizing of SMBHs. However, considering the current systematic and statistical uncertainties in the empirical results, it is difficult to place rigorous empirical constraints on the models such that certain models may be ruled out. Because of this, we simply summarize some of the different recent models that have been developed which predict the BHMF.
Granato et al. [144] developed a model incorporating feedback from AGN and supernovae, where the feeding of the SMBH is driven by stellar radiation drag on gas. Their predicted local BHMF agrees with that estimated by Shankar et al. [57]. Cattaneo et al. [30] used halo merger trees constructed from -body simulations to track the growth of SMBHs. In their model, the black hole fueling rate was proportional to the star formation rate of the host galaxy burst component and the density of the cold gas in the starburst component. Their model predicted SMBH downsizing, with the most massive part of the BHMF being built up first, in agreement with the subsequent empirical studies.
Hopkins et al. [145] describe a model for the coevolution of SMBHs and galaxies whereby all major mergers of gas-rich galaxies trigger a quasar. In this model, the final black hole mass is assumed to be on average proportional to the host spheroidal mass, in agreement with the local scaling relationships between SMBHs and their host galaxies. Hopkins et al. [145] estimated the merger rate of gas-rich galaxies by combining theoretical constraints of the halo and subhalo mass functions with empirical constraints on halo occupation models. Their model also predicts SMBH downsizing, and their predicted BHMF matches the local BHMF derived by Marconi et al. [53]. Similarly, Shen [146] also assumed that quasars are triggered by major mergers of gas-rich galaxies, with the SMBHs growing via accretion in these quasar phases. Shen [146] used a halo merger rate based on theoretical expectations from -body simulations and assumed a universal quasar lightcurve shape having an exponential increase followed by a power law decay (see also [62]). The BHMF predicted by Shen [146] broadly agrees with the local one estimated by Shankar et al. [64] and predicts that most SMBHs with were in place by but only 50% of them were assembled by .
Most recently, Fanidakis et al. [147, 148] extended the model of [23], which includes AGN feedback, to also follow the spin distribution of SMBHs. In their model, SMBHs are fueled through accretion of cold gas from mergers, disk instabilities, and cooling flows from hot halos. However, the inclusion of SMBH spin enabled them to include different radiative efficiencies, which dictates how much accreted material actually grows the black hole, and to provide an improved model for the amount of mechanical feedback imparted through an AGN jet, both of which depend on the spin of the black hole. Their model predicts that the present-day universe is dominated by SMBHs with –, and that the BHMF at was largely built up at due to an increase in both lower-accretion rate “radio-mode” growth and mergers of SMBHs.
Almost all models for the cosmological coevolution of SMBHs and galaxies that predict a BHMF have been of an analytical or semi-analytical nature. An exception is the study done by Di Matteo et al. [31], who present the results from cosmological hydrodynamic simulations of the ΛCDM model that follow the growth of galaxies and SMBHs, including their feedback processes, at . Direct cosmological simulations such as these should, in principle, provide the most accurate results as to the predicted BHMF, and for identifying the relevant physical processes that are important in shaping the BHMF. However, current cosmological hydrodynamic simulations suffer from the fact that they cannot resolve processes on physical scales corresponding to the SMBH accretion flow. In fact, Di Matteo et al. [31] use a gravitational softening length of h−1 kpc. Instead, Di Matteo et al. [31] employ a subresolution model where the accretion onto the SMBH is estimated using a Bondi-Hoyle-Lyttleton parameterization [149–151] with a correction factor to account for the fact that the Bondi radius is not resolved. They assume a radiative feedback energy efficiency of [20], which is the only free parameter in their model and required in order to match the normalization of the observed local - relationship. Their calculated BHMF at matches the local BHMF for . In addition, Di Matteo et al. [31] also find downsizing in their model, in agreement with observations, with the high-mass end of the BHMF being largely in place by .
In Figure 6, we compile predicted BHMFs from several recent models for SMBH formation and growth [31, 145–147, 152]. In general, the models tend to agree to within a factor of a few with regards to the BHMF. However, they diverge at , where their predicted SMBH number densities can differ by over an order of magnitude.
(a)
(b)
(c)
(d)
6.2. Modeling the BHMF of SMBH Seeds
Recent work has made improvement to models for the BHMF by focusing on theoretical modeling of the distribution of seed SMBHs. The discovery of quasars at –7 with [153, 154] places strong constraints on the formation of SMBH seeds due to the very limited amount of time available at that redshift to grow SMBHs (e.g., [155]). Lodato and Natarajan [156] derive the BHMF of SMBH seeds at that are the result of the collapse of pregalactic disks which have not yet been enriched by metals [157]. Black holes formed through such a mechanism have masses , while black holes which are the remnants of Pop III stars have . A similar seed black hole formation mechanism is through “quasistars” [158–160], which are also able to produce seed black holes with .
Volonteri et al. [161] describe a model for the growth of SMBHs seeded according to the direct collapse model of Lodato and Natarajan [157] with varying formation efficiencies. In addition, they also compared the results from this model using SMBHs seeded from Pop III remnants. Volonteri et al. [161] grow SMBHs through major mergers, and force the black hole mass after the galaxy merger to scale with the circular velocity of the host halo; additional growth is also provided through black hole mergers. Their merger trees are based on a Monte Carlo algorithm based on the extended Press-Schechter formalism. They find that most significant differences in the local BHMF with respect to black hole formation efficiency occur at , with the number density of SMBHs with increasing with increasing formation efficiency. Volonteri and Begelman [152] performed a similar analysis as that of Volonteri et al. [161] but instead used SMBH seeds formed via quasistars. The BHMFs calculated by Volonteri and Begelman [152] match those of Merloni and Heinz [63] at the high-mass end, at least at .
Natarajan and Volonteri [162] used a growth and seeding model which is very similar to that employed by Volonteri et al. [161]. However, they also predict the BHMF for broad line quasars, assuming that 20% of quasars are unobscured. They compare the BHMF derived from their model at to the BHMF for broad line quasars reported by Kelly et al. [115], and to the BHMF for all SMBHs reported by Merloni and Heinz [63]. Natarajan and Volonteri [162] concluded that seeds from Pop III stars have difficulty reproducing the BLAGN BHMF, especially at high redshift, while seeds resulting from the direct collapse of pregalactic disks do better at fitting the high mass end of the BLAGN BHMF at .
Lippai et al. [163] studied the impact that the fraction of halos that form a SMBH seed at has on the local BHMF. They are able to reproduce the observed quasar luminosity function and local BHMF for a suitable range of radiative efficiency and quasar lifetime, so long as at least ≈10% of high- halos contained SMBH seeds. Tanaka and Haiman [164] present BHMFs at based on a comprehensive exploration of the parameter space governing the buildup of the BHMF, including models for SMBH seeds that form from Pop III remnants and through direct collapse, variations in the occupation fraction, and a detailed treatment of gravitational recoil. Comparing their predicted BHMFs to observations of quasars from the SDSS, they concluded that ~100 seeds can grow into the SMBHs observed at , so long as they are nearly continuously embedded in dense gas, form at , have low occupation fractions, and stop forming by .
7. Directions for Future Work and Improvement
Before concluding this paper, we present a discussion of possible future empirical and theoretical work relevant to BHMF studies. These include the following.
7.1. Better Characterization of the SMBH-Host Galaxy Scaling Relationships
Currently, the local BHMF is estimated from the distribution of host galaxy properties assuming that has a constant log-normal scatter about a single-power law scaling relationship. As discussed in Section 3, recent observations have provided reason to doubt this assumption, suggesting that the correlations break down at the highest and lowest masses. This will create biases in the BHMF determined from the scaling relationships, which in turn will also affect the BHMF estimated from the continuity equation. Further direct estimates from dynamical and kinematic modeling should be obtained for a variety of galaxy types, especially at the high- and low- end. The next class of 25+ m telescopes should provide a significantly improved picture of the scaling relationships, thus providing us with more accurate estimates of the local BHMF.
7.2. Improvements to Techniques Based on the Continuity Equation
Most studies that have invoked the continuity equation to link the local BHMF to the AGN luminosity function have assumed a single radiative efficiency, which is equivalent to assuming a single black hole spin, and a universal AGN lightcurve. Neither of these assumptions are likely to be true, and improvements to this type of modeling should include a distribution of SMBH spin and AGN lightcurves. In addition, we need to better characterize the bolometeric corrections, which remain a significant source of systematic uncertainty. The continuity equation techniques should also be extended to map the evolution of the full joint 3-dimensional distribution of black hole mass, accretion rate, and spin. While this will not necessarily have a direct effect on estimating the BHMF, it will provide insight into the dominant accretion modes experienced by active SMBH and into the dominant fueling mechanism for AGN activity, as the spin distribution traces the SMBH fueling history [165].
7.3. Better Characterization of Scaling Relationships Involving for AGN
The dominant scaling relationship for estimating in AGN involves the broad emission lines. However, as discussed in Section 5.1, a number of systematic uncertainties remain. In order to reduce these systematics, we need to better understand the broad-line region geometry for the different emission lines and how it scales with luminosity, which will provide us with a more accurate conversion from line width to velocity dispersion. Moreover, accurate characterization of the broad-line region geometry should remove the need to calibrate the scaling relationships to the local - relationship, which has its own set of systematics. Improvements to reverberation mapping campaigns and modeling [166–170], as well as increasing the number of AGN monitored for reverberation mapping, will be needed in order to really understand the systematics involved in the broad line mass estimates.
There is also the need to better characterize the black hole fundamental plane. Because the BHFP describes how the emission mechanisms responsible for the radio and X-ray flux scale with , the BHFP coefficients depend on these emission mechanisms. However, these emission mechanisms depend on the geometry of the accretion state and the existence of a jet, which in turn depends on the accretion rate [171], and therefore the BHFP coefficients will be different for different classes of AGN [172–174]. In particular, the BHFP is currently poorly constrained for “soft-state” galactic black holes and radio-quiet AGN. Therefore, to reduce the systematics involved with the BHFP, it will be necessary to characterize the scaling relationships and their scatter for radio-quiet objects. A correlation between the radio and X-ray luminosity has been observed for radio-quiet objects [175], implying that a BHFP should also exist for these objects. In order to better characterize the BHFP for radio-quiet objects it will be necessary to obtain radio detections for a well-defined sample of radio-quiet AGN with reliable estimates and X-ray detections.
Finally, there has recently been the discovery of scaling relationships involving and the optical [176, 177] and X-ray [42, 45, 49] variability properties of AGN. Mass estimators based on these scaling relationships have not been rigorously developed yet, nor have they seen widespread use. However, the existence of these scaling relationships implies that the variability properties may offer another avenue for estimating and BHMFs, which may become increasingly valuable in the era of current and future large time-domain surveys, such as Pan-Starrs and LSST.
7.4. Understanding the Redshift Evolution of Scaling Relations
From the theoretical point of view, it is clear that high-redshift scaling relations (or the lack thereof) between SMBH and their hosts provide unique and powerful constraints to models for AGN feeding and feedback, which cannot be otherwise distinguished (see, e.g., Merloni et al. [74] and references therein).
In practical terms, a better understanding of the evolution of scaling relations may also be very advantageous for BHMF studies. As we discussed above, current technique for BH mass estimation at involve unobscured, broad emission line QSOs. One can argue that, as long as we are restricted to just this class of QSOs, we will have to make critical assumptions about the properties of a significant part of the population to draw conclusions about the full BHMF (if we wanted, e.g., to compare with “continuity-equation-based” methods). On the other hand, large multiwavelength surveys do and will provide a wealth of information on the host galaxies of obscured AGN at high redshift, that represent the numerically dominant part of the growing black holes population [178–181]. Therefore, if we had an independent way to put constraints on the nature of the BH-host relation for these objects, we could explore the uncharted territory of BHMF for obscured AGN (and for the entire population). Such an independent information could come, for example, from IR studies of broad emission lines which could act as probes of the BH potential less affected by obscuration. The first exploratory works pursuing this line of research have recently been published, for example, [182, 183].
From the technical point of view, a lot of work of course is needed to better understand how reliable these estimators are. Another big “technical" challenge of all studies of the evolution of scaling relations is the fact that they require a thorough assessment of the many observational biases one encounters in studying high redshift AGN and their hosts [184].
7.5. Accounting for Black Hole Kicks in Theoretical Models
Many theoretical models for the BHMF do not include recoiling effects caused by the merger of two black holes. However, recent theoretical work on black hole recoils suggests that black holes can spend a significantly large enough amount of time offset from the central region of the host galaxy to alter their growth, thereby increasing the scatter about the scaling relationships and decreasing the final black hole mass [185–188]. On the other hand, Volonteri et al. [189] and Volonteri et al. [190] find that ejected SMBHs are rare at , especially for massive SMBHs, suggesting that accounting for ejected black holes will not make a significant difference in the BHMF. Moreover, Tanaka and Haiman [164] find that recoiling black holes only modestly effect the BHMF. Further improvement in our understanding of the effects and frequency of black hole recoil will ensure the accurate implementation of black hole recoil into models for SMBH growth.
7.6. Improvements to Our Understanding of AGN Feedback
Most current theoretical models for SMBH growth involve AGN feedback and assume a single efficiency for coupling feedback energy to the gas; this feedback efficiency is usually treated as a free parameter. An improved physical understanding of AGN feedback will improve theoretical models for the BHMF, as the feedback efficiency affects the dynamics of the SMBH’s fuel supply and therefore the amount that the SMBH accretes as a function of redshift. Recent high-resolution hydrodynamic simulations in one dimension [191–193] and two dimensions [194, 195] have concluded that AGN feedback efficiency increases with the Eddington ratio and that the values are below the value of ~5% assumed in many current theoretical models for SMBH growth. Further improvements to simulations developed for studying AGN feedback will lead to a better physical understanding of AGN feedback, which will improve theoretical models for SMBH growth and the BHMF.
On the observational side, future X-ray observations should provide considerable improvement in our understanding of AGN feedback. X-ray spectra are needed in order to determine the total column density of the gas, and thus its kinetic energy flux, which can be compared to the energy output of the SMBH. Current X-ray observations from Chandra have found evidence that AGN feedback exists in the local universe [196]. However, X-ray calorimeters on future X-ray satellites will be needed for further improvement as they provide the high throughput and spectral resolution needed to measure column densities and velocities of ionized gas, and consequently the kinetic energy flux, in a large sample of AGN across a broad redshift range.
7.7. Improvements in Resolution and Subresolution Modeling for Direct Hydrodynamic Simulations
Full hydrodynamic cosmological simulations offer the most promising avenue for providing a physically motivated BHMF without free parameters and for unambiguously identifying the relevant physical processes in building up the BHMF. However, they currently cannot resolve scales relevant to the accretion flow onto the SMBH. Numerical codes based on adaptive mesh refinement techniques will provide improvement in resolution, but it will likely be a while before hydrodynamic cosmological simulations are able to follow SMBH growth in large cosmological volumes while simultaneously resolving the scales relevant for individual black holes. In the meantime, further improvement can be made to the subresolution modeling employed by current hydrodynamic simulations.
One way of improving current sub-resolution models may be to implement the results on AGN feedback based on the type of work described in the previous bullet point. Another improvement is in modifying the sub-resolution model for the SMBH accretion flow. Current methods assume the Bondi rate combined with a correction factor to account for the fact that the temperature and density of the gas are not resolved at the Bondi radius. Not surprisingly, the growth of the SMBH is sensitive to how this correction factor is modeled [33]. Moreover, sub-resolution models based on the Bondi rate neglect the angular momentum of the gas, and thus the Bondi rate may not be representative of the actual accretion rate onto the SMBH. Hopkins and Quataert [197] used high-resolution hydrodynamic simulations to conclude that the Bondi rate was a poor estimate of the actual accretion rate onto the SMBH, and describe a sub-resolution model which accurately estimated the actual accretion rate in their simulations. In addition, Power et al. [198] suggest an alternative sub-resolution model based on an “accretion particle” to provide a more accurate estimate of the black hole accretion rate. The implementation of improved sub-resolution models for accretion rate and feedback into cosmological hydrodynamic simulations, as well as further improvements to the sub-resolution models, will result in more accurate predicted BHMFs, allowing a more insightful comparison with empirical BHMFs.
In this paper we have reviewed current estimates of the SMBH mass function, as well as theoretical models for the BHMF. As discussed above, each of the methods for estimating the BHMF has their own set of systematics. In Figure 7, we compare the empirical estimates of the local BHMF (defined by the shaded region in Figure 2) with the BHMFs predicted by the theoretical models compiled in Figure 6. In addition, in Figure 7, we also compare the empirical BHMFs at , as estimated using the lightcurve method (shown in Figure 3) and the black hole fundamental plane (shown in Figure 4), with BHMFs predicted by the theoretical models. In both Figures we also include the BHMFs for broad-line AGN as estimated by Schulze and Wisotzki [121] and Kelly et al. [115]. In general, the theoretical models are consistent to within a factor of a few with the empirical estimates of the BHMF, although there is a large spread in the models and empirical estimates at . Moreover, the estimated number densities of broad-line AGN are significantly lower than those of all SMBHs, suggesting that only a small fraction of SMBHs are active across a broad range in , except for possibly SMBHs at with .
(a)
(b)
Despite the differences in the methods for estimating the BHMF and the theoretical models, they have lead to a number of common conclusions. In particular, the empirical results have presented a picture whereby SMBHs grow primarily via accretion in active phases (Eddington ratios ), that quasar activity is a relatively short-lived phenomenon relative to the lifetime of the SMBH and host galaxy (i.e., small “duty cycles” for AGN activity) and that SMBH growth is anti-hierarchical with the most massive end of the BHMF being built up first. These empirical results are qualitatively in agreement with the steadily improving theoretical models.
Acknowledgments
The authors would like to thank Tommaso Treu, Priya Natarajan, Marta Volonteri, Aneta Siemiginowska, Xiaohui Fan, Marianne Vestergaard, Alister Graham, and Zoltán Haiman for helpful comments. In addition, they would like to thank two anonymous referees whose comments improved our paper. B. C. Kelly acknowledges support by NASA through Hubble Fellowship Grant no. HF-51243.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under Contract NAS 5-26555.