Mass Functions of Supermassive Black Holes Across Cosmic Time

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 review 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.


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 [93,108,109], stellar velocity dispersion (the M BH -σ * relationship, e.g., [47,116,162]), concentration or Sersic index [50,49], bulge mass [101,104,65], and binding energy [4,72]. 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 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 [145,40,6,118]. 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. [128,78].
Feedback-driven 'self-regulated' growth of black holes has been able to reproduce the local M BH -σ * relationship in smoothed particle hydrodynamics simulations [38,159,80]. Moreover, AGN feedback has also been invoked as a means of quenching the growth of the most massive galaxies [20,33]. There have been numerous models linking SMBH growth, the quasar phase, and galaxy evolution [63,82,62,182,175,26,37,158,19]. 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 [106], and accretion of ambient gas via internal galactic processes [29,68], may fuel these fainter, lower M BH AGN at lower z. This is supported by the fact that many AGN are observed to live in late-type galaxies out to z ≈ 1 [61,46], and the X-ray luminosity function of AGN hosted by late-type galaxies suggests that fueling by minor interactions or internal instabilities represents a non-negligible contribution to the accretion history of the Universe [48].
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 [113,42,107,35,156,86]. 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 review we discuss the current status of BHMF estimation and theoretical modeling. In § 2 we discuss the non-trivial task of estimating the BHMF. In § 3 we discuss current estimates of the local SMBH BHMF. In § 4 we discuss BHMF estimates derived by combining the local BHMF with the AGN luminosity function via a continuity equation. In § 5 we discuss BHMFs estimated for AGN only. In § 6 we review theoretical models for SMBH growth that predict the SMBH BHMF. Finally, in § 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.

Estimating the Black Hole Mass Function
The black hole mass function, denoted as φ(M BH , z)dM BH , is the number of sources per comoving volume V (z) with black hole masses in the range M BH , M BH + dM BH . The black hole mass function is related to the joint probability distribution of M BH and z, p(M BH , z), as The normalization of the BHMF is N , the total number of SMBHs in the observable universe, and is given by the integral of φ over M BH and V (z).

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 [87]. 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, e.g., 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 be 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 M BH , 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 M BH and the properties of the host galaxy bulge or the luminosity and the width of the broad emission lines for AGN [178,170]. It has also recently been suggested that the X-ray variability properties of AGN may also provide another scaling relationship for estimating M BH [123,188,86], but further work is needed for developing this. While these scaling relationships enable one to estimate M BH for large numbers of SMBHs, they also contain a significant intrinsic statistical scatter. Gültekin, K., et al. [60] find that for early type galaxies there is an intrinsic scatter in M BH of 0.31 ± 0.06 dex and 0.38 ± 0.09 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 & Peterson [170] estimate the scatter in M BH 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 M BH . This is not the case for M BH estimated through dynamical modeling, however independence between M BH and its error is likely to be a good approximation for M BH estimated using scaling relationships. Because scaling relationships are the only feasible manner in which to estimate M BH for a large sample of SMBHs, which is necessary for any estimate of the BHMF, we will assume that M BH and its error are independent. Under the assumption of independence between the estimated M BH 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 than the intrinsic mass function. In particular, the distribution of the mass estimates falls off more slowly with increasing M BH , and overpredicts the number of SMBHs at the high M BH end of the mass function. The bias is worse when the dispersion in the scatter in the mass estimates becomes larger.

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 1 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 M BH and the observable quantities. We can make the above procedure more quantitative by deriving the likelihood function for the SMBH mass function. Kelly, Vestergaard, & Fan [87] 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 asM BH . In addition, denote as X 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 n SMBHs is where the BHMF is related to N and the probability distribution of M BH and z via Equation (1). Here, C N n is the binomial coefficient, θ denotes the parameters for the BHMF, ψ denotes the parameters for the distribution in X at fixed M BH and z, and s(θ, ψ) is the probability of including a SMBH in one's sample as a function θ and ψ. Here, we have assumed that the dis-tribution in the mass estimates at fixed M BH , z, and X, p(M BH |M BH , z, X), 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 2 with N 'trials' and probability of success s(θ, ψ). The probability of including a SMBH in one's survey as a function of the BHMF, s(θ, ψ), is calculated from the survey selection function s(X, z) as (3) It is up to the researcher to choose the particular parameteric form for the SMBH mass function, the distribution in the mass estimates at fixed M BH , z, and X, and the distribution of the observable that the sample is selected on (e.g., flux) at fixed M BH and z. Typical choices are log-normal distributions, Schechter functions, and mixtures of log-normal distributions. Once one has done this, one can use Equation (2) to compute a maximumlikelihood 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, Y , and the intrinsic distribution of Y is known. This is commonly used to estimate the local mass function using host-galaxy scaling relationships [105]. In this case, the mass function is where φ(Y ) is the comoving number density of SMBHs as a function of the quantity Y . When both p(M BH |Y ) and φ(Y ) are known then the BHMF follows directly from Equation (4). As an example, if the mass function is derived from the scaling between M BH and host galaxy spheroidal luminosity, L sph , then Y = L sph , p(M BH |Y ) is the M BH -L sph relationship, and φ(Y ) is the luminosity function of stellar bulges hosting SMBHs. As with BHMFs determined from a mass estimator, improper treatment of the intrinsic scatter in M BH at fixed Y will lead to a biased estimate of the BHMF. However, when calculating the BHMF from Equation (4), ignoring 2 Often in the luminosity function literature the likelihood is assumed to be a Poisson distribuiton. A poisson distribution is an approximation to the binomial distribution when N → ∞ and s(θ, ψ) → 0, so Equation (2) converges to the Poisson distribution when n ≪ N . See [85] for further details.
the intrinsic scatter results in an estimated BHMF that is too narrow, underpredicting the number of SMBHs at the high mass end of φ(M BH ). This is opposite to the case when one estimates the BHMF directly from mass estimates.

