format: text

text:
# Observational constraints on growth of massive black holes

## 1 Introduction

Most nearby galaxies contain massive dark objects at their centers (e.g. Kormendy & Richstone 1995, Kormendy & Gebhardt 2001), which are presumably black holes (BHs). The existence of these objects was predicted by arguments based on quasi-stellar object (QSO) energetics and demography (e.g. Sołtan 1982, Rees 1984). Studies of central BHs in nearby galaxies have also revealed a tight correlation between BH mass and galactic velocity dispersion (Ferrarese & Merritt, 2000; Gebhardt et al., 2000), and a less tight correlation between BH mass and the luminosity (or mass) of the hot stellar component of the host galaxy (e.g. Kormendy & Gebhardt 2001 and references therein; by “hot”component we mean either an elliptical galaxy or the bulge of a spiral or S0 galaxy). These correlations strongly suggest a close link between the formation and evolution of galaxies and their central BHs.

The simplest and most elegant argument relating the properties of QSOs to BHs in nearby galaxies is due to Sołtan (1982), who pointed out that the luminosity function of QSOs as a function of redshift traces the accretion history of these BHs. The BH mass density due to bright QSO phases can be estimated directly from the local energy density in QSO photons and an assumed mass-to-energy conversion efficiency (Sołtan 1982; Chokshi & Turner 1992). Comparison of this density with the local BH mass density provides considerable insight into the formation and growth of massive BHs (e.g. Haehnelt & Kauffmann 2001 and references therein).

In § 2, we provide new estimates of the local BH mass density and the BH mass density accreted during bright QSO phases. There are several reasons why a re-analysis of this classic problem is timely. (i) The past several years have seen dramatic improvements in optical QSO surveys, both in terms of the total numbers of QSOs and the parameter ranges (i.e., luminosities and redshifts) over which the QSO luminosity function is reliably determined. For example, the 2dF QSO redshift survey (Boyle et al., 2000) and the Large Bright QSO Survey (Hewett, Foltz & Chaffee, 1995) have found $`6000`$ QSOs with absolute magnitude $`26<M_B<23`$ ($`\mathrm{\Omega }_\mathrm{m}=1,\mathrm{\Omega }_\mathrm{\Lambda }=0,h=0.5`$), at redshifts $`0.35<z<2.3`$; the Sloan Digital Sky Survey (SDSS) has found several hundred QSOs with redshifts $`3.6<z<5.0`$ (e.g. Fan et al. 2001; see also other surveys at high redshifts, e.g. Schmidt, Schneider & Gunn 1995 etc.) (ii) The SDSS has provided a sample of $`9000`$ nearby early-type galaxies (elliptical galaxies and S0 galaxies; Bernardi et al. 2001) with accurate luminosities and velocity dispersions. Using the tight correlation between BH mass and galactic velocity dispersion, we can estimate the BH masses in these galaxies and thereby obtain the mass function of BHs at $`z=0`$. (iii) A large number of obscured active galactic nuclei (AGNs) (e.g. $``$4–10 times the number of unobscured AGNs) is required to explain the observed X-ray background spectrum (e.g. Fabian & Iwasawa 1999; Gilli, Salvati & Hasinger 2001). The BH mass density inferred from the X-ray background by Fabian & Iwasawa (1999) is larger than the one derived from optically bright QSOs by Chokshi & Turner (1992) by a factor of 3–4. BHs in local galaxies are the remnants of both obscured and unobscured AGNs. Thus, an accurate estimate of the local BH mass density and the accreted mass density from optically bright QSOs can also constrain the fraction of obscured AGNs for comparison to models of the X-ray background.

If BH growth occurs mainly during optically bright QSO phases, then not only should the BH mass density in local galaxies be consistent with the total energy density in QSO light, but also the BH mass distribution in local galaxies should be consistent with the QSO luminosity function (e.g. Small & Blandford 1992; Ciotti, Haiman & Ostriker 2001; Marconi & Salvati 2001). In § 3, by studying the continuity equation for the BH mass distribution, we test this hypothesis and give more general observational constraints (e.g. on the mass-to-energy conversion efficiency and lifetime of QSOs) on the growth of massive BHs in galactic nuclei. The effects of BH mergers (caused by galaxy mergers) on the BH mass distribution are also considered. Since the physics of BH mergers is poorly understood, we consider two extreme cases: one is the classical case (mergers do not emit gravitational radiation and the total BH mass is conserved during BH mergers), which is implicitly assumed in Sołtan’s argument; the other is the adiabatic case (total BH area and entropy is conserved during BH mergers and mass is radiated away as gravitational radiation; see also Ciotti & van Albada 2001). Our conclusions are summarized in § 4.

In this paper, the Hubble constant is written as $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$; and if not otherwise specified, the cosmological model used is $`(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_\mathrm{\Lambda },h)=(0.3,0.7,0.65)`$ (Wang et al., 2000).

## 2 Total BH mass density

### 2.1 Early-type galaxies in the SDSS

The SDSS will image $`10^4`$ square degrees of the sky in five bands $`(u,g,r,i,z)`$ with central wavelengths (3560, 4680, 6180, 7500, 8870)Å, and take spectra of $`10^6`$ galaxies and $`10^5`$ QSOs (York et al., 2000). Among the spectra to be taken by SDSS, there will be roughly $`2\times 10^5`$ spectra of early-type galaxies. A sample of nearly 9000 of these with well-measured velocity dispersions ($`\delta \mathrm{log}\sigma `$0.02–$`0.06\mathrm{dex}`$ with a median value $`0.03\mathrm{dex}`$), in the redshift range $`0.01z0.3`$, has been selected from early SDSS observations by Bernardi et al. (2001). The joint distribution of luminosities, sizes and velocity dispersions for this sample is well described by a tri-variate Gaussian distribution in the variables $`M=2.5\mathrm{log}L`$, $`R=\mathrm{log}R_o`$, and $`V=\mathrm{log}\sigma `$, where $`L`$ is the total luminosity, $`R_o`$ is the effective radius and $`\sigma `$ is the line-of-sight velocity dispersion within a circular aperture of radius $`R_o/8`$. A maximum-likelihood fit to the sample gives the mean values of the above variables, $`M=19.65+5\mathrm{log}h1.15z`$ (the $`1.15z`$ term represents evolution), $`R=0.36`$, and $`V=2.20`$ in the $`g^{}`$ band (here $`10^R`$ and $`10^V`$ are in units of $`h^1\mathrm{kpc}`$ and $`\mathrm{km}\mathrm{s}^1`$, respectively; $`g^{}`$ is used rather than $`g`$, similarly for other bands, because the photometric calibration is preliminary; for details, see Stoughton et al. 2002) and their standard deviations $`\sigma _M=0.84`$, $`\sigma _R=0.25`$ and $`\sigma _V=0.11`$ (cf. table 2 in Bernardi et al. 2001). By extrapolating this luminosity distribution at both faint and bright ends, the average comoving number density of early-type galaxies is found to be

$$\varphi _{}=(5.8\pm 0.3)\times 10^3h^3\mathrm{Mpc}^3.$$
(1)
For comparison to other studies, we must convert magnitudes in the SDSS $`g^{}`$ band to the Johnson-Morgan $`B`$ band. We take the color transformation to be:

$$B=g^{}+0.47(g^{}r^{})+0.17,$$
(2)
which is derived from the preliminary SDSS calibration (Smith et al., 2001) of the Landolt stars (Landolt, 1992a, b). The standard deviation in this relation is $`0.03\mathrm{mag}`$. The mean color of the early-type galaxies in Bernardi et al.’s SDSS sample is $`g^{}r^{}0.740.30z`$ with standard deviation $`0.06\mathrm{mag}`$ in the colors of individual galaxies (cf. table 6 in Bernardi et al. 2001).

### 2.2 Local BH mass density

The SDSS data provide us with the distributions of velocity dispersion and luminosity in nearby early-type galaxies. We denote the luminosity function and the velocity-dispersion distribution of early-type galaxies as $`n_L^{\mathrm{early}}(L,z)`$ and $`n_\sigma ^{\mathrm{early}}(\sigma ,z)`$, where $`n_L^{\mathrm{early}}(L,z)\mathrm{d}L`$ \[or $`n_\sigma ^{\mathrm{early}}(\sigma ,z)\mathrm{d}\sigma `$\] represents the comoving number density of early-type galaxies in the range $`LL+\mathrm{d}L`$ (or $`\sigma \sigma +\mathrm{d}\sigma `$) at redshift $`z`$. The distributions $`n_L^{\mathrm{early}}(L,z)`$ and $`n_\sigma ^{\mathrm{early}}(\sigma ,z)`$ are obtained by integrating the tri-variate distribution in $`M=2.5\mathrm{log}L`$, $`R=\mathrm{log}R_o`$ and $`V=\mathrm{log}\sigma `$ over two of these variables. Thus, the average number density of these galaxies is (eq. 1)

$$\varphi _{}=_0^{\mathrm{}}n_L^{\mathrm{early}}(L,z=0)dL=_0^{\mathrm{}}n_\sigma ^{\mathrm{early}}(\sigma ,z=0)d\sigma ,$$
(3)
and the luminosity density is

$$j=_0^{\mathrm{}}Ln_L^{\mathrm{early}}(L,z=0)𝑑L.$$
(4)
The luminosity density derived from the Bernardi et al. sample is $`j_g^{}=6.1\times 10^7h\mathrm{L}_{\mathrm{g}^{},}\mathrm{Mpc}^3`$. For comparison, the total luminosity density derived from SDSS is $`j_g^{}=(2.8\pm 0.4)\times 10^8h\mathrm{L}_{\mathrm{g}^{},}\mathrm{Mpc}^3`$ (Blanton et al. 2001). This total luminosity density is a factor of 1.5 to 2 times higher than earlier estimates, but Blanton et al. argue that the difference arises because earlier surveys neglect low surface-brightness regions that are included in SDSS. The luminosity density in the Blanton et al. sample can be compared directly to the Bernardi et al. sample, as follows. Strateva et al. (2001) show that early- and late-type galaxies in SDSS can be separated by a cut in concentration index at $`C=2.6`$; applying this cut to the Blanton et al. sample we derive $`j_g^{}=7.4\times 10^7h\mathrm{L}_{\mathrm{g}^{},}\mathrm{Mpc}^3`$ for early-type galaxies. This result differs by only 20% from the estimate of Bernardi et al., which uses several parameters including the concentration index to separate early- from late-type galaxies. Thus the normalization of the sample isolated by Bernardi et al. appears to be consistent with other early-type galaxy samples in SDSS.

The relation between BH mass and galactic velocity dispersion is (Tremaine et al., 2002):