Black Hole Mass Functions Derived from Host Galaxy Scaling Relationships
The observed scaling between M BH and the properties of the SMBH host galaxy bulge have motivated several groups to estimate the local BHMF [136,187,3,105,146,96,49,165,185,112,147,171], 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 § 4. Typically, the local BHMF is estimated using the local M BH -σ * relationship or the local M BH -L sph relationship, combined with the local number density of galaxies as a function of stellar velocity dispersion or bulge luminosity. The scaling relationships between M BH 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. An exception is Tamura, Ohta, & Ueda [161], who estimated the BHMF out to z ≈ 1 assuming that evolution in the M BH -L sph relation is driven only by passive evolution in L sph . 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 z [163,129,164,180,115,9], at least for active SMBHs. However, there are still concerns regarding potential biases due to selection effects [97], but see Treu et al. [163] and Bennert et al. [9] 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 z [151]. 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 M BH -σ * 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. [153], with an additional component representing the brightest cluster galaxies [96]. Sheth et al. [153] estimate the velocity dis-  [147]. Based on this estimate the local universe is dominated by SMBHs with M BH < 10 7 M ⊙ . persion 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 σ = v c / √ 2. When the M BH -L sph relation is used it is typical to estimate the distribution of L sph seperately 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 M BH 4 × 10 7 M ⊙ [185]. Shankar, Weinberg, & Miralda-Escudé [147] 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 M BH -σ * , M BH -L sph , and M BH -M star relationships, as presented in Shankar et al. [147]. In general, estimates of the local mass density of SMBHs span the range ρ BH = (3.2-5.4) × 10 5 M ⊙ Mpc −3 for h = 0.7 [147].
While the procedure for estimating the local BHMF is, in theory, straight-forward, a number of significant systematics remain. First, there is the observational difficulty that most BHMFs derived from the M BH -σ * 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 M BH < 10 7 M ⊙ . Another concern is that the local BHMF is derived by assuming that the M BH -σ * or M BH -L sph relations are single power-laws with a constant scatter in M BH at fixed σ * or L sph . However, recent work has shown these assumptions to be incorrect. For one, the M BH -σ * and M BH -L sph relations diverge at the high M BH end, which Lauer et al. [96] suggest implies that the M BH -σ * relation is not a single power-law. This divergence creates an inconsistency in the BHMFs derived from these two scaling relationships [96,165]. Similarly, the σ-L relationships for the SDSS and dynamical M BH SMBH samples are inconsistent, suggesting a possible selection bias in the estimated BHMFs [187,13]. The scatter in the M BH -σ * relation is larger for spirals [60], and appears to increase at low M BH such that most SMBHs lie below the M BH -σ * relation, [57, e.g.]. Several authors have found differences in the slope and scatter of the scaling relations for pseudobulges [77,56,51]; 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. Kormendy, Bender, & Cornell [92] argue that M BH does not correlate with galaxy disks, and only correlates weakly, if at all, with pseudobulges. Although there is still much that we do not understand about the M BH and host galaxy scaling relationship, these recent results suggest that the scaling relationships are not a single powerlaw with constant intrinsic dispersion in M BH , representing a significant source of systematic uncertainty in the estimated local BHMF, especially at the low-mass end.