$$M_{}=(1.5\pm 0.2)\times 10^8\mathrm{M}_{}\left(\frac{\sigma _\mathrm{e}}{200\mathrm{km}\mathrm{s}^1}\right)^{\alpha _1},$$
(5)
where $`\alpha _1=4.02\pm 0.32`$, $`M_{}`$ is the BH mass, and $`\sigma _\mathrm{e}`$ is the luminosity-weighted line-of-sight velocity dispersion within a slit extending to the effective radius. Note that $`\sigma _\mathrm{e}`$ in equation (5) is the velocity dispersion within a slit extending to the effective radius $`R_o`$ (Gebhardt et al., 2000) while $`\sigma `$ in the SDSS is the velocity dispersion within a circular aperture extending to $`R_o/8`$; however, replacing $`\sigma _\mathrm{e}`$ with $`\sigma `$ in equation (5) will not cause much difference as the two definitions should give very similar results (Tremaine et al. 2002). Merritt & Ferrarese (2001a) give an alternative version of the correlation (5) that is based on the velocity dispersion within $`R_o/8`$, and find a steeper slope, $`\alpha _1=4.72`$. The reasons for the difference in slopes are discussed by Tremaine et al. (2002). We have performed all of the calculations below using both versions of the correlation and the differences are negligible.

The mean relation between the $`B`$-band luminosity of the hot stellar component and the BH mass is (Kormendy & Gebhardt, 2001)

$$M_{}=0.7\times 10^8\mathrm{M}_{}\left(\frac{L_{B,\mathrm{hot}}}{10^{10}L_B_{}}\right)^{\alpha _2},$$
(6)
where $`\alpha _2=1.08`$ (Kormendy & Gebhardt, 2001). The scatter in relation (6) is substantially larger than the scatter in relation (5). The statistical uncertainties in $`\alpha _2`$ and in the normalization are about $`\pm 0.15`$–0.2 and 30%, respectively. A similar relation to equation (6) was derived by Magorrian et al. (1998), but with normalization larger by a factor of 3–4, presumably because of the limitations of ground-based spectroscopy and isotropic dynamical models.

Most of the galaxies in Tremaine et al. (2002) and Kormendy & Gebhardt (2001) have distances derived from Tonry et al. (2001), which implies a Hubble constant $`h=0.73\pm 0.04\pm 0.11`$ (Blakeslee et al., 2002), and the normalizations in equations (5) and (6) have been adjusted to our assumed Hubble constant $`h=0.65`$.

Because the mass–dispersion relation is the tighter of the two relations, we will first use equation (5) to estimate the local BH mass density. Thus, with the distribution of velocity dispersion $`n_\sigma ^{\mathrm{early}}(\sigma ,z=0)`$ determined by Bernardi et al. (2001), we have the local BH mass density in early-type galaxies:

$`\rho _{}^{\mathrm{early}}(z=0)`$ $`=`$ $`{\displaystyle M_{}n_M_{}^{\mathrm{early}}(M_{},z=0)dM_{}}={\displaystyle M_{}n_\sigma ^{\mathrm{early}}(\sigma ,z=0)\frac{\mathrm{d}\sigma }{dM_{}}dM_{}}`$ (7)
$`=`$ $`(1.6\pm 0.2)\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3,`$
where $`\sigma `$ has been equated to $`\sigma _\mathrm{e}`$ and $`\mathrm{d}\sigma /\mathrm{d}M_{}`$ is obtained from equation (5). The error estimate represents only the effects of the uncertainties in equations (1) and (5).

We may also use the mass–luminosity relation (6) to estimate $`\rho _{}^{\mathrm{early}}`$. We note first that the luminosity provided in Bernardi et al. (2001) is the total luminosity. Early-type galaxies include both elliptical galaxies and S0 galaxies. For S0 galaxies, the bulge and disk both contribute to the total luminosity. Therefore, when applying equation (6), we need to correct the total luminosity of an S0 galaxy to its bulge luminosity. Since the BH mass is approximately proportional to the luminosity of the hot stellar component (i.e., $`\alpha _2=1.08`$ in eq. 6 is approximately 1), we may simply correct the total luminosity density of early-type galaxies in Bernardi et al. (2001) by the ratio of the luminosity density of hot components to the total luminosity density. Fukugita, Hogan & Peebles (1998) estimated the relative luminosity densities $`f_i`$ (in the $`B`$ band) for elliptical galaxies, S0 galaxies, spiral and irregular galaxies, and hot components. They found $`(f_\mathrm{E},f_{\mathrm{S0}},f_{\mathrm{Spiral}+\mathrm{Irr}})=(0.12,0.23,0.65)`$ and $`f_{\mathrm{hot}}=0.385`$. They also estimate that in S0 galaxies, the bulge luminosity is a fraction $`f_{\mathrm{bulge},\mathrm{S0}}/f_{\mathrm{S0}}=0.64`$ of the total luminosity. Thus, the total luminosity in hot components of early-type galaxies occupies a fraction

$$\frac{f_\mathrm{E}+f_{\mathrm{bulge},\mathrm{S0}}}{f_\mathrm{E}+f_{\mathrm{S0}}}0.76$$
(8)
of the total luminosity in early-type galaxies. The local BH mass density in early-type galaxies is therefore given by:

$`\rho _{,L}^{\mathrm{early}}(z=0)`$ $`=`$ $`{\displaystyle M_{}n_M_{}^{\mathrm{early}}(M_{},z=0)dM_{}}0.76{\displaystyle M_{}n_L^{\mathrm{early}}(L,z=0)\frac{\mathrm{d}L}{\mathrm{d}M_{}}dM_{}}`$ (9)
$`=`$ $`2.0\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3,`$
where the color transformation between $`g^{}`$ and $`B`$ (eq. 2) has been used.

Massive BHs exist not only in early-type galaxies but also in spiral bulges. Hence, we must augment the estimates (7) and (9) by the contribution from BHs in spiral bulges. As above, we assume that the luminosity fraction is approximately equal to the BH mass fraction, and thus the total BH mass density is $`f_{\mathrm{hot}}/(f_\mathrm{E}+f_{\mathrm{bulge},\mathrm{S0}})1.44`$ times the BH mass density in early-type galaxies, i.e., from equation (7)

$$\rho _{}(z=0)1.44\rho _{}^{\mathrm{early}}=(2.3\pm 0.3)\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3,$$
(10)
and from equation (9)

$$\rho _{,L}(z=0)1.44\rho _{,L}^{\mathrm{early}}2.9\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3.$$
(11)
One additional correction to these estimates must be considered. In deriving equations (7) and (9), we have ignored any intrinsic dispersion in the mass–dispersion relation and the mass–luminosity relations (eqs. 5 and 6). Note that equations (5) and (6) are fitted in the $`\mathrm{log}M_{}`$$`\mathrm{log}\sigma `$ and $`\mathrm{log}M_{}`$$`\mathrm{log}L_{B,\mathrm{hot}}`$ spaces. Let us assume that the distribution in $`\mathrm{log}M_{}`$ at a given value of $`\sigma `$ or $`L_{B,\mathrm{hot}}`$ is Gaussian, with mean given by the relations (5) or (6), and standard deviation $`\mathrm{\Delta }_{\mathrm{log}M_{}}`$ that is independent of $`\sigma `$ or $`L_{B,\mathrm{hot}}`$. Then our estimates of the mass density must be increased by a factor

$$\mathrm{exp}\left[\frac{1}{2}(\mathrm{\Delta }_{\mathrm{log}M_{}}\mathrm{ln}10)^2\right]=1+\frac{1}{2}(\mathrm{ln}10)^2\mathrm{\Delta }_{\mathrm{log}M_{}}^2+\text{O}(\mathrm{\Delta }_{\mathrm{log}M_{}}^4).$$
(12)
The intrinsic scatter $`\mathrm{\Delta }_{\mathrm{log}M_{}}`$ in the mass-dispersion relation is small enough that it is difficult to distinguish from measurement errors, but is less than $`0.27\mathrm{dex}`$ (Tremaine et al., 2002). The corresponding correction factor is therefore between 1 and 1.2, and we shall adopt 1.1. The intrinsic dispersion in the mass–luminosity relation is $`\mathrm{\Delta }_{\mathrm{log}M_{}}0.5\mathrm{dex}`$ (note that only the intrinsic dispersion, not measurement error, is included in the definition of $`\mathrm{\Delta }_{\mathrm{log}M_{}}`$), so the corresponding correction is a factor of 2.0. Thus the density estimates (10) and (11) become

$$\rho _{}(z=0)=(2.5\pm 0.4)\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3,$$
(13)
and

$$\rho _{,L}(z=0)=5.8\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3.$$
(14)
The error estimate in equation (13) represents the combination of the $``$10% uncertainty in the correction factor 1.1 due to intrinsic scatter, and the error estimate in equation (10). The error in equation (14) is difficult to estimate but is certainly much larger than the error in equation (13). Earlier estimates include $`\rho _{,L}(z=0)=3.5\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3`$ (Salucci et al., 1999), $`\rho _{,L}(z=0)=3.7\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3`$ (Merritt & Ferrarese, 2001b), and $`\rho _{,L}(z=0)=(4\pm 2)\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3`$ (Marconi & Salvati, 2001).

The estimates (13) and (14) differ by a factor of more than two. We believe that equation (13) is substantially more reliable than equation (14), since the correlation between $`M_{}`$ and $`\sigma `$ is tighter and the result does not depend on the uncertain bulge-disk decomposition. Nevertheless, the discrepancy between the two equations is worth investigating. The problem cannot arise from different estimates of the local galaxy density, since both results employ the same galaxy sample (Bernardi et al., 2001); nor can it plausibly arise from errors in the BH mass determinations, since the mass-dispersion and mass-luminosity relation are based on the same set of mass estimates (see Tremaine et al. 2002 and Kormendy & Gebhardt 2001, respectively). To isolate the problem, it is useful to examine the Faber-Jackson relation between dispersion and luminosity. The mean relation derived from a large sample of early-type galaxies (Prugniel & Simien, 1996) by Forbes & Ponman (1999) is

$$\mathrm{log}(\sigma /\mathrm{km}\mathrm{s}^1)=2.250.102(M_B+20),$$
(15)
while the mean relation derived from the Bernardi et al. (2001) sample is

$$\mathrm{log}(\sigma /\mathrm{km}\mathrm{s}^1)=2.180.100(M_B+20).$$
(16)
Thus the mean absolute magnitude at a given dispersion in the Bernardi et al. (2001) sample is displaced by $`\mathrm{\Delta }_1M_B=0.7`$ from equation (15). On the other hand, the mean residual in absolute magnitude from the relation (15) for the 31 early-type galaxies with BH mass determinations in Tremaine et al. (2002) is $`\mathrm{\Delta }_2M_B=+0.54`$. A component of the latter difference arises because equation (15) includes both bulge and disk luminosity in S0 galaxies, while the luminosities quoted in Tremaine et al. (2002) refer only to the bulges in S0 galaxies. Correcting crudely for this using equation (8) reduces $`\mathrm{\Delta }_2M_B`$ to +0.24. Together, $`\mathrm{\Delta }_1M_B`$ and $`\mathrm{\Delta }_2M_B`$ imply that the mean absolute magnitude at a given dispersion is about 1 magnitude brighter in the $`10^4`$ SDSS galaxies used to determine the velocity-dispersion and luminosity distribution of nearby galaxies than in the $`30`$ galaxies used to determine the mass-luminosity and mass-dispersion relation for BHs. This difference is sufficient to explain the difference in BH density estimates in equations (13) and (14). What is less clear is the origin of this difference. Possible explanations include selection effects in either the Bernardi et al. (2001) sample or the sample of galaxies with BH masses, or systematic differences in the magnitude or dispersion measurements between SDSS and pointed galaxy surveys.

We close this section with comments on the contribution to the estimated BH density from the faint end of the galaxy luminosity function. The SDSS galaxy sample is almost complete at the bright end, but the sample is cut off at the faint end ($`M_g^{}5\mathrm{log}h19.0`$, cf. Fig. 6 in Bernardi et al. 2001). In our calculation we have extrapolated Bernardi et al.’s distribution of magnitude (or velocity dispersion) to arbitrarily faint galaxies. This uncertain extrapolation should not significantly affect our estimate of BH mass density, as seen from the following arguments. (i) If we assume that the distribution of absolute magnitude obtained in the SDSS is constant for all magnitudes fainter than $`M_g^{}=19.65+5\mathrm{log}h`$ (rather than obeying a Gaussian distribution in $`M_g^{}`$), the total luminosity density is increased by less than 15%. (ii) The 2dF Galaxy Survey extends to considerably fainter absolute magnitudes ($`M_{\mathrm{b}_\mathrm{J}}14+5\mathrm{log}h`$) \[$`M_{\mathrm{b}_\mathrm{J}}=M_g_{}+0.16+0.15(g^{}r^{})`$ in Norberg et al. 2001\] than the SDSS sample. The type 1 galaxies in the 2dF Galaxy Survey include all the early-type galaxies and some spiral galaxies. The difference between the total luminosity density obtained by assuming that the distribution of absolute magnitude of the type 1 galaxies is constant for all magnitudes fainter than $`M_{\mathrm{b}_\mathrm{J}}=M_g^{}+0.16+0.15(g^{}r^{})`$ and that obtained by extrapolating the distribution of absolute magnitude of the type 1 galaxies (cf. Madgwick et al. 2001) is less than 4%. In addition, in § 3 below (see also Figs. 2 and 3), we will focus on comparing the distribution of local high-mass BHs ($`10^8\mathrm{M}_{}`$) with the luminosity function of luminous QSOs, where we need not consider uncertainties of local BH mass function at the low-mass ends.

### 2.3 BH mass density accreted during bright QSO phases

Most astronomers believe that QSOs are powered by gas accretion onto BHs, and in this case the luminosity function of QSOs as a function of redshift reflects the gas accretion history of local remnant BHs (Sołtan, 1982; Chokshi & Turner, 1992; Small & Blandford, 1992). We denote the QSO bolometric luminosity produced by a mass accretion rate $`\dot{M}_{\mathrm{acc}}`$ as $`L_{\mathrm{bol}}=ϵ\dot{M}_{\mathrm{acc}}c^2=ϵ\dot{M}_{}c^2/(1ϵ)`$, where $`ϵ`$ is the mass-to-energy conversion efficiency (up to 31% for thin-disk accretion onto a maximally rotating Kerr BH; Thorne 1974), $`\dot{M}_{}=\dot{M}_{\mathrm{acc}}(1ϵ)`$ is the growth rate of BH mass, and $`c`$ is the speed of light. The comoving BH mass density accreted during bright QSO phases is given by

$$\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z)=_z^{\mathrm{}}\frac{\mathrm{d}t}{\mathrm{d}z}dz_0^{\mathrm{}}\frac{(1ϵ)L_{\mathrm{bol}}}{ϵc^2}\mathrm{\Psi }(L,z)dL,$$
(17)
where $`\mathrm{\Psi }(L,z)`$ is the luminosity function of QSOs, defined so that $`\mathrm{\Psi }(L,z)dL`$ gives the comoving number density of QSOs at redshift $`z`$ ($`L`$ can be either the bolometric luminosity or the luminosity at a given band). In the $`B`$ band, the bolometric correction $`C_B`$, defined by $`L_{\mathrm{bol}}=C_B\nu _BL_\nu (B)`$, is about 11.8 (Elvis et al., 1994). Here $`\nu _BL_\nu (B)`$ is the energy radiated at the central frequency of the $`B`$ band per unit time and logarithmic interval of frequency.

At $`z3`$, the QSO luminosity function is often fitted with a double power law:

$$\mathrm{\Psi }_M(M_B,z)=\frac{\mathrm{\Psi }_M^{}}{10^{0.4(\beta _1+1)[M_BM_B^{}(z)]}+10^{0.4(\beta _2+1)[M_BM_B^{}(z)]}},$$
(18)
where $`\mathrm{\Psi }_M(M_B,z)\mathrm{d}M_B`$ is the comoving number density of QSOs with absolute magnitude in the range \[$`M_B,M_B+\mathrm{d}M_B`$\] and $`\mathrm{\Psi }_M(M_B,z)=\mathrm{\Psi }(L_B,z)|\mathrm{d}L_B/\mathrm{d}M_B|=0.92L_B\mathrm{\Psi }(L_B,z)`$. Boyle et al. (2000) use this functional form to fit the data sets from the 2dF QSO Redshift Survey (Boyle et al., 2000) and the Large Bright QSO Survey (Hewett, Foltz & Chaffee, 1995), which contain over 6000 QSOs. For example, in our standard model, $`(\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_\mathrm{\Lambda })=(0.3,0.7)`$, the luminosity function of QSOs with absolute magnitudes $`26<M_B<23`$ and redshifts $`0.35<z<2.3`$ is described by the following parameters:

$`\mathrm{\Psi }_M^{}=2.9\times 10^5h^3\mathrm{Mpc}^3\mathrm{mag}^1,`$ (19)
$`M_B^{}(z)=M_B^{}(0)2.5(k_1z+k_2z^2),`$ (20)
$`M_B^{}(0)=21.14+5\mathrm{log}h,k_1=1.36,k_2=0.27,`$ (21)
$`\beta _1=1.58\mathrm{and}\beta _2=3.41.`$ (22)
The quadratic dependence of the characteristic magnitude $`M_B^{}(z)`$ on $`z`$ (eqs. 20 and 21) shows an increasing characteristic luminosity with increasing redshift at low redshift ($`z2.5`$) and a decline of the characteristic luminosity at higher redshift, which is suggested by observations (e.g. Shaver et al. 1996), although the luminosity function is not yet accurate enough to confirm the decline at $`z>2.5`$. The QSO luminosity function over the range $`3.6<z<5.0`$ provided in Fan et al. (2001) gives a flatter bright-end slope ($`\beta _2=2.5`$) than equation (18); however in our calculations below, we will simply extrapolate equation (18) to high redshift ($`z>3`$) because the detailed QSO luminosity function at $`z>3`$ does not affect our final results much (cf. Figure 1 below or the discussion in Chokshi & Turner 1992; if we use the luminosity function in Fan et al. 2001, our results change by less than one percent).

Using equation (17) we may now obtain the accretion history of the BH mass density arising from optically bright QSO phases, $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z)`$. The result is shown in Figure 1, where the mass-to-energy conversion efficiency $`ϵ`$ is assumed constant. About 56% of the present density $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)`$ is accreted when $`z<2`$ and 90% when $`z<3`$, consistent with our claim that $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)`$ is insensitive to uncertainties in the QSO luminosity function at high redshift. The accreted BH mass density during optically bright QSO phases is found to be

$$\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)=2.1\times 10^5(C_B/11.8)[0.1(1ϵ)/ϵ]\mathrm{M}_{}\mathrm{Mpc}^3,$$
(23)
which is independent of the Hubble constant and remarkably close to the value obtained by Chokshi & Turner (1992) using similar methods, $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)=2.2\times 10^5(C_B/16.5)(0.1/ϵ)\mathrm{M}_{}\mathrm{Mpc}^3`$. If the QSOs have a mass-to-energy conversion efficiency $`ϵ0.1`$, $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)`$ would be close to the local BH mass density $`\rho _{}`$ (eq. 13), which implies that the local BH mass density comes mainly from accretion during optically bright QSO phases.

The accreted BH mass density can also be inferred from the X-ray and infrared backgrounds, for example, $`\rho _{}^\mathrm{X}`$(6–9)$`\times 10^5(0.1/ϵ)\mathrm{M}_{}\mathrm{Mpc}^3`$ in Fabian & Iwasawa (1999), $`\rho _{}^\mathrm{X}`$ (7.5–16.8)$`\times 10^5(0.1/ϵ)\mathrm{M}_{}\mathrm{Mpc}^3`$ in Elvis, Risaliti & Zamorani (2002), $`\rho _{}^{\mathrm{IR}}7.5\times 10^5(0.1/ϵ)`$ $`\mathrm{M}_{}`$$`\mathrm{Mpc}^3`$ in Haehnelt & Kauffmann (2001). Multi-waveband (including optical and far infrared) observations of Chandra hard X-ray sources give $`\rho _{}^{\mathrm{HX}+\mathrm{OPT}+\mathrm{FIR}}9\times 10^5(0.1/ϵ)`$ $`\mathrm{M}_{}\mathrm{Mpc}^3`$ in Barger et al. (2001). Fabian & Iwasawa (1999) conclude that 85% of the emitted energy from AGNs has been absorbed, so that a large fraction of the local BH mass density must be due to obscured QSOs, which do not contribute to the estimate in equation (23). Thus, their estimate $`\rho _{}^\mathrm{X}`$ exceeds $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)`$ by a factor 3–4, and then the comparison of $`\rho _{}^\mathrm{X}`$ with the local BH mass density (eq. 13) implies the efficiency $`ϵ`$0.3–0.5, at or beyond the upper limit of plausible accretion processes, $`ϵ=0.31`$ (Thorne, 1974). Elvis, Risaliti & Zamorani (2002) reached a similar conclusion, $`ϵ0.15`$, by comparing $`\rho _{}^\mathrm{X}`$ with the local BH mass density \[they used an early estimate of the local BH mass density $`3.5\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3`$ from Salucci et al. (1999), a factor of 1.4 larger than our estimate in eq. (13)\]. However, there exist many uncertainties in the estimation of the BH mass density from the X-ray or infrared background since the population of AGNs contributing to those backgrounds is not yet clearly understood, and is not necessarily the same as the optically bright QSOs.

If the local BH mass density comes mainly from optically bright QSOs, not only the local BH mass density but also the local BH mass distribution should be consistent with the accretion history determined from QSOs. The latter comparison is the subject of the next section.

## 3 Some clues from the BH mass distribution

In this section, we employ a continuity equation for the evolution of the BH distribution to relate the local BH mass function and the QSO luminosity function.

### 3.1 Expectations from the continuity equation for BH distributions