Black Hole Mass Functions Derived from the Local Mass Function and the AGN Luminosity Function
By employing the argument of Soltan [157], 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 [187,184,105,111,73,112,24,23]. These methods em-ploy a continuity equation describing the evolution of the number density of SMBHs [27,155]: Here, Ṁ (M BH , t) is the average growth rate of SMBHs as a function of M BH and cosmic age, t(z), and n merge (M BH , t) 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, n merge (M BH , t) can also include a contribution from SMBHs which are created, but this has not been thought to occur over the redshift range in which Equation (5) is typically applied, i.e., z 5. Because the merger rate of black holes is currently unknown, many studies that have employed Equation (5) Under the assumption that SMBHs grow during phases of AGN activity, AGN demographics in combination with the local BHMF may be used to compute φ M (M BH , t). This is because the AGN luminosity function maps the accretion history onto SMBHs, and the local BHMF acts as a boundary condition on Equation (5); it is also possible in principle to include the BHMF for AGN, which provides more information. Studies that have used Equation (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 seperately.

Methods that Assume an AGN Lightcurve
Most authors employing Equation (5) have assumed a parameteric form for Ṁ (M BH , t) . The accretion rate is related to the bolometric luminosity output of the accretion flow onto the SMBH as L = ǫ rṀacc c 2 , where ǫ r is the radiative efficiency of the accretion flow,Ṁ acc is the accretion rate of matter onto the SMBH, and c is the speed of light. The growth rate of the SMBH isṀ = (1 − ǫ r )Ṁ acc , due to the fact that a fraction ǫ r 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 an time given the local BHMF, an assumed average accretion flow lightcurve as a function of M BH , L(M BH , t) , and an assumed radiative efficiency. Because φ M (M BH , z) and L(M BH , t(z)) imply a luminosity function, the local BHMF and AGN luminosity function can be used to place constraints on ǫ r and L(M BH , t) . This means that, in practice, one also has to assume a bolometeric correction, which itself likely depends on both black hole mass [84] and L/L Edd [166,167]. In addition, an estimate of L(M BH , t) 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 M BH and z.
A variety of lightcurve models have been used when employing Equation (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 [136,187,105,146,147] 3 to study the build-up of the local black hole mass function, although [147] also considered models where the average accretion rate relative to Eddington falls off toward lower z and higher M BH . Raimundo & Fabian [133] employed a variation on the constant L/L Edd models, assuming three different populations of AGN with their own Eddington ratio: a population of obscured low L/L Edd AGN, a population of obscured AGN with higher L/L Edd , and a population of unobscured AGN. Yu & Lu [185] modeled the quasar lightcurve as radiating at the Eddington limit for a period of time, and then transitioning into a power-law decay. Cao [23] 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 [68,69] or from viscous evolution of the accretion disk [186,89].
In Figure 3 we compare the BHMF calculated by Shankar et al. [147] with that calculated by Cao [23]. For Shankar et al. [147] we show their reference model, which assumes a radiative efficiency of ǫ r = 0.065, an accretion rate relative to Eddington ofṀ /Ṁ Edd = 0.6, and that half of all SMBHs are active at z = 6. We show the model from Cao [23] which assumes a radiative efficiency of ǫ r = 0.11 and a quasar lifetime of 7.5 × 10 8 yr, as it better matches the Shankar et al. [147] estimates. The two estimates of the BHMF agree fairly well, despite the different quasar lightcurve models.
In general, most of the studies that have used Equation (6) in combination with an assumed quasar lightcurve have concluded the following: • Most SMBH growth occurs in periods when the quasar is radiating near the Eddington limit.
• 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.
• SMBH growth is anti-hierarchical, with the most massive black holes growing first. This has also been termed 'down-sizing' of active SMBHs.
• The lifetime of AGN activity is ∼ a few × 10 8 yr.
• Most SMBHs have non-zero spin, as implied by inferred radiative efficiencies of ǫ r 0.06.
However, while Equation (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 use of Equation (6) often entails some strong assumptions. These methods rely on the assumed form of the quasar lightcurve, distribution of radiative efficiencies, and bolometeric 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 § 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 ǫ r 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 Equation (6), and further work is needed in reducing these systematics.

Methods that Include the Distribution of Active Supermassive Black Holes
An alternative to the methods described in § 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 L(M BH , t) is derived directly from the estimated distribution of L/L Edd . 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 [111] and Merloni & Heinz [112], who employed the black hole 'fundamental plane' (BHFP) [113,42]. The BHFP is a scaling relationship between M BH , radio luminosity, and X-ray luminosity, that exists for low-accretion rate black holes (i.e., M /Ṁ Edd 0.01), extending from galactic black holes to supermassive ones. The BHFP likely reflects the connection between M BH 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 M BH and Xray luminosity at a given redshift for the active SMBH population, Merloni [111] and Merloni & Heinz [112] then convert this to a joint distribution  [112], showing the redshift evolution of their estimated BHMF. The dashed line is the local BHMF and the shaded regions reflect the uncertainty in the BHMF that is due to uncertainties in the AGN luminosity function. The high mass end of their estimated BHMF is built up faster than the low mass end, a phenomenon that has been called 'downsizing'.
of M BH andṀ acc assuming a conversion from X-ray luminosity toṀ acc which depends on the Eddington ratio. The joint distribution of M BH anḋ M acc at a given redshift for active SMBHs therefore enables calculation of the average growth rate Ṁ (M BH , t(z)) , 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 & Heinz [112] concluded that SMBHs grow antihierarchically; however, in contrast to the lightcurve methods, Merloni & Heinz [112] concluded that most SMBHs have low spin as inferred from their derived radiative efficiency. In addition, Merloni & Heinz [112] 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 [111] and Merloni & Heinz [112] 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 & Heinz [112], the BHFP is only defined for low-accretion rate objects, i.e., objects with L/L Edd 10 −2 . Merloni & Heinz [112] 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 L/L Edd 10 −2 and lacking a jet). Unfortunately, the AGN in the radiatively efficient mode make a significant contribution to the X-ray luminosity function, from which Ṁ (M BH , t) is derived. Moreover, most studies, including those based on the BHFP have concluded that most SMBH growth occurs at L/L Edd 10 −2 , 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 M BH and the AGN SED.

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.

Methods Based on Scaling Relationships Involving the Broad Emission Lines
The steady improvement in reverberation mapping of AGN [130,10] has revealed a correlation between the luminosity of AGN and the broad line region radius [81,11]. It is therefore possible, in principle, to obtain an estimate of M BH for broad line AGN (BLAGN) by combining a luminositybased 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 [178]. These virial mass estimates are then calibrated to the estimates of M BH obtained from reverberation mapping, which themselves are calibrated to be consistent with the local M BH -σ * relationship [125,181]. Currently, calibrations exist for Hα [53], Hβ [170, e.g.], Mg II [110,169,149], and C IV [170]. The statistical scatter in the virial mass estimates is currently estimated to be ∼ 0.4 dex [170], although there are indications that the scatter may be smaller, at least for the most luminous quasars [90,150,160,88,152]. 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 [170]. 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 M BH , 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 1/V max correction [179,54,168,169], a technique borrowed from luminosity function estimation. Greene & Ho [54] estimated the local BHMF for BLAGN from the SDSS DR4, while Vestergaard et al. [168] estimated the BHMF for BLAGN over 0.3 < z < 5 using the uniformly-selected quasar sample from the SDSS DR3 [134]. Vestergaard & Osmer [169] 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 § 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 1/V max correction to correct for incompleteness in M BH . Subsequent attempts have further improved in their methodology, providing more accurate BHMFs.
Shen et al. [150] 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 & Wisotzki [139] employed a maximum-likelihood technique for estimating the local BHMF for BLAGN. Their method corrects for incompleteness in M BH but does not correct the BHMF for the broadening caused by the statistical scatter in the virials mass estimates. Kelly et al. [87] 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 [138]. Kelly et al. [88] used the method of [87] to estimate the BHMF of BLAGN at 1 < z < 4.5 from the mass estimates in the SDSS DR3 quasar sample [168]. 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 z = 4.5. More recently, Shen & Kelly [152] extended the Bayesian method of [87] to include a possible luminosity-dependent bias in virial mass estimates derived from the emission line F W HM , the existence of which was suggested by Shen & Kelly [151]. Shen & Kelly [152] 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, i.e., 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. [88] used the BLAGN BHMF to estimate the lifetime of broad line quasar activity to be t BL ∼ 150 Myr among SMBHs with M BH ∼ 10 9 M ⊙ , which is similar to quasar lifetimes inferred from the continuity equation. Kelly et al. [88] also used their estimated BHMF to estimate the maximum mass of a SMBH to be M BH ∼ 3 × 10 10 M ⊙ , which is in agreement with theoretical expectations [119,143].
Mass functions estimated from scaling relationships for BLAGN have the advantage that they are derived from estimates of M BH 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 & Ho [55]. 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 F W HM -based scaling relationships [32,152], 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 [135]. Along these lines, unaccounted for radiation pressure on broad line clouds may also bias the virial masses, especially among those AGN radiating near the Eddington limit, [102]; however, its importance is still debated [121,103,122]. In addition, the reliability of line width measurements can rapidly deteriorate for low S/N data [36]. And finally, the BLAGN virial mass estimates are calibrated to the reverberation mapping derived masses, which themselves are calibrated to lie on the local M BH -σ * relationship. Most of the AGN that are used to calibrate the reverberation mapping masses to the M BH -σ * relationship have lower masses and are hosted by late type galaxies, for which there is evidence that the M BH -σ * relationship begins to break down [57]. Greene et al. [57] 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 [57]. However, dynamical mass estimates exist for two reverberation mapped AGN: NGC 3227 [34,67] and NGC 4151 [124,67]. 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.  [45] found a tight correlation between M BH 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. [45] 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 SMBH 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. [45] 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.

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.

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 occuring in quasar phases initiated by mergers. In general, early studies that predicted a BHMF used various perscriptions to relate M BH to the mass of the host halo [64,82,75,28,39]. More recent models for the coevolution of SMBHs and galaxies have 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. [52] 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. [146]. Cattaneo et al. [26] used halo merger trees constructed from N -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. [70] 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. [70] 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. [105]. Similarly, Shen [148] also assumed that quasars are triggered by major mergers of gas-rich galaxies, with the SMBHs growing via accretion in these quasar phases. Shen [148] used a halo merger rate based on theoretical expectations from N -body simulations, and assumed a universal quasar lightcurve shape having an exponential increase followed by a power-law decay (see also [185]). The BHMF predicted by Shen [148] broadly agrees with the local one estimated by Shankar et al. [147] and predicts that most SMBHs with M BH > 3 × 10 8 M ⊙ were in place by z = 1, but only 50% of them were assembled by z = 2.
Most recently, Fanidakis et al. [43,44] extended the model of [20], 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 M BH ∼ 10 7 -10 8 M ⊙ , and that the BHMF at M BH > 10 9 M ⊙ was largely built up at z < 2 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. [37], 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 z > 1. 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 [37] use a gravitational softening length of ǫ = 2.73h −1 kpc. Instead, Di Matteo [37] employ a subresolution model where the accretion onto the SMBH is estimated using a Bondi-Hoyle-Lyttleton parameterization [76,18,17] with a correction factor to account for the fact that the Bondi radius is not resolved. They assume a radiative feedback energy efficiency of 5% [38], which is the only free parameter in their model and required in order to match the normalization of the observed local M BH -σ * relationship. Their calculated BHMF at z = 1 matches the local BHMF for M BH > 2 × 10 8 M ⊙ . In addition, Di Matteo [37] also find downsizing in their model, in agreement with observations, with the high mass end of the BHMF being largely in place by z ∼ 2.
In Figure 6 we compile predicted BHMFs from several recent models for SMBH formation and growth [70,148,43,37,173]. In general, the models tend to agree to within a factor of a few with regards to the BHMF. However, they diverge at M BH 10 9 M ⊙ , where their predicted SMBH number densities can differ by over an order of magnitude.

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 z ≈ 6-7 with M BH ∼ 10 9 M ⊙ [79,117] places strong constraints In general the number densities predicted by the models agree to within a factor of a few, although they diverge at M BH 10 9 M ⊙ . Some of these authors did not report BHMFs at each redshift shown, so we only show those available at each redshift. on the formation of SMBH seeds due to the very limited amount of time available at that redshift to grow SMBHs. Lodato & Natarajan [99] derive the BHMF of SMBH seeds at z ∼ 15 that are the result of the collapse of pregalactic disks which have not yet been enriched by metals [98]. Black holes formed through such a mechanism have masses M BH ∼ 10 5 M ⊙ , while black holes which are the remnants of Pop III stars have M BH ∼ 10 3 M ⊙ . A similar seed black hole formation mechanism is through 'quasi-stars' [8,7,5], which are also able to produce seed black holes with M BH ∼ 10 5 M ⊙ .
Volonteri, Lodato, & Natarajan [176] describe a model for the growth of SMBHs seeded according to the direct collapse model of Lodato & Natarajan [98] with varying formation efficiencies. In addition, they also compared the results from this model using SMBHs seeded from Pop III remnants.
Volonteri et al. [176] 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 M BH < 10 7 M ⊙ , with the number density of SMBHs with M BH < 10 7 M ⊙ increasing with increasing formation efficiency. Volonteri & Begelman [173] performed a similar analysis as that of Volonteri et al. [176] but instead used SMBH seeds formed via quasi-stars. The BHMFs calculated by Volonteri & Begelman [173] match those of Merloni & Heinz [112] at the high mass end, at least at z > 2.
Natarajan & Volonteri [120] used a growth and seeding model which is very similar to that employed by Volonteri et al. [176]. 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 1 < z < 4.5 to the BHMF for broad line quasars reported by Kelly et al. [88], and to the BHMF for all SMBHs reported by Merloni & Heinz [112]. Natarajan & Volonteri [120] 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 z > 2.

Directions for Future Work and Improvement
Before concluding this review, we present a discussion of possible future empirical and theoretical work relevant to BHMF studies. These include: • Better characterization of the SMBH-host galaxy scaling relationships. Currently, the local BHMF is estimated from the distribution of host galaxy properties assuming that M BH has a constant log-normal scatter about a single-power law scaling relationship. As discussed in § 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 M BH estimates from dynamical and kinematic modeling should be obtained for a variety of galaxy types, especially at the high and low M BH 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.
• 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 [14].
• Better characterization of scaling relationships involving M BH for AGN. The dominant scaling relationship for estimating M BH in AGN involves the broad emission lines. However, as discussed in § 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 M BH -σ * relationship, which has its own set of systematics. Improvements to reverberation mapping campaigns and modeling [74,12,127,189,21], 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 M BH , 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 [114], and therefore the BHFP coefficients will be different for different classes of AGN [91,59,131]. 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 [95], 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 M BH estimates and X-ray detections.
Finally, there has recently been the discovery of scaling relationships involving M BH and the optical [83,100] and X-ray [123,107,86] 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 M BH and BHMFs, which may become increasingly valuable in the era of current and future large time-domain surveys, such as Pan-Starrs and LSST.
• 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. [115] 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 z¿0 involve un-obscured, 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, for example, to compare with "continuity-equation-based" methods).
On the other hand, large multi-wavelength 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 [22,25,183,1]. 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, e.g., [2,137].
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 [140].
• Accounting for black hole kicks in theoretical models. Most 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 [172,16,144,15]. On the other hand, Volonteri, Gültekin, & Dotti [174] and Volonteri, Natarajan, Gültekin [177] find that ejected SMBHs are rare at z < 5, especially for massive SMBHs, suggesting that accounting for ejected black holes will not make a significant difference in 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.
• 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 usally 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 [30,154,31] and two dimensions [94,126] have concluded that AGN feedback efficiency in-creases 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. Xray spectra are needed in order to determine the toal 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 [41]. 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.
• Improvements in resolution and sub-resolution 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 [19]. Moreover, Figure 7: Comparison of empirical estimates of the BHMF (grey shaded region) with BHMFs predicted by theoretical models for SMBH formation and growth (red shaded region), in both the local universe (left) and at z = 2 (right). Also shown are the estimated BHMF of broad line AGN only (solid green line), from Schulze & Wisotzki [139] (left) and Kelly et al. [88] (right). The empirical estimates of the BHMF are those shown in Figures 3 and 4, while the theoretical estimates are those shown in Figure  6. The shaded regions define the spread in the estimates and models, and may be considered to be a crude estimate of their uncertainty. In general, the theoretical number densities are consistent with the empirical ones to within a factor of a few.
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 & Quataert [71] 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, Nayakshin, & King [132] 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 z = 2, 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 & Wisotzki [139] and Kelly et al. [88]. 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 z = 2. 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 M BH , except for possibly SMBHs at z ∼ 2 with M BH 10 9 M ⊙ .
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 L bol /L Edd > 0.01), 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.
We would like to thank Tommaso Treu, Priya Natarajan, Marta Volonteri, Aneta Siemiginowska, Xiaohui Fan, and Marianne Vestergaard for helpful comments on an earlier draft of this review. In addition, we would like to thank two anonymous referees whose comments improved our review. BK acknowledges support by NASA through Hubble Fellowship grant #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.