We define $`S`$ to be any physical variable associated with a BH (e.g. mass, entropy etc.) and $`n_S(S,t)`$ so that $`n_S(S,t)\mathrm{d}S`$ is the number density of BHs with $`SS+\mathrm{d}S`$ at time $`t`$. The evolution of the distribution $`n_S(S,t)`$ is determined by both BH coalescence (caused by galaxy mergers) and other physical processes (e.g. gas accretion). We can use the following continuity equation (cf. Small & Blandford 1992) to describe the evolution of the $`S`$-distribution:

$$\frac{n_S(S,t)}{t}+\frac{[n_S(S,t)\dot{S}]}{S}=\gamma _{\mathrm{merge}}(S,t),$$
(24)
where $`\gamma _{\mathrm{merge}}(S,t)`$ represents the variation of the BH distribution caused by mergers, and $`\dot{S}(S,t)`$ is the mean rate of change of $`S`$ due to physical processes that conserve the number of BHs. In equation (24), we ignore the generation of seed BHs, which are usually assumed to have small masses (e.g. $`10^6\mathrm{M}_{}`$). Since the process of BH mergers is not understood, we shall consider two extreme cases: one is the classical case, in which gravitational radiation can be neglected and the total BH mass is conserved during a merger; and the other is the adiabatic case, in which gravitational radiation cannot be ignored and the total BH area (or entropy) is conserved. We will consider the evolution of the BH mass distribution in the classical case (i.e., set $`S`$ in eq. 24 to be the reducible<sup>1</sup><sup>1</sup>1The dynamically measured BH mass $`M_{}`$ (in eqs. 5 and 6) is the reducible mass. For a BH with reducible mass $`M_{}`$ and dimensionless spin $`a`$ ($`0|a|1`$), the irreducible mass is given by $`M_{\mathrm{ir}}=M_{}(1+\sqrt{1a^2})^{1/2}/\sqrt{2}`$ (Misner, Thorne & Wheeler, 1973). Thus, $`M_{}/\sqrt{2}M_{\mathrm{ir}}M_{}`$. BH mass $`M_{}`$), and the evolution of the area distribution in the adiabatic case (i.e., set $`S`$ in eq. 24 to be the BH area $`A_{}`$). In either case, mergers conserve the total value of $`S`$, so we may consider a time interval $`\mathrm{d}t`$ in which two BHs parameterized by $`S_1`$, $`S_2`$ merge to form a single BH parameterized by $`S_1+S_2`$, i.e.,

$$\gamma _{\mathrm{merge}}(S^{},t)\mathrm{d}t=\delta (S^{}S_1S_2)\delta (S^{}S_1)\delta (S^{}S_2).$$
(25)
Thus

$$_0^{\mathrm{}}S\gamma _{\mathrm{merge}}(S,t)dS=0.$$
(26)
If $`\gamma _{\mathrm{merge}}(S,t)`$ were known, we could integrate equation (24) to obtain $`n_S(S,t)`$, which could be compared to the local BH mass distribution. Unfortunately the BH merger rate is very uncertain. We therefore ask a more restricted question: can we find some function of the mass distribution that can only increase during mergers, so that even without knowing the merger rate we can establish some inequalities on the current BH mass distribution and QSO luminosity function? For example, we would like to choose a function $`f(S,S^{})`$ so that

$`\mathrm{\Gamma }_{\mathrm{merge}}(S,t){\displaystyle _0^{\mathrm{}}}f(S,S^{})\gamma _{\mathrm{merge}}(S^{},t)dS^{}`$ (27)
is positive-definite. By inserting equation (25) in equation (27), we have

$$\mathrm{\Gamma }_{\mathrm{merge}}(S,t)\mathrm{d}t=f(S,S_1+S_2)f(S,S_1)f(S,S_2).$$
(28)
A sufficient condition that equation (28) is positive definite for $`S_1,S_2>0`$ is that $`f(S,S^{})/S^{}`$ is a monotonically increasing function of $`S^{}`$ (Hardy, Littlewood & P$`\stackrel{´}{\mathrm{o}}`$lya, 1934). We shall choose a specific function $`f(S,S^{})`$ below. By integrating equation (24) over time $`t`$ and setting the initial BH $`S`$-distribution $`n_S(S,0)=0`$, we may obtain the BH distribution $`n_S(S,t_0)`$ at the present time $`t_0`$ as follows:

$$n_S(S,t_0)=_0^{t_0}\frac{[n_S(S,t)\dot{S}]}{S}dt+_0^{t_0}\gamma _{\mathrm{merge}}(S,t)dt.$$
(29)
Now integrate equation (27) from $`t=0`$ to $`t_0`$, and use equation (29) to eliminate $`\gamma _{\mathrm{merge}}`$. Assuming $`f(S,0)=0`$, we have

$$G_{\mathrm{local}}(S,t_0)=G_{\mathrm{acc}}(S,t_0)+_0^{t_0}\mathrm{\Gamma }_{\mathrm{merge}}(S,t)dt,$$
(30)
where

$$G_{\mathrm{local}}(S,t_0)_0^{\mathrm{}}f(S,S^{})n_S(S^{},t_0)dS^{}$$
(31)
is related to the BH $`S`$-density in local galaxies, and

$$G_{\mathrm{acc}}(S,t_0)_0^{\mathrm{}}dS^{}_0^{t_0}n_S(S^{},t)\dot{S}^{}\frac{f(S,S^{})}{S^{}}dt$$
(32)
is related to the accretion history of local BHs. Since $`\mathrm{\Gamma }_{\mathrm{merge}}(S,t)`$ is non-negative, we have

$$G_{\mathrm{local}}(S,t_0)G_{\mathrm{acc}}(S,t_0).$$
(33)
Inequality (33) must hold for every $`f(S,S^{})`$ such that $`f(S,S^{})/S^{}`$ is monotonically increasing in $`S^{}`$. We have assumed that $`f(S,0)=0`$; it is also useful for $`f(S,S^{})`$ to be continuous, since its derivative appears in equation (32). A simple function satisfying these conditions, which we shall use for the remainder of this paper, is

$$f(S,S^{})=\{\begin{array}{cc}0\hfill & \text{for }S^{}<S\text{,}\hfill \\ S^{}S\hfill & \text{for }S^{}S\text{.}\hfill \end{array}$$
(34)
In this case,

$$G_{\mathrm{local}}(S,t_0)=_S^{\mathrm{}}(S^{}S)n_S(S^{},t_0)dS^{},$$
(35)
$$G_{\mathrm{acc}}(S,t_0)=_S^{\mathrm{}}dS^{}_0^{t_0}n_S(S^{},t)\dot{S}^{}dt.$$
(36)
With this definition, equation (26) implies that $`\mathrm{\Gamma }_{\mathrm{merge}}(0,t)=0`$, so

$$G_{\mathrm{local}}(0,t_0)=G_{\mathrm{acc}}(0,t_0).$$
(37)
Using this result we can define a normalized version of the fundamental inequality (33),

$$g_{\mathrm{local}}(S,t_0)\frac{G_{\mathrm{local}}(S,t_0)}{G_{\mathrm{local}}(0,t_0)}\frac{G_{\mathrm{acc}}(S,t_0)}{G_{\mathrm{acc}}(0,t_0)}g_{\mathrm{acc}}(S,t_0).$$
(38)
In the classical case where $`S=M_{}`$, equation (37) is identical to Sołtan’s argument (1982) relating the total local energy density in QSOs to the total BH mass density in nearby galaxies, while the inequality (38) provides a constraint on the normalized BH mass distribution that is independent of Sołtan’s. Inequalities (33) and (38) will be our principal tools to constrain the accretion history of local BHs.

#### 3.1.1 Relating the local BH mass function to the QSO luminosity function

To proceed further, we now assume that the luminosity $`L`$ during bright QSO phases is an (increasing) function of only the central BH mass \[$`\mathrm{d}L(M_{})/\mathrm{d}M_{}>0)`$\], which is plausible since usually we have $`LL_{\mathrm{bol}}`$ and, especially for luminous QSOs, $`L_{\mathrm{bol}}(M_{})`$ is usually assumed to be near the Eddington luminosity. In our two extreme cases, we have

* Classical case: $`SM_{}`$. We have

$$G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)=_{L(M_{})}^{\mathrm{}}dL^{}_0^{t_0}\mathrm{\Psi }(L^{},t)\frac{(1ϵ)L_{\mathrm{bol}}^{}}{ϵc^2}dt.$$
(39)
which represents the time-integrated mass density accreted onto BHs with masses larger than $`M_{}`$ during optically bright QSO phases. Equation (39) is independent of the duty cycle of bright QSO phases, so long as the luminosity during these phases is determined by $`M_{}`$. If the local BH mass density arises entirely by mergers and by accretion during optically bright QSO phases, then by combining inequality (38) with equation (39) we have

$$G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)=G_{\mathrm{acc}}(M_{},t_0),$$
(40)
$$g_{\mathrm{local}}(M_{},t_0)\frac{G_{\mathrm{local}}(M_{},t_0)}{G_{\mathrm{local}}(0,t_0)}\frac{G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)}{G_{\mathrm{acc}}^{\mathrm{QSO}}(0,t_0)}g_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0).$$
(41)
* Adiabatic case: $`SA_{}=16\pi G^2M_{\mathrm{ir}}^2/c^4`$, where $`G`$ is the gravitational constant and $`M_{\mathrm{ir}}`$ is the irreducible mass of a BH with reducible mass $`M_{}`$ (see footnote 1). We assume that the dimensionless spin of BHs in QSOs is in an equilibrium state (i.e., $`\dot{a}=0`$), which we denote by $`a_{\mathrm{QSO}}`$. For example, in the thin-disk accretion model for QSOs <sup>2</sup><sup>2</sup>2We assume that any QSO is due to accretion onto a single massive BH, rather than accretion onto two or more BHs in the galactic nucleus, although binary BHs are likely to be present in some galactic centers (Begelman, Blandford & Rees 1980; Yu 2002). If there exists more than one BH in a QSO and each BH is accreting with the Eddington luminosity, the mass $`M_{}`$ obtained from the total QSO luminosity will represent the total mass of the BHs, and our arguments in the classical case are unaffected. However, in the adiabatic case, the entropy or area increase associated with a given mass accretion rate depends on the number of BHs., the dimensionless spin $`a`$ increases from 0 to reach the equilibrium value $`a_{\mathrm{QSO}}=0.998`$ after the BH accretes $`1.5`$ times its initial mass (Thorne 1974). Thus, in the equilibrium spin state, we have the area (or entropy) variation rate of a BH as follows:

$$\dot{A}_{}=32\pi G^2M_{\mathrm{ir}}\dot{M}_{\mathrm{ir}}/c^4=16\pi G^2M_{}\dot{M}_{}\left(1+\sqrt{1a_{\mathrm{QSO}}^2}\right)/c^4(\dot{a}=0).$$
(42)
If the local BH entropy density arises entirely by mergers and by accretion during optically bright QSO phases, then equation (36) can be written as

$$G_{\mathrm{acc}}(A_{},t_0)=_{(M_{\mathrm{ir}})}^{\mathrm{}}dL^{}_0^{t_0}\mathrm{\Psi }(L^{},t)\dot{A}_{}(M_{}^{},a_{\mathrm{QSO}})dt,$$
(43)
where the function $`(M_{\mathrm{ir}})`$ is defined by $`(M_{\mathrm{ir}})L(M_{})`$. It proves useful to use the irreducible mass rather than the area as the independent variable, so we define

$$𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)G_{\mathrm{acc}}(A_{},t_0),\text{where}A_{}=16\pi G^2M_{\mathrm{ir}}^2/c^4.$$
(44)
In the adiabatic case, equation (35) may be written as

$$G_{\mathrm{local}}(A_{},t_0)=_{M_{}(A_{},a)}^{\mathrm{}}[A_{}(M_{}^{},a)A_{}]n_M_{}(M_{}^{},t_0)dM_{}^{}$$
(45)
where for the moment we assume that all the local BHs have the same spin $`a`$. Unfortunately, the BH spin parameter $`a`$ in dead quasars is unknown (in particular, there is no reason to suppose that it equals the equilibrium spin parameter $`a_{\mathrm{QSO}}`$ for accreting BHs). However, for a given value of the reducible mass $`M_{}`$, both the irreducible mass $`M_{\mathrm{ir}}`$ and the area $`A_{}`$ decrease with increasing $`a`$. Thus we have the inequality

$$G_{\mathrm{local}}(A_{},t_0)_{M_{}(A_{},a)}^{\mathrm{}}[A_{}(M_{}^{},a=0)A_{}]n_M_{}(M_{}^{},t_0)dM_{}^{}.$$
(46)
Once again introducing the irreducible mass rather than the area as the independent variable, we define

$$𝒢_{\mathrm{local}}^A_{}(M_{\mathrm{ir}},t_0)_{M_{\mathrm{ir}}}^{\mathrm{}}[A_{}(M_{}^{},a=0)A_{}(M_{\mathrm{ir}})]n_M_{}(M_{}^{},t_0)dM_{}^{}.$$
(47)
The integral on the right side of equation (46) differs from the integral on the right side of equation (47) only in the lower limit of the integrand. Since $`M_{\mathrm{ir}}M_{}`$ (see footnote 1), the second integral has a larger integration range and therefore is larger. Thus

$$G_{\mathrm{local}}(A_{},t_0)𝒢_{\mathrm{local}}^A_{}(M_{\mathrm{ir}},t_0).$$
(48)
The inequality (48) is an equality if and only if BHs in local galaxies are all Schwarzschild BHs ($`a=0`$). By using inequalities (33) and (48) and equations (43) and (44), we have finally:

$$𝒢_{\mathrm{local}}^A_{}(M_{\mathrm{ir}},t_0)𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0).$$
(49)
Note that inequality (49) still holds even if we discard the assumption made after equation (45) that all local BHs have the same spin parameter.

To estimate $`G_{\mathrm{local}}`$ and $`𝒢_{\mathrm{local}}^A_{}`$, we proceed as follows. High-luminosity galaxies are mainly early-type galaxies (see, for example, Figure 4.14 in Binney & Merrifield 1998 or Figure 7 in Madgwick et al. 2001); moreover, the mass of the central BH is approximately proportional to the luminosity of the hot stellar component (cf. eq. 6), which is less than $``$30% of the total luminosity for spiral galaxies (Fukugita, Hogan & Peebles, 1998); hence high-mass BHs ($`M_{}10^8\mathrm{M}_{}`$, corresponding to $`M_{B,\mathrm{hot}}<20`$) are mostly in early-type galaxies. The local BH mass distribution in early-type galaxies can be obtained from their velocity dispersions and the mass–dispersion correlation. By analogy with equations (35) and (45), we define

$$G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)_M_{}^{\mathrm{}}(M_{}^{}M_{})n_M_{}^{\mathrm{early}}(M_{}^{},t_0)dM_{}^{}$$
(50)
and

$$𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)_M_{}^{\mathrm{}}[A_{}(M_{}^{},a=0)A_{}(M_{\mathrm{ir}})]n_M_{}^{\mathrm{early}}(M_{}^{},t_0)dM_{}^{}.$$
(51)
For high-mass BHs ($`M_{}10^8\mathrm{M}_{}`$), we have $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)G_{\mathrm{local}}(M_{},t_0)`$ and $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{},t_0)`$ $`𝒢_{\mathrm{local}}^A_{}(M_{},t_0)`$. In addition, as seen from inequalities (33), (41) and (49), we have

$$G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0),$$
(52)
$$g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)\frac{G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)}{G_{\mathrm{local}}(0,t_0)}\frac{G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)}{G_{\mathrm{acc}}^{\mathrm{QSO}}(0,t_0)}g_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)(M_{}10^8\mathrm{M}_{}),$$
(53)
in the classical case (with approximate equality if and only if mergers are neglected) and

$$𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)(M_{\mathrm{ir}}10^8\mathrm{M}_{})$$
(54)
in the adiabatic case (with approximate equality if and only if mergers are neglected and BHs in local galaxies are all Schwarzschild BHs). The approximate nature of the inequalities arises solely from our restriction to early-type galaxies.

### 3.2 Comparison with observations

We now compare inequalities (52)–(54) with the observations. Throughout this section, all the variables related to the distribution of BH mass and entropy density, $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$, $`g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ and $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)`$, are obtained by using the mass–dispersion relation (eq. 5) and ignoring the small intrinsic scatter in this relation.

#### 3.2.1 The classical case

In the classical case, the total BH mass is conserved during mergers, which is implicitly assumed in Sołtan’s argument and many subsequent discussions of the growth of BHs. In § 2, we have obtained that the total BH mass density in local galaxies is consistent with that accreted during optically bright QSO phases, if the mass-to-energy efficiency $`ϵ0.1`$, i.e., $`G_{\mathrm{local}}(0,t_0)G_{\mathrm{acc}}^{\mathrm{QSO}}(0,t_0)`$ (cf. eqs. 13 and 23). In contrast, inequality (53) involves normalized quantities $`g(M_{},t_0)=G(M_{},t_0)/G(0,t_0)`$ and therefore provides a constraint independent of Sołtan’s argument. \[For example, note that the normalized quantity $`g_{\mathrm{acc}}(M_{},t_0)`$ in equation (53) is independent of the efficiency $`ϵ`$ if $`ϵ`$ is a constant.\]

Figure 2(a) shows $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ derived from the SDSS sample of early-type galaxies (dotted and dashed lines, eq. 50) and $`G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$ derived from optically bright QSOs (thick solid line, eq. 39). The QSOs are assumed to have the Eddington luminosity and the efficiency $`ϵ`$ is given by $`ϵ/(1ϵ)=0.1`$. Figure 2(b) shows the normalized quantity $`g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ derived from the SDSS sample of early-type galaxies (dotted and dashed lines, eq. 53) and $`g_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$ derived from optically bright QSOs (thick solid line, eq. 53) using the same assumptions. The dotted line in panel (a) \[or panel (b)\] is used for $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ \[or $`g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$\] in the region $`M_{}10^8\mathrm{M}_{}`$ where spiral bulges can make a substantial contribution, so that $`G_{\mathrm{local}}^{\mathrm{early}}`$ (or $`g_{\mathrm{local}}^{\mathrm{early}}`$) is only a lower bound to $`G_{\mathrm{local}}`$ (or $`g_{\mathrm{local}}`$). For $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ and $`g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$, we show separately the results obtained by using the mass–dispersion relations in Tremaine et al. (2002) (thick dotted and dashed line) and Merritt & Ferrarese (2001a) (thin dotted and dashed line). Figure 2 shows that for high BH masses ($`M_{}10^8\mathrm{M}_{}`$), $`G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$ \[or $`g_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$\] (thick solid line) is larger than $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ \[or $`g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$\] (dashed lines) by a factor of more than two, which is inconsistent with inequalities (52) (or 33) and (53).

This inconsistency is robust to a number of changes in our assumptions. For example:

* We have assumed that QSOs radiate at the Eddington limit; if they have sub-Eddington luminosity (e.g. one tenth of the Eddington luminosity in Ciotti, Haiman & Ostriker 2001), the solid lines in Figures 2(a) and (b) would shift to higher mass and the discrepancy would become worse.
* The QSO luminosity function is uncertain at the bright end, and the mass–dispersion relation is uncertain at the high-mass end. However, even if we cut off the QSO luminosity function at $`L_{\mathrm{bol}}L_{\mathrm{Edd}}(10^9\mathrm{M}_{})10^{47}\mathrm{erg}\mathrm{s}^1`$ ($`L_{\mathrm{Edd}}`$: Eddington luminosity), the new $`G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$ \[or $`g_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$\] (thin solid line in Figure 2) is still significantly higher than $`G_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$ \[or $`g_{\mathrm{local}}^{\mathrm{early}}(M_{},t_0)`$\] for $`10^8\mathrm{M}_{}M_{}10^9\mathrm{M}_{}`$.
* If there exist mechanisms other than accretion during bright QSO phases \[e.g. accretion via advection dominated accretion flow (ADAF), see discussion in Haehnelt, Natarajan & Rees (1998); or accretion of non-baryonic dark matter\] that contribute to the local BH mass density, then $`G_{\mathrm{acc}}^{\mathrm{QSO}}(M_{},t_0)`$ is only a lower limit to $`G_{\mathrm{acc}}(M_{},t)`$, and the inconsistency with inequality (33) for high-mass BHs would become worse.
* Study of the X-ray background suggests that a majority of the total BH mass density is contributed by obscured accretion (Fabian & Iwasawa, 1999). If the number ratio of obscured QSOs to optically bright QSOs is independent of BH mass (e.g., in the standard unification model for AGNs, obscuration is a purely geometrical effect and obscured fraction has no relation to central BH mass), $`g_{\mathrm{acc}}(M_{},t_0)`$ would be unchanged, and the inconsistency in Figure 2(b) would remain.

One possible resolution to the inconsistency shown in Figure 2 is that the mass-to-energy conversion efficiency $`ϵ`$ depends on QSO luminosity. In particular, the average efficiency in luminous QSOs (e.g. $`L_{\mathrm{bol}}10^{46}\mathrm{erg}\mathrm{s}^1`$, which corresponds to the Eddington luminosities of BHs with $`M_{}10^8\mathrm{M}_{}`$) must be $`0.2`$ to make the BH mass density accreted during optically bright QSO phases as large as that in local galaxies (cf. Fig. 2a). Note that the required average efficiency is close to the maximum efficiency ($`0.3`$ in Thorne 1974) allowed in thin-disk accretion models. This result suggests that other accretion mechanisms are less significant than accretion during optically bright QSO phases, and we may expect that obscured accretion (Fabian & Iwasawa, 1999) suggested by the X-ray background is not important for the growth of high-mass BHs, which is consistent with the scarcity of Type II QSOs.

If the growth of low-mass BHs also occurs mainly during optically bright QSO phases, then less luminous QSOs ($`L_{\mathrm{bol}}10^{46}\mathrm{erg}\mathrm{s}^1`$) must accrete with an efficiency less than $`0.1`$ to maintain the consistency between the local BH mass density and the BH mass density accreted during bright QSO phases (cf. eq. 37). If on the other hand, obscured accretion contributes significantly to the total BH mass density (Fabian & Iwasawa, 1999), less luminous QSOs may have efficiency $`0.1`$ or even greater. In this case, the obscured fraction must be larger for low-mass BHs than for high-mass BHs, which might suggest that obscured accretion occurs mainly for small BHs in young growing galaxies as described in the model by Fabian (1999).

Other possible resolutions to the inconsistency shown in Figure 2 include:

* Not all massive BHs may reside in galactic centers; for at least three reasons: (i) BHs may be left in the halo after galaxy mergers; however, Yu (2002) shows that the orbital decay time from dynamical friction for BHs in the galactic halo is generally much less than a Hubble time, so long as the orbiting BHs remain part of the tidally stripped remnant of their original host galaxy. (ii) Incomplete merging of binary BHs formed in mergers is expected to result in binary BHs with semimajor axes in the range $`10^3`$$`10\mathrm{pc}`$ depending on the velocity dispersion and shape of the host galaxy, and the masses of the BHs (Yu, 2002); however, in most cases the separation is so small that the binary should look like a single merged BH at current telescope resolutions. (iii) BHs can be ejected from galaxy centers through either three-body interactions with other massive BHs (e.g. Valtonen 1996) or gravitational radiation reaction during BH coalescence (e.g. Rees 2001).
* The bolometric correction that we have used, $`C_B=11.8`$, may be too large; however, the uncertainty in this estimate of the bolometric correction (Elvis et al., 1994) is only $`\pm 4.3`$, and the smallest bolometric correction in the Elvis et al. sample is 5.5, so it is unlikely that errors in the bolometric correction can lower the required efficiency by more than a factor of two.
* Luminous QSOs may accrete with super-Eddington luminosities (Begelman, 2001, 2002).

#### 3.2.2 The adiabatic case

BH mass may not be conserved during BH mergers, since mass can be radiated as gravitational waves. According to the second law of BH thermodynamics, BH mergers never decrease total BH entropy. Here, we study an extreme case in which the BH entropy is conserved during BH mergers.

Figure 3 shows the variables related to the BH entropy density as a function of BH irreducible mass that we introduced in § 3.1: $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)`$ obtained from the SDSS sample of early-type galaxies (dotted and dashed lines, see eq. 51) and $`𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)`$ obtained from optically bright QSOs (thick solid lines, see eq. 44). As usual QSOs are assumed to have the Eddington luminosity and $`ϵ`$ is given by $`ϵ/(1ϵ)=0.1`$ (eqs. 1822). The spin of BHs in QSOs is assumed to be in an equilibrium state, $`\dot{a}=0`$, $`aa_{\mathrm{QSO}}`$. The upper and lower thick solid lines represent the results obtained by assuming that all the BHs in QSOs are Schwarzschild ($`a_{\mathrm{QSO}}=0`$) and maximally rotating Kerr BHs ($`a_{\mathrm{QSO}}=1`$), respectively. The dotted lines are used for $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)`$ in the region $`M_{\mathrm{ir}}10^8\mathrm{M}_{}`$ where spiral bulges can make a substantial contribution, so that $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}`$ is only a lower bound to $`𝒢_{\mathrm{local}}^A_{}`$. For $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)`$, we show separately the results obtained by using the mass–dispersion relations in Tremaine et al. (2002) (thick dotted and dashed line) and in Merritt & Ferrarese (2001a) (thin dotted and dashed line). The functions $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)`$ obtained from these two relations differ by a factor of two, much more than the analogous curves in Figure 2. The difference is mainly caused by BHs with $`M_{}10^9\mathrm{M}_{}`$ because the BH area (or entropy) $`A_{}M_{\mathrm{ir}}^2`$ so the results are quite sensitive to the mass–dispersion relation at the high-mass end. Figure 3 shows that the thick solid lines \[i.e., $`𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)`$\] are much higher than the dashed lines \[i.e., $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(M_{\mathrm{ir}},t_0)`$\], no matter whether we choose $`a_{\mathrm{QSO}}=0`$ or $`a_{\mathrm{QSO}}=1`$ and for both mass–dispersion relations, which is inconsistent with inequality (54). This inconsistency may be partly caused by uncertainties in the QSO luminosity function at the bright end. To investigate the effects of uncertainties in the QSO luminosity function, we cut off the luminosity function at $`L_{\mathrm{bol}}L_{\mathrm{Edd}}(10^9\mathrm{M}_{})`$ (thin solid lines in Figure 3). For consistency, we must then use the Tremaine et al. mass–dispersion relation since it gives $`𝒢_{\mathrm{local}}^{A_{},\mathrm{early}}(10^9\mathrm{M}_{},t_0)0`$. We see from Figure 3 that the new $`𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)`$ for $`a_{\mathrm{QSO}}=0`$ (upper thin solid line) is still significantly higher than the thick dashed line, which is inconsistent with inequality (54). Considering that the maximum efficiency for thin-disk accretion onto Schwarzschild BHs is less than 0.1, a more realistic $`𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)`$ for $`a_{\mathrm{QSO}}=0`$ (i.e., $`ϵ0.06`$) should be even higher than that shown in Figure 3 and the discrepancy between the upper thin solid line and the thick dashed line would become worse. However, the new $`𝒢_{\mathrm{acc}}^{A_{},\mathrm{QSO}}(M_{\mathrm{ir}},t_0)`$ for $`a_{\mathrm{QSO}}=1`$ (lower thin solid line) is close to or only slightly higher than the thick dashed line, which is not seriously inconsistent with inequality (54), and any inconsistency can be reduced further because rapidly rotating BHs can have efficiencies substantially larger than 0.1. Thus, Figure 3 suggests that: it is impossible that all BHs in luminous QSOs ($`L_{\mathrm{bol}}10^{46}\mathrm{erg}\mathrm{s}^1`$) are Schwarzschild BHs; BHs in most luminous QSOs are rapidly rotating Kerr BHs; and luminous QSOs are accreting at an efficiency $`0.1`$ (cf. Elvis, Risaliti & Zamorani 2002).

Possible solutions to the inconsistency (mainly for the comparison of local BHs with QSOs having $`a_{\mathrm{QSO}}=0`$) shown in Figure 3 have been discussed in § 3.2.1. An additional possibility is that there exists more than one BH accreting materials in a QSO. In this case, the entropy or area increase associated with a given mass accretion rate depends on the number of BHs (see footnote 2) and equation (42) is an overestimate of the area increase.

In the adiabatic case, we do not discuss the roles of obscured accretion (or other accretion mechanisms) on the BH entropy density since the total BH entropy obtained from the X-ray background is not clear and also the conclusion may be affected by the uncertainty of the mass–dispersion relations at the high BH-mass end or the QSO luminosity function at the bright end.

### 3.3 Mass-to-energy conversion efficiency and mean lifetime of QSOs

In § 3.2, we have argued that the growth of BHs with $`M_{}>10^8\mathrm{M}_{}`$ mainly occurs during optically bright QSO phases and the mass-to-energy conversion efficiency of luminous QSOs should be larger than 0.1. In this subsection, we will ignore BH mergers \[the term $`\gamma _{\mathrm{merge}}(M_{},t)`$ in eq. 24\] and estimate in more detail the mass-dependent efficiency of QSOs that is required to match the local BH mass function. We will also give an estimate of the mean lifetime of QSOs.

#### 3.3.1 Mass-to-energy conversion efficiency

We assume that the mass-to-energy conversion efficiency is a function of BH mass $`ϵ(M_{})`$, or a function of QSO bolometric luminosity $`ϵ(L_{\mathrm{bol}})`$ since $`L_{\mathrm{bol}}`$ is assumed to be a function of only BH mass $`M_{}`$. By replacing the variable $`S`$ in equation (24) by $`M_{}^{}`$ and ignoring the term $`\gamma _{\mathrm{merge}}`$, we have

$$\frac{n_M_{}(M_{}^{},t)}{t}+\frac{}{M_{}^{}}\left[\mathrm{\Psi }(L,t)\frac{\mathrm{d}L}{\mathrm{d}M_{}^{}}\frac{(1ϵ)L_{\mathrm{bol}}(M_{}^{})}{ϵc^2}\right]=0.$$
(55)
By integrating equation (55) over $`M_{}^{}`$ from $`M_{}`$ to $`\mathrm{}`$ and over $`t`$ from 0 to $`t_0`$, we find

$$ϵ(L_{\mathrm{bol}})=\eta (L_{\mathrm{bol}})/[1+\eta (L_{\mathrm{bol}})],$$
(56)
where

$$\eta (L_{\mathrm{bol}})\left(\frac{\mathrm{d}L}{\mathrm{d}M_{}}\right)\frac{L_{\mathrm{bol}}}{c^2}\frac{_0^{t_0}\mathrm{\Psi }(L,t)dt}{_{M_{}(L_{\mathrm{bol}})}^{\mathrm{}}n(M_{}^{},t_0)dM_{}^{}}.$$
(57)
By inserting the QSO luminosity function and the BH mass function in local early-type galaxies in equation (56), and assuming that QSOs have a fixed fraction of the Eddington luminosity (we examine two cases, 0.1 and 1), we obtain the efficiency $`ϵ(L_{\mathrm{bol}})`$ in Figure 4. The solid lines in the interval $`10^{46}\mathrm{erg}\mathrm{s}^1L_{\mathrm{bol}}10^{47}\mathrm{erg}\mathrm{s}^1`$ \[corresponding to the BH mass in the range $`10^8(L_{\mathrm{Edd}}/L_{\mathrm{bol}})\mathrm{M}_{}M_{}10^9(L_{\mathrm{Edd}}/L_{\mathrm{bol}})\mathrm{M}_{}`$\] denote the region where the observational constraints are fairly strong. The dotted lines at the luminous and faint ends are obtained by extrapolating the BH distribution in early-type galaxies or the QSO luminosity function. We also show separately the results obtained by using the mass–dispersion relations in Tremaine et al. (2002) (thick dotted and solid lines) and in Merritt & Ferrarese (2001a) (thin dotted and solid lines). The lower thick/thin solid and dotted lines represent the results by assuming $`L_{\mathrm{bol}}=L_{\mathrm{Edd}}`$; and the upper ones represent the results by assuming $`L_{\mathrm{bol}}=0.1L_{\mathrm{Edd}}`$. The fact that the efficiency is required to be higher than 0.3 (the maximum efficiency for thin-disk accretion onto BHs) or even approach 1 when $`L_{\mathrm{bol}}=0.1L_{\mathrm{Edd}}`$ suggests that the assumption of sub-Eddington luminosity for luminous QSOs is unrealistic, a conclusion derived by different arguments in § 3.2.1. Figure 4 shows that even for $`L_{\mathrm{bol}}=L_{\mathrm{Edd}}`$, when $`L_{\mathrm{bol}}`$ approaches $`10^{47}\mathrm{erg}\mathrm{s}^1`$ the mass-to-energy efficiency may be slightly larger than 0.3 (the maximum efficiency for thin-disk accretion onto BHs). As described in § 3.2.1, this unrealistic $`ϵ`$ may be caused by (i) ignoring BH mergers, (ii) ignoring ejection of BHs from galactic nuclei or BHs left in the halo after a galaxy merger, (iii) uncertainties in the mass–dispersion relation at the high-mass end or the QSO luminosity function at the bright end, or (iv) ignoring the possibility of super-Eddington luminosities.

#### 3.3.2 Mean lifetime

If the duty cycle of QSOs with BH mass $`M_{}`$ is $`\delta (M_{},t)`$ at time $`t`$, we may define the mean lifetime of QSOs as follows:

$$\tau (>M_{})\frac{_M_{}^{\mathrm{}}dM_{}^{}_0^{t_0}\delta (M_{}^{},t)n_M_{}(M_{}^{},t)dt}{_M_{}^{\mathrm{}}n_M_{}(M_{}^{},t_0)dM_{}^{}}.$$
(58)
With this definition, we have

$`\tau (>M_{})={\displaystyle \frac{_{L(M_{})}^{\mathrm{}}dL_0^{t_0}\mathrm{\Psi }(L,t)dt}{_M_{}^{\mathrm{}}n_M_{}(M_{}^{},t_0)dM_{}^{}}}{\displaystyle \frac{_{L(M_{})}^{\mathrm{}}dL_0^{t_0}\mathrm{\Psi }(L,t)dt}{_M_{}^{\mathrm{}}n_M_{}^{\mathrm{early}}(M_{}^{},t_0)dM_{}^{}}};`$ (59)
the approximate equality holds for $`M_{}10^8\mathrm{M}_{}`$, where early-type galaxies dominate the number density. The lifetime $`\tau (>M_{})`$ in equation (59) is independent of any assumed value of the efficiency $`ϵ`$. By inserting the QSO luminosity function and the BH mass function in local early-type galaxies in equation (59), and assuming that QSOs have the Eddington luminosity, we obtain the lifetime $`\tau (>M_{})`$ as a function of $`M_{}`$ in Figure 5. As seen from Figure 5, the mean lifetime obtained by using the mass–dispersion relation in Tremaine et al. (2002) (thick solid line) is $`\tau (>M_{})`$ (3–13)$`\times 10^7\mathrm{yr}`$, i.e., the average duty cycle of QSOs is about $`\tau /t_{\mathrm{Hubble}}`$ 3–13$`\times 10^3`$, where $`t_{\mathrm{Hubble}}`$ is the Hubble time $`10^{10}\mathrm{yr}`$. This range is comparable to the Salpeter time $`4.5\times 10^7ϵ/[0.1(1ϵ)]\mathrm{yr}`$, the time for a BH accreting with the Eddington luminosity to e-fold in mass (for a BH with spin in an equilibrium state, the time for a BH radiating at the Eddington luminosity to e-fold in entropy is half of the Salpeter time). Hence, accretion during bright QSO phases can significantly increase the mass of BHs, which is once again consistent with the argument that growth of high-mass BHs comes mainly from accretion during optically bright QSO phases. The lifetimes of QSOs derived here are also consistent with some theoretical models for the luminosity function of QSOs (e.g. Haehnelt, Natarajan & Rees 1998; Haiman & Loeb 1998) and with the spatial clustering of QSOs (e.g. Haiman & Hui 2001; Martini & Weinberg 2001).

## 4 Conclusions

In this paper, we have studied the observational constraints on the growth of massive BHs in galactic nuclei. We use the velocity dispersions of early-type galaxies obtained by the SDSS and the relation between BH mass and velocity dispersion to estimate the local BH mass density to be $`\rho _{}(z=0)(2.5\pm 0.4)\times 10^5h_{0.65}^2\mathrm{M}_{}\mathrm{Mpc}^3`$ (eq. 13). We also use the QSO luminosity function from the 2dF QSO Redshift Survey to estimate the BH mass density accreted during optically bright QSO phases to be $`\rho _{,\mathrm{acc}}^{\mathrm{QSO}}(z=0)=2.1\times 10^5(C_B/11.8)[0.1(1ϵ)/ϵ]\mathrm{M}_{}\mathrm{Mpc}^3`$ (eq. 23). These two results are consistent if the QSO mass-to-energy conversion efficiency $`ϵ`$ is $`0.1`$.

By studying a continuity equation for the BH distribution and including the effect of BH mergers in two extreme cases (the classical case and the adiabatic case), we predict relations between the local BH mass function and the QSO luminosity function. In the classical case, the predicted relation is not consistent with the observations at high BH masses ($`M_{}10^8\mathrm{M}_{}`$) unless luminous QSOs ($`L_{\mathrm{bol}}10^{46}\mathrm{erg}\mathrm{s}^1`$) have an efficiency higher than 0.1 (e.g. $`ϵ`$0.2, which is possible for thin-disk accretion onto a Kerr BH; see Elvis, Risaliti & Zamorani 2002 for a similar conclusion). Possible ways to evade this conclusion include super-Eddington luminosities (Begelman, 2001, 2002) or ejection of BHs from galactic centers, but if these are not important then we will come to the following conclusions. Accretion other than what is traced by optically bright QSOs is not important for the growth of high-mass ($`10^8\mathrm{M}_{}`$) BHs (e.g. accretion by obscured QSOs, accretion via advection dominated accretion flow, or accretion of non-baryonic dark matter). If the growth of low-mass BHs also occurs mainly during optically bright QSO phases, less luminous QSOs ($`L_{\mathrm{bol}}10^{46}\mathrm{erg}\mathrm{s}^1`$) must accrete with a low efficiency $`0.1`$; alternatively, other mechanisms such as obscured accretion may contribute to the mass density of low-mass BHs, and their efficiency could be similar to high-mass BHs, $`ϵ`$0.2. The comparison of the observations with the predicted relation in the adiabatic case also suggests that BHs in most luminous QSOs should be rapidly rotating Kerr BHs with an efficiency $`0.1`$, or there should exist more than one BH in a QSO. We may expect that these constraints on BH demography are useful for estimating the gravitational wave signal from BH mergers that may be detected by Laser Interferometer Space Antenna (LISA).

We estimate the mean lifetime of luminous QSOs with BH mass in the range $`10^8\mathrm{M}_{}M_{}10^9\mathrm{M}_{}`$, which is (3–13)$`\times 10^7\mathrm{yr}`$ and comparable to the Salpeter time if $`ϵ`$0.1–0.3, which is once again consistent with the argument that the growth of high-mass BHs ($`M_{}10^8\mathrm{M}_{}`$) comes mainly from accretion during optically bright QSO phases.

We are grateful to Amy Barger, Mitchell Begelman, Xiaohui Fan, Jeremy Goodman, Martin Haehnelt, Zoltán Haiman, Youjun Lu, David Schlegel, David Spergel and Michael Strauss for helpful discussions. This research was supported in part by NSF grant AST-9900316 and by NASA grant HST-GO-09107.09-A.


text:
# An extrapolation method for shell model calculations

## Abstract

We propose a new shell model method, combining the Lanczos digonalization and extrapolation method. This method can give accurate shell model energy from a series of shell model calculations with various truncation spaces, in a well-controlled manner. Its feasibility is demonstrated by taking the $`fp`$ shell calculations.

PACS numbers: 21.60.-n, 21.60.Cs

Quantum many-body phenomena frequently reveal novel and emerging nature beyond a simple extrapolation from the single-particle picture, mean-field predictions, and exact diagonalization of Hamiltonian matrices by such as the Lanczos method. Solving such many-body problems has been a common challenging issue in elementary particle physics, nuclear structure physics, condensed matter physics and others. In these decades, diverse numerical methods, such as quantum Monte Carlo, stochastic diagonalization/variation, density matrix renormalization group (DMRG), were proposed. Among them, as a method combining the Path-integral formalism and diagonalization/renormalization group, in nuclear structure physics, the Quantum Monte Carlo diagonalization (QMCD) method  was proposed and was developed. In condensed matter physics, the Path-integral renormalization group (PIRG) method was proposed. Though both share common procedure, in the latter method, an extrapolation method is utilized. This extrapolation method is significant by itself and can be applied, independently of other parts of the PIRG method. In the present paper, we combine this extrapolation method with the Lanczos diagonalization and propose a new shell model method.

The shell model is one of the fundamental frameworks for nuclear many-body problem. In the conventional shell model calculation, Hamiltonian matrix elements of the effective shell model interaction are calculated by a complete basis set in a given shell model space, and then, diagonalization is carried out by the Lanczos method. Such Lanczos shell model calculations have a long history over a half century and have extended its realm of the application due to the progress of computer technologies. Nowadays, part of full $`fp`$ shell problem can be solved, for instance, the yrast states of 52Fe, the $`M`$ scheme dimension of which in the complete $`fp`$ shell model space is 109,954,620, was studied, while there still exist many nuclei not yet explored in the full $`fp`$ shell due to their huge shell model dimensions.

For such nuclei, truncation scheme is utilized for reducing the shell model space. In the $`fp`$ shell calculations, due to relatively large gap of spherical single particle energy between the $`f_{7/2}`$ orbit and others, we restrict the shell model space by specifying the number of excited nucleons across this shell gap, that is, $`_{st}(f_{7/2})^{A40s}(r)^s`$ where $`r`$ means the set of the $`f_{5/2}`$, $`p_{3/2}`$ and $`p_{1/2}`$ orbits and $`t`$ is the maximum number of nucleons allowed to be excited. As an example, we consider 48Cr within $`fp`$ shell model space, of which shell model calculation needs about 2 million dimension. This was one of the state-of-the-art large-scale shell model calculations in mid 1990’s . In Fig. 1 (a) and 1(b), the ground state energies of the truncated space are plotted as functions of the $`t`$ value and the corresponding $`M`$ scheme dimensions, respectively. The definite behavior of the convergence is not seen in Fig. 1, because energy eigenvalue as a function of the $`t`$ or $`M`$ scheme dimension does not have a well-defined scaling property. Therefore, in general, we cannot predict the ground state energy in the entire shell model space by extrapolating shell model results within smaller truncated spaces unless the truncation calculations give well-converged results. Hence, it is crucial to find a useful extrapolation scheme to extract the correct results from truncation calculations.

In order to extrapolate the shell model results within smaller truncated spaces into the complete shell model result, a well-defined scaling property for energy eigenvalues is needed. We define the difference $`\delta E`$ between the energy eigenvalue $`\widehat{H}`$ in a given truncated space and true energy eigenvalue $`\widehat{H}_g`$, that is, $`\delta E=\widehat{H}\widehat{H}_g`$. The energy variance $`\mathrm{\Delta }E`$ in a given truncated spaces also defined as, $`\mathrm{\Delta }E=\frac{\widehat{H}^2\widehat{H}^2}{\widehat{H}^2}.`$ The difference $`\delta E`$ vanishes linearly as a function of the energy variance $`\mathrm{\Delta }E`$. This relation was utilized in Refs.  and  and the summary of its proof can be seen in Ref. . An approximate ground state $`|\psi `$ can be decomposed by the true eigenstate $`|\psi _g`$ and rest component $`|\psi _e`$ as $`|\psi =c|\psi _g+d|\psi _e`$ where $`c^2+d^2=1`$, and $`|\psi _g`$ and $`|\psi _e`$ are orthnormalized states. Up to $`O(d^3)`$, $`\delta E\mathrm{\Delta }E`$ is satisfied. Therefore, this relation shows a good scaling if $`|\psi `$ is a good approximation of the true eigenstate.

As an example, we reconsider the shell model results of 48Cr by the energy variance $`\mathrm{\Delta }E`$ as the $`x`$-axis. In Fig. 2, energy eigenvalues for various truncation spaces are plotted as a function of the energy variance for the $`0_1^+`$ and $`2_1^+`$ states of 48Cr with the KB3 interaction. Linear relation between $`\delta E`$ and $`\mathrm{\Delta }E`$ is seen from $`t=1`$. By the results corresponding to the $`t=1,2,3`$, $`t=3,4,5`$, $`t=4,5,6`$, we can extrapolate the ground state energy as -33.148 $`\pm `$ 0.15, -33.047 $`\pm `$ 0.02, -32.975 $`\pm `$ 0.02 MeV, respectively, where these error-bars are given by the $`\chi ^2`$ fitting. As the exact ground state energy is -32.954 MeV, the extrapolation from $`t=1,2,3`$ results is already a good value. Moreover, as the truncation space becomes larger, the error due to the extrapolation procedure becomes smaller. Thus the extrapolated energy is given with error-bar, and its error can be well-controlled. In the SMMC calculation, a different extrapolation method is utilized for overcoming the minus sign problem . The extrapolated energy of this nucleus is -32.3 $`\pm `$ 0.4 MeV. This error comes from two reasons. One is a statistical uncertainties due to the Monte Carlo sampling, the other is an systematic error due to the extrapolation and finite temperature corrections. Compared with the SMMC calculation, the present method has only a small systematic error. Here, we comment that, extrapolated energies may become lower than the exact energy because this process does not necessarily guarantee the variational principle.

Next we consider the expectation value of various operators by the present extrapolation method. In Fig. 3, we show the quadrupole moment, occupation number and B(E2,$`2_1^+0_1^+`$) for 48Cr. As the first two values are expectation values by the same wavefunction, the $`\mathrm{\Delta }E`$ can be uniquely defined. We can extrapolate them into zero energy variance for quadrupole moments. In Fig. 3(a), unlike energy, $`t=1,2`$ results seem to be rather poor. From the $`t=1,2,3`$, $`t=3,4,5`$ results, the extrapolated quadrupole moments are -24.6 $`\pm `$ 4.1, -29.0 $`\pm `$ 0.3 $`efm^2`$, respectively. The latter extrapolated value reaches almost exact one (-29.1 $`efm^2`$). On the other hand, for the B(E2), the initial and final states are different. Therefore, initial and final states have different $`\mathrm{\Delta }E`$’s, which makes extrapolation procedure difficult. In this case, however, the shell model spaces with the same $`t`$ value but different total magnetic quantum number have similar $`\mathrm{\Delta }E`$. Therefore, average $`\mathrm{\Delta }E`$ works quite well in Fig. 3(c). The extrapolated B(E2) from $`t=2,3,4`$ results is 220 $`\pm `$ 7 $`e^2fm^4`$, which very well reproduces the exact one ( 220 $`e^2fm^4`$ ).

By this example, we conclude that energy vs. energy variance plot is useful to get an extrapolated energy from a series of shell model results with various truncated spaces. In addition to the energy estimates, the same extrapolation is also shown to be possible for quadrupole moments, B(E2) and other quantities.

Moreover we comment on the present extrapolation from the structure of 48Cr. Low-lying yrast states of 48Cr is deformed, and for its description, the particle-hole excitation across the shell gap between the $`f_{7/2}`$ orbit and others is essential. Therefore it is interesting that the extrapolation from the truncated shell model results with few particle-hole excitations works quite well.

We proceed further large-scale shell model problem. The $`M`$ scheme dimension of 54Fe is about 0.35 billion and several recent methods solved this shell model problem. For astrophysical interest, the Gamow Teller transition of this nucleus is known to be important. The Gamow Teller transition also offers good benchmark test for the present extrapolation. For the following calculation of 54Fe, the FPD6 interaction is taken as the effective interaction. Fig. 4(a) shows the ground state energies of $`t=`$1,2,3,4 and 5 shell model spaces. The extrapolated ground state energy from $`t=3,4,5`$ is -174.49 $`\pm `$ 0.02 MeV. Here we compare this value to the energies of other methods. The ground state energy of the $`t=7`$ shell model spaces, of which $`M`$ scheme dimension is 91,848,462, is -174.217 MeV. The extrapolated energy of the SMMC calculation  is -172.9 $`\pm `$ 0.6 MeV. In Ref., from the statistical point of view, the $`JT`$ scheme bases are selected and then the Hamiltonian matrix is diagonalized. This energy is about -172.9 MeV . The ground state energy of the QMCD method with ten $`J`$-projected bases  is -173.9 MeV. The present extrapolated energy is consistent with the QMCD and $`t=7`$ truncated shell model energies.

In Fig. 4(b), the total GT+ strength is plotted. The extrapolated total GT+ from $`t=3,4,5`$ result is 4.98 $`\pm `$ 0.06. The total GT+ strength of $`t=7`$ calculation is 5.03. In Fig.2 of Ref. , the total GT+ strength is plotted as a function of the $`t`$ value and it suggests that its total strength is about 5. On the other hand, in Ref. , the extrapolated value of the SMMC calculation is 4.32 $`\pm `$ 0.24. Our extrapolated value is consistent with the former two values.

The 56Ni is a key nucleus for understanding the shell structure around the $`fp`$ shell region because it has a doubly magic structure, while its closed shell is suggested to be rather soft. Moreover, recently two deformed bands are discovered . This nucleus is quite interesting but the $`M`$ scheme dimension is quite large ( about $`1.1\times 10^9`$ ). Then this nucleus was a target of a state-of-the-art shell model calculation in the end of 1990’s . Therefore it is also a good touchstone for the present extrapolation method.

In Fig. 5, the energies of the $`0_1^+`$, $`2_1^+`$, $`2_2^+`$ and $`2_3^+`$ states with $`t=26`$ shell model spaces are plotted. In principle, the present extrapolation method can be applied to the excited states with the same quantum numbers without any difficulty. Here we test its numerical feasibility. The extrapolated ground state energy from $`t=46`$ is -203.26 $`\pm `$ 0.06 MeV. In Ref. , the ground state energy of the QMCD calculation is -203.152 MeV. The ground state energy of the $`t=8`$ truncated space is -203.178 MeV, of which $`M`$ scheme dimension is 255,478,309. The extrapolated ground state energy is consistent with the results of the QMCD calculation and the $`t=8`$ truncated shell model calculation. The extrapolated $`2_1^+`$ energy is -200.46 $`\pm `$ 0.05 MeV from $`t=46`$ results, which is also consistent with the QMCD result and the $`t=8`$ truncated shell model result (-200.075 MeV).

Next we investigate the excited $`2^+`$ energies. In the left part of the Fig. 5, the QMCD energies are shown. These $`2^+`$ states are characterized by the quadrupole moments and the occupation numbers in the $`f_{7/2}`$ orbit. The quadrupole moments are 21.6, -41.9 and 9.5 $`efm^2`$ for the $`2_1^+`$, $`2_2^+`$ and $`2_3^+`$ states, respectively. The occupation numbers are 12.9, 11.6 and 12.1. In Ref. , the $`2_2^+`$ state is concluded to be a prolate deformed state. As the truncation scheme is based on the spherical single particle energies, the deformed $`2_2^+`$ state is considered to be difficult to describe within small $`t`$ truncated spaces. For $`t6`$, the states with negative quadrupole moment appear as the third $`2^+`$ state. By plotting the shell model energies for the $`t=26`$ truncated spaces, we can find the linear relations for $`2_1^+`$, $`2_2^+`$ and $`2_3^+`$ energies, respectively, as shown in Fig. 5. Therefore, for excited states, the present extrapolation method can be applied. Moreover Fig. 5 shows crossing between the $`2_2^+`$ and $`2_3^+`$ energies as a function of the energy variance, which means that, the present extrapolation method can correctly handle these excited states with different nature. The extrapolated energy of the $`2_2^+`$ state from $`t=36`$ results is -198.35 $`\pm `$ 0.35 MeV and extrapolated quadrupole moment and occupation number are -28.2 $`\pm `$ 8.5 $`efm^2`$ and 11.5 $`\pm `$ 0.6, respectively, which are consistent with the QMCD values. For more presice description of this prolate deformed state, $`t=7`$ results are needed, while it is quite surprising that such a deformed excited state can be described by the present extrapolation method based on the small truncated shell model spaces.

Finally we discuss a computational aspect. For the present extrapolation method, in addition to the usual Lanczos procedure, the evaluation of the matrix element $`\widehat{H}^2`$ is required. To obtain this matrix element, we operate the Hamiltonian $`H`$ to the $`|\psi `$ of the $`t`$ truncated space in the ($`t+2`$) truncated space, then we can evaluate the norm of the $`H|\psi `$. Therefore, the implementation of this extrapolation method is quite easy for a usual shell model program.

In summary, we have proposed a new shell model method based on the Lanczos diagonalization and an extrapolation method. By taking the state-of-the-art $`fp`$ shell calculations in the last decade, we have shown its numerical feasibility. The present extrapolation method can give accurate energies, quadrupole moments and other quantities from a series of shell model results with smaller truncated shell model spaces.

The present extrapolation method was successfully applied to condensed matter problems as a part of the PIRG method. Therefore, it is intriguing to examine whether this extrapolation method is useful also in the QMCD method for the nuclear shell model. For instance, we expect more precise energy estimates from the QMCD calculations with smaller number of bases. Moreover, it might become quite useful in other methods of nuclear structure physics. Such researches are in progress.

We thank Professor N. Yoshinaga for reading the manuscript and useful comments. This work was supported from ’Research for the Future Progam’ by Japan Society for the Promotion of Science under grant RFTF97P01103, and was supported in part by Grant-in-Aid for Scientific Research (A)(2) (10304019), and one for Specially Promoted Research (13002001) from the Ministry of Education, Science and Culture.