Supermassive Black Holes in the Hierarchical Universe: A General Framework and Observational Tests
Abstract
We present a simple framework for the growth and evolution of supermassive black holes (SMBHs) in the hierarchical structure formation paradigm, adopting the general idea that quasar activity is triggered in major mergers. In our model, black hole accretion is triggered during major mergers (mass ratio ) between host dark matter halos. The successive evolution of quasar luminosities follows a universal light curve form, during which the growth of the SMBH is modeled selfconsistently: an initial exponential growth at constant Eddington ratio of order unity until it reaches the peak luminosity, followed by a powerlaw decay. Assuming that the peak luminosity correlates with the postmerger halo mass, we convolve the light curve with the triggering rate of quasar activity to predict the quasar luminosity function (LF). Our model reproduces the observed LF at for the full luminosity ranges probed by current optical and Xray surveys. At , our model underestimates the LF at , allowing room for AGN activity triggered by secular processes instead of major mergers. At , in order to reproduce the observed quasar abundance, the typical quasar hosts must shift to lower mass halos, and/or minor mergers can also trigger quasar activity. Our model reproduces both the observed redshift evolution and luminosity dependence of the linear bias of quasar/AGN clustering. Due to the scatter between instantaneous luminosity and halo mass, quasar/AGN clustering weakly depends on luminosity at low to intermediate luminosities; but the linear bias rises rapidly with luminosity at the high luminosity end and at high redshift. In our model, the Eddington ratio distribution is roughly lognormal, which broadens and shifts to lower mean values from high luminosity quasars () to low luminosity AGNs (), in good agreement with observations. The model predicts that the vast majority of SMBHs were already in place by , and of them were in place by ; but the less massive () SMBHs were assembled more recently, likely more through secular processes than by major mergers – in accordance with the downsizing picture of SMBH assembly since the peak of bright quasar activity around .
Subject headings:
black hole physics – galaxies: active – cosmology: observations – largescale structure of universe – quasars: general – surveys1. Introduction
The cosmic assembly and evolution of supermassive black holes
(SMBHs) is a central topic in modern cosmology: in the local
universe, SMBHs reside in almost every bulge dominant galaxy
(e.g., Kormendy & Richstone, 1995; Richstone et al., 1998); and
they likely played important roles during their coevolution with
galaxy bulges
(e.g., Silk & Rees, 1998; Wyithe & Loeb, 2002, 2003; King, 2003; Di Matteo et al., 2005; Croton et al., 2006), as inferred from
observed scaling relations between bulge properties and the mass
of the central BH
(e.g., Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Graham et al., 2001; Tremaine et al., 2002; Marconi & Hunt, 2003). It is also generally
accepted that the local dormant SMBH population was largely
assembled via gas accretion during a luminous QSO
phase
One of the leading hypotheses for the QSO triggering mechanism is galaxy mergers (e.g., Hernquist, 1989; Carlberg, 1990; Kauffmann & Haehnelt, 2000; Hopkins et al., 2008, and references therein). In the hierarchical CDM paradigm, small structures merge to form large structures, and the merger rate of dark matter halos peaks at early times, in broad consistency with the observed peak of bright quasar activities. Gasrich major mergers (i.e., those with comparable mass ratios) between galaxies provide an efficient way to channel a large amount of gas into the central region to trigger starbursts and possibly feed rapid black hole growth. Observationally this merger hypothesis is supported by the ULIRGquasar connection (e.g., Sanders & Mirabel, 1996; Canalizo & Stockton, 2001), the signature of recent mergers in quasar hosts (e.g., Bennert et al., 2008), and the smallscale overdensities of galaxies around luminous quasars or quasar binaries (e.g., Fisher et al., 1996; Bahcall et al., 1997; Serber et al., 2006; Hennawi et al., 2006; Myers et al., 2008; Strand et al., 2008). These observations do not necessarily prove that mergers are directly responsible for triggering QSO activity, but they at least suggest that QSO activity is coincident with mergers in many cases.
The last decade has seen a number of models for QSO evolution based on the merger hypothesis (e.g., Haiman & Loeb, 1998; Kauffmann & Haehnelt, 2000; Wyithe & Loeb, 2002, 2003; Volonteri et al., 2003; Hopkins et al., 2008), which agree with the bulk of QSO observations reasonably well. Motivated by the success of these mergerbased QSO models, we revisit this problem in this paper with improved observational data on SMBHs and QSOs from dedicated large surveys, and with updated knowledge of the merger rate as inferred from recent numerical simulations. The present work is different from previous studies in many aspects in both methodology and the observational tests used. Our goal is to check if a simple mergerbased cosmological QSO framework can reproduce all the observed statistics of SMBHs (both active and dormant), and to put constraints on the physical properties of SMBH growth. Our model is observationally motivated, therefore we do not restrict ourselves to theoretical predictions of SMBH/QSO properties from either cosmological or merger event simulations.
The paper is organized as follows. In §1.1 we review some aspects of local SMBH demographics and QSO luminosity function; in §1.2 we briefly review the current status of quasar clustering observations; §§1.3 and 1.4 discuss the halo and subhalo merger rate from numerical simulations. Our model formulism is presented in §2 and we present our fiducial model and compare with observations in §3. We discuss additional aspects of our model in §4 and conclude in §5.
We use friendsoffriends (fof) mass with a link length to define the mass of a halo (roughly corresponding to spherical overdensity mass with times the mean matter density), since the fitting formulae for the halo mass function are nearly universal with this mass definition (Sheth et al., 2001; Jenkins et al., 2001; Warren et al., 2006; Tinker et al., 2008). For simplicity we neglect the slight difference between fof mass and virial mass (Appendix A; and see White, 2002, for discussions on different mass definitions). Throughout the paper, always refers to the bolometric luminosity, and we use subscripts or to denote band or Xray luminosity when needed. We adopt a flat CDM cosmology with , , , , , . We use the Eisenstein & Hu (1999) transfer function to compute the linear power spectrum, and use the fitting formulae for halo abundance from Sheth & Tormen (1999) and for halo linear bias from Sheth et al. (2001) based on the ellipsoidal collapse model and tested against numerical simulations.
1.1. Supermassive Black Hole Demographics
In the local universe, the dormant BH mass function (BHMF) can be estimated by combining scaling relations between BH mass and galaxy bulge properties, such as the relation (e.g., Gebhardt et al., 2000; Ferrarese & Merritt, 2000; Tremaine et al., 2002) or the relation (e.g., Magorrian et al., 1998; Marconi & Hunt, 2003), with bulge luminosity or stellar mass/velocity dispersion functions (e.g., Salucci et al., 1999; Yu & Tremaine, 2002; Marconi et al., 2004; Shankar et al., 2004; McLure & Dunlop, 2004; Yu & Lu, 2004, 2008; Tundo et al., 2007; Shankar et al., 2009b). Although some uncertainties exist on the usage of these scaling relations at both the high and low BH mass end (e.g., Lauer et al., 2007; Tundo et al., 2007; Hu, 2008; Greene et al., 2008; Graham, 2008; Graham & Li, 2009), the total BH mass density is estimated to be with an uncertainty of a factor (e.g. Yu & Lu, 2008; Shankar et al., 2009b) – but the shape of the local BHMF is more uncertain (see the discussion in Shankar et al., 2009b).
It is now generally accepted that local SMBHs were once luminous QSOs (Salpeter, 1964; Zel’dovich & Novikov, 1964; LyndenBell, 1969). An elegant argument tying the active QSO population to the local dormant SMBH population was proposed by Soltan (1982): if SMBHs grow mainly through a luminous QSO phase, then the accreted luminosity density of QSOs to should equal the local BH mass density:
(1) 
where is the bolometric luminosity function (LF) per interval. Given the observed QSO luminosity function, a reasonably good match between and can be achieved if the average radiative efficiency (e.g., Yu & Tremaine, 2002; Shankar et al., 2004; Marconi et al., 2004; Hopkins et al., 2007; Shankar et al., 2009b). The exact value of , however, is subject to some uncertainties from the luminosity function and local BH mass density determination.
Assuming that SMBH growth is through gas accretion, an extended
version of the Soltan argument can be derived using the continuity
equation
(2) 
where is the BH mass function per , and is the mean accretion rate of BHs with mass averaged over both active and inactive BHs. Given the luminosity function, and a model connecting luminosity to BH mass, one can derive the BH mass function at all times by integrating the continuity equation (e.g., Marconi et al., 2004; Merloni, 2004; Shankar et al., 2009b).
Since the local BHMF and the QSO luminosity function together determine the assembly history of the cosmic SMBH population, any cosmological model of AGN/quasars must first reproduce the observed luminosity function (e.g., Haiman & Loeb, 1998; Kauffmann & Haehnelt, 2000; Wyithe & Loeb, 2002, 2003; Volonteri et al., 2003; Lapi et al., 2006; Hopkins et al., 2008; Shankar et al., 2009b). This is also one of the central themes of this paper.
There has been great progress in the measurements of the QSO luminosity function over wide redshift and luminosity ranges for the past decade, mostly in the optical band (e.g., Fan et al., 2001, 2004; Wolf et al., 2003; Croom et al., 2004; Richards et al., 2005, 2006; Jiang et al., 2006; Fontanot et al., 2007), and the soft/hard Xray band (e.g., Ueda et al., 2003; Hasinger et al., 2005; Silverman et al., 2005; Barger et al., 2005). While widefield optical surveys provide the best constraints on the bright end of the QSO luminosity functions, deep Xray surveys can probe the obscured faint end AGN population. Although it has been known for a while that the spatial density of opticallyselected bright quasars peaks at redshift , the spatial density of fainter AGNs selected in the Xray seems to peak at lower redshift (e.g., Steffen et al., 2003; Ueda et al., 2003; Hasinger et al., 2005), a trend now confirmed in other observational bands as well (e.g., Bongiorno et al., 2007; Hopkins et al., 2007, and references therein) – the manifestation of the socalled downsizing of the cosmic SMBH assembly.
Hopkins et al. (2007) compiled QSO luminosity function data from surveys in various bands (optical, Xray, IR, etc). Assuming a general AGN spectral energy distribution (SED) and a column density distribution for obscuration, Hopkins et al. (2007) were able to fit a universal bolometric luminosity function based on these data. We adopt their compiled bolometric LF data in our modeling. The advantage of using the bolometric LF data is that both unobscured and obscured SMBH growth are counted in the model; but we still suffer from the systematics of bolometric corrections. The nominal statistical/systematic uncertainty level of the bolometric LF is (cf., Shankar et al., 2009b).
1.2. Quasar Clustering
A new ingredient in SMBH studies, due to dedicated largescale optical surveys such as the 2QZ (Croom et al., 2004) and SDSS (York et al., 2000), is quasar clustering. While quasar clustering studies can be traced back to more than two decades ago (e.g., Shaver, 1984), statistically significant results only came very recently (e.g., Porciani et al., 2004; Croom et al., 2005; Porciani & Norberg, 2006; Myers et al., 2006, 2007a, 2007b; Shen et al., 2007; da Ângela et al., 2008; Padmanabhan et al., 2008; Shen et al., 2008b, 2009; Ross et al., 2009).
Within the biased halo clustering picture (e.g., Kaiser, 1984; Bardeen et al., 1986; Mo & White, 1996), the observed QSO twopoint correlation function implies that quasars live in massive dark matter halos and are biased tracers of the underlying dark matter. For optically luminous quasars (), the host halo mass inferred from clustering analysis is a few times . Thus quasar clustering provides independent constraints on how SMBHs occupy dark matter halos, where the abundance and evolution of the latter population can be well studied in analytical theories and numerical simulations (e.g., Press & Schechter, 1974; Bond et al., 1991; Lacey & Cole, 1993; Mo & White, 1996; Sheth & Tormen, 1999; Sheth et al., 2001; Springel et al., 2005b). By comparing the relative abundance of quasars and their host halos, one can constrain the average quasar duty cycles or lifetimes (e.g., Cole & Kaiser, 1989; Martini & Weinberg, 2001; Haiman & Hui, 2001) to be yr, which means bright quasars are short lived.
More useful constraints on quasar models come from studies of quasar clustering as a function of redshift and luminosity. The redshift evolution of quasar clustering constrains the evolution of duty cycles of active accretion, while the luminosity dependence of quasar clustering puts constraints on quasar light curves (LC). Successful cosmological quasar models therefore not only need to account for quasar abundances (LF), but also need to explain quasar clustering properties, both as function of redshift and luminosity. We will use quasar clustering observations as additional tests on our quasar models, as has been done in recent studies (e.g., Hopkins et al., 2008; Shankar et al., 2009a; Thacker et al., 2009; Bonoli et al., 2009; Croton, 2009).
The luminosity dependence of quasar bias can be modeled as follows. Assume at redshift the probability distribution of host halo mass given bolometric luminosity is , then the halo averaged bias factor is:
(3) 
where is the linear bias factor for halos with mass at redshift . The derivation of the probability distribution is described in later sections [Eqn. (2.1)].
1.3. Dark Matter Halo Mergers
In the hierarchical universe small density perturbations grow under gravitational forces and collapse to form virialized halos (mostly consisting of dark matter). Smaller halos coalesce and merge into larger ones. Early work on the merger history of dark matter halos followed the extended PressSchechter (EPS) theory (e.g., Press & Schechter, 1974; Bond et al., 1991; Lacey & Cole, 1993), which is based on the spherical collapse model (Gunn & Gott, 1972) with a constant barrier for the collapse threshold (the critical linear overdensity is independent on mass, albeit it depends slightly on cosmology in the CDM universe). Although the (unconditional) halo mass function predicted by the EPS theory agrees with numerical body simulations reasonably well, it is well known that it overpredicts (underpredicts) the halo abundance at the low (high) mass end (Sheth & Tormen, 1999, and references therein). Motivated by the fact that halo collapses are generally triaxial, Sheth & Tormen (1999), based on earlier work by Bond & Myers (1996), proposed the ellipsoidal collapse model with a moving barrier where the collapse threshold also depends on mass. By imposing this mass dependence on collapse barrier, the abundance of small halos is suppressed and fitting formulae for the (unconditional) halo mass function are obtained, which agree with body simulations much better than the spherical EPS predictions (e.g., Sheth & Tormen, 1999; Sheth et al., 2001; Jenkins et al., 2001; Warren et al., 2006; Tinker et al., 2008).
In addition to the unconditional halo mass function, the conditional mass function gives the mass spectrum of progenitor halos at an earlier redshift of a descendant halo at redshift . This conditional mass function thus contains information about the merger history of individual halos and can be used to generate halo merger trees. A simple analytical form for can be obtained in the spherical EPS framework (e.g., Lacey & Cole, 1993); for the ellipsoidal collapse model, an exact but computationally consuming solution of the conditional mass function can be obtained by solving the integral equation proposed by Zhang & Hui (2006). Recently, Zhang et al. (2008) derived analytical formulae for in the limit of small lookback times for the ellipsoidal collapse model.
Alternatively, the halo merger rate can be retrieved directly from large volume, high resolution cosmological body simulations, which bypasses the inconsistencies between the spherical EPS theory and numerical simulations. Using the product of the Millennium Simulation (Springel et al., 2005b), Fakhouri & Ma (2008) quantified the mean halo merger rate per halo for a wide range of descendant (at ) halo mass , progenitor mass ratio and redshift , and they provided an almost universal fitting function for the mean halo merger rate to an accuracy within the numerical simulation results:
(4) 
Here is the instantaneous merger rate at redshift , for halos with mass and with progenitor mass ratio in the range [ is in units of number of mergers ], is the halo mass function (in units of ), where are the masses of the two progenitors, and is the linear density threshold for spherical collapse with the linear growth factor. We adopt Eqn. (1.3) to estimate the halo merger rate in our modeling. The mass definition used here is friendsoffriends mass with a link length . For the range of mass ratios we are interested in (i.e., major mergers), the halo merger rate derived from the spherical EPS model can overpredict the merger rate by up to a factor of (Fakhouri & Ma, 2008).
1.4. Galaxy (subhalo) Mergers
Galaxies reside in the central region of dark matter halos where the potential well is deep and gas can cool to form stars. When a secondary halo merges with a host halo, it (along with its central galaxy) becomes a satellite within the virial radius of the host halo. It will take a dynamical friction time for the subhalo to sink to the center of the primary halo, where the two galaxies merge. The subsequent galaxy (subhalo) merger rate is therefore different from the halo merger rate discussed in §1.3, and a full treatment with all the dynamics (tidal stripping and gravitational shocking of subhalos) and baryonic physics is rather complicated.
The simplest argument for the galaxy (subhalo) merger timescale is given by dynamical friction (Chandrasekhar, 1943; Binney & Tremaine, 1987):
(5) 
where and are the masses for the host
and satellite halos respectively, is the Coulomb logarithm,
is a function of the orbital energy and angular momentum of the
satellite, is an adjustable parameter and is the halo dynamical timescale, usually
estimated at the halo virial radius . This formula
(5) is valid at the small satellite mass limit ; although it is also used in cases of
in many semianalytical models
(SAMs) with the replacement of (i.e., the original definition of the
Coulomb logarithm) or . However, in recent years deviations from the predictions
by Eqn. (5) in numerical simulations have been
reported for both the and the
regimes
(e.g., Taffoni et al., 2003; Monaco et al., 2007; BoylanKolchin et al., 2008; Jiang et al., 2008).
In particular even in the regime
where the original Chandrasekhar formula is supposed to work,
Eqn. (5) substantially underestimates the merger
timescale from simulations by a factor of a few (getting worse at
lower mass ratios). This is because Eqn. (5) is
derived by treating the satellite as a rigid body, while in
reality the satellite halo loses mass via tidal
stripping
Given the limitations of analytical treatments, we therefore retreat to numerical simulation results. Fitting formulae for the galaxy merger timescales within merged halos have been provided by several groups (Taffoni et al., 2003; Monaco et al., 2007; BoylanKolchin et al., 2008; Jiang et al., 2008), and the merger rate of subhalos has also been directly measured from simulations (Wetzel et al., 2009; Stewart et al., 2008).
The fitting formula for in Jiang et al. (2008), who used hydro/body simulations, has the following form:
(6) 
where is the circularity parameter [, and is the orbit eccentricity]. In deriving this equation Jiang et al. have removed the energy dependence of the orbit, i.e., they fix where is the circular orbit that has the same energy as the satellite’s orbit. In their simulation the ratio of ranges from and has a median value . The distribution of circularity in their simulation (see their eqn. 7) has a median value , so in what follows we take . The redshift dependence of is thus the same as that of the halo dynamical time (see Appendix A).
Alternatively, the fitting formula of in BoylanKolchin et al. (2008) is given by:
(7) 
Combining the halo merger rate (1.3) and galaxy (subhalo) merger timescale (6) or (7) we derive the instantaneous merger rate of galaxies within a merged halo of mass for progenitor mass ratio and at redshift :
(8) 
where is a function of and is determined by:
(9) 
where is the cosmic time at redshift . Again, here is the number of mergers per volume per halo mass per mass ratio per redshift. We find that is almost constant at all redshifts for and monotonically decreases as increases. This constant can be approximated by its asymptotic value at large when :
(10) 
for the fitting formula of in Jiang et al. (2008), or
(11) 
for the fitting formula of in BoylanKolchin et al. (2008).
Unfortunately, current studies on the galaxy merger timescale have not converged yet, and different groups do report similar but quantitatively different results, at least partly caused by the different definitions and setups in their numerical simulations (e.g., hydro/body versus pure body simulations, halo finding algorithms and definition of mergers, etc). These fitting formulae are generally good within a factor of uncertainty. This uncertainty in the galaxy merger timescale within DM halos will lead to quite substantial differences in the galaxy merger rate at high redshift. To see this, we show in Fig. 1 the comparison of different fitting formulae for and their consequences. In the upper panel of Fig. 1 we show the ratio of for the most likely orbit with circularity and using the fitting formula in BoylanKolchin et al. (2008) (solid line), Jiang et al. (2008) (dashed line), and a commonlyused SAM prescription: (dotted line). For the major merger regime the fitting formula in BoylanKolchin et al. (2008) agrees with the SAM prescription well, and they both approach the dynamical time at the high mass ratio end. However, the fitting formula in Jiang et al. (2008) yields a factor of longer than dynamical time at the high mass ratio end.
In the lower panel of Fig. 1 we show the halo and galaxy major merger rates per unit time (integrated over ), as function of redshift. The black solid lines show the halo major merger rate, the black dashed lines show the galaxy merger rate using the fitting formula of from BoylanKolchin et al. (2008), and the red dashed lines show the galaxy merger rate adopting the fitting formula from Jiang et al. (2008). In general the galaxy merger rates fall below the halo merger rate at high redshift and take over at lower redshift. At redshift , the galaxy merger rate using the Jiang et al. (2008) formula is lower by almost an order of magnitude than that using the BoylanKolchin et al. (2008) formula (as well as the SAM prescription, since it agrees with eqn. 7 for ), which makes it difficult to produce the quasar abundance at high redshift (see later sections).
It is worth noting that although the merger rate peaks at some redshift, this trend is hierarchical such that it peaks at later time for more massive halos. This is somewhat in contradiction to the downsizing of the QSO luminosity function. This apparent discrepancy is likely caused by the fact that at low reshift, it becomes progressively more difficult for the black hole to accrete efficiently in massive halos because of the global deficit of cold gas due to previous mergers experienced by massive halos and/or possible feedback quenches (e.g., Kauffmann & Haehnelt, 2000; Croton et al., 2006). We will treat this fraction of QSOtriggering merger event as function of redshift and halo mass explicitly in our modeling.
If we choose to normalize the total merger rate () by , the abundance of descendant halos with mass , at the same redshift for halo mergers and galaxy mergers, then falls below at high redshift, then catches up and exceeds a little, and evolves more or less parallel to afterwards. Therefore the evolution in is shallower than that in , consistent with the findings by Wetzel et al. (2009) once the different definitions of and halo mass are taken into account.
To summarize §1.3 and §1.4, we have compared the halo merger rate and galaxy merger rate as functions of redshift, (post)merger halo mass and mass ratio with different prescriptions for the dynamical friction timescale. It is, however, unclear on which rate we should link to the QSO triggering rate. It is true that it will take a dynamical friction time for the two galaxies to merge after their host halo merged. But the black hole accretion might be triggered very early during their first orbital crossing, well before the galaxies merge. In what follows, we adopt the halo merger rate to model the QSO triggering rate, and we discuss the consequences of adopting the alternative subhalo merger rates in §4.2.
Symbol  Description  Value/Units 

Halo mass ratio  
Postmerger halo mass  
Black hole mass  
Halo(subhalo) merger rate  
QSO luminosity function  
QSO triggering rate  
QSO peak luminosity  
QSO triggering redshift  
QSO triggering fraction  
The relation  
Slope  
Normalization at  
Scatter  dex  
Evolution in normalization  
The light curve model  
Radiative efficiency  
Eddington luminosity per  
folding time  
Fraction of seed BH mass to peak BH mass  
Time to reach the peak luminosity  
Eddington ratio before  
Powerlaw slope of the decaying phase  
The QSOtriggering rate  
Minimum mass ratio  0.25  
Exponential lower mass cut  
Exponential upper mass cut 
Note. – Parameter values are for the fiducial model.
1.5. BH Fueling During Mergers
Modeling black hole fueling during mergers from first principles is not trivial: aside from the lack of physical inputs, most hydrodynamical simulations do not yet have the dynamical range to even resolve the outer Bondi radius. Although semianalytical treatments of BH accretion are possible (e.g., Granato et al., 2004; Monaco et al., 2007), we do not consider such treatments in this paper because of the poorly understood accretion physics.
Throughout this paper we adopt the following ansatz:

QSO activity is triggered by a gasrich major merger event ().

The buildup of the cosmic SMBH population is mainly through gas accretion during the QSO phase.

The QSO light curve follows a universal form , where the peak luminosity is correlated with the mass of the merged halo .
Once we have specified the QSO light curve model, and estimated the QSO triggering rate from the major merger rate, we can convolve them to derive the QSO luminosity function. Within this evolutionary QSO framework, we can derive the distributions of host halo masses and BH masses for any given instantaneous luminosity and redshift, allowing us to compare with observations of quasar clustering and Eddington ratio distributions. The details of this framework are elaborated in §2.
We note that, any correlations established with our model are only for hosts where a gasrich major merger is triggered and selfregulated BH growth occurs. Host halos which experience many minor mergers or dry mergers may deviate from these relations. We will come back to this point in the discussion section.
2. Model Formalism
In this section we describe our model setup in detail. Some notations and model parameters are summarized in Table 1 for clarity.
2.1. Determining the Luminosity Function
Major mergers are generally defined as events with mass ratio
, but our framework can be easily generalized to
other values of . Because it takes more than just a major
merger event to trigger QSO activity, we shall model the QSO
triggering rate as a redshift and massdependent fraction of halo
merger rate. The QSO luminosity function at a given redshift can be
obtained by combining the triggering rate and the evolution of QSO
luminosity
(12) 
where is the QSO number density in the luminosity range at redshift ; is the QSOtriggering rate (merger number density per redshift per peak luminosity) at redshift with peak bolometric luminosity , which is determined by the specific form of the light curve; and are the cosmic time at and . In integrating this equation we impose an upper limit of redshift , but we found that the integral converges well below this redshift since the merger rate decays rapidly at high redshift. Eqn. (12) implies that the distribution of triggering redshift given at redshift is:
(13) 
where again is tied to via the light curve model. Obviously only quasars with can contribute to .
The QSOtriggering rate is related to the halo merger rate by:
(14) 
where we parameterize the fraction as
(15) 
where describes how decreases towards lower redshift, i.e., the overall reduction due to the consumption of cold gas. We also introduce exponential cutoffs of at both high and low mass such that at each redshift, halos with too small cannot trigger QSO activity; on the other hand, overly massive halos cannot cool gas efficiently and BH growth halts, and the gasrich fraction may become progressively smaller at higher halo masses (especially at low redshift). We discuss choices for these parameters later in §3.
To proceed further we must specify the relation between and , and choose a light curve model. We assume that the correlation is a powerlaw with lognormal scatter, as motivated by some analytical arguments and hydrodynamical simulations (e.g., Wyithe & Loeb, 2002, 2003; Springel et al., 2005b; Hopkins et al., 2005; Lidz et al., 2006):
(16) 
where the mean relation is:
(17) 
with normalization and powerlaw slope . The log scatter around this mean relation is denoted as . When the scatter between and is incorporated, Eqn. (14) becomes:
from which we obtain the probability distribution of postmerger halo mass at fixed peak luminosity and redshift :
(19) 
We derive the probability distribution of host halo mass at given instantaneous luminosity and redshift, , as:
(20) 
which can be convolved with the halo linear bias to derive the QSO linear bias at instantaneous luminosity and redshift .
Since we have assumed that a second QSOtriggering merger event has not occurred, at redshift the postmerger halo should maintain most of its identity as it formed sometime earlier, i.e., we neglect mass added to the halo via minor mergers or accretion of diffuse dark matter between the halo formation time and the time when the QSO is observed at redshift .
Given the halo mass distribution (2.1), we can further derive the “active” halo mass function hosting a QSO with :
(21) 
which can be used to compute the halo duty cycles.
2.2. The Light Curve Model
For the light curve model there are several common choices in
mergerbased cosmological QSO models
(e.g., Kauffmann & Haehnelt, 2000; Wyithe & Loeb, 2002, 2003):
a) a light bulb model in which the QSO shines at a constant value
for a fixed time ; b) an exponential decay model
; c) a powerlaw decay model
with . Once
a light curve model is chosen, the total accreted mass during the
QSO phase is simply
(22) 
The evolution of the luminosity Eddington ratio is (neglecting the seed BH mass):
(23) 
where is the Eddington luminosity per .
Unfortunately none of the three LC models is selfconsistent without having at the very beginning of BH growth, simply because they neglect the rising part of the light curve. To remedy this, we must model the evolution of luminosity and BH growth selfconsistently (e.g., Small & Blandford, 1992; Yu & Lu, 2004, 2008). We consider a general form of light curve where the BH first grows exponentially at constant luminosity Eddington ratio (e.g., Salpeter, 1964) to at , and then the luminosity decays monotonically as a powerlaw (e.g., Yu & Lu, 2008):
(24) 
where is determined by:
(25) 
where is the seed BH mass at the triggering time . In eqn. (24), determines how rapidly the LC decays. Larger values lead to more rapid decay. Thus this model accommodates a broad range of possible light curves. With this light curve model (24) we have:
(26) 
where we require so that a BH cannot grow infinite mass. The evolution of the luminosity Eddington ratio is therefore:
(27) 
where during the decaying part of the LC the Eddington ratio monotonically decreases to zero. We assume that the seed BH mass is a fraction of the peak mass , and we have
(28) 
In what follows we set . Thus the seed BH is negligible compared to the total mass accreted, and where is the (folding) Salpeter timescale (Salpeter, 1964). Although is the formal definition of Eddingtonlimited accretion, it assumes the electron scattering cross section, and in practice superEddington accretion cannot be ruled out (e.g., Begelman, 2002). So we consider within the range .
The relic BH mass, given eqn. (26), is simply
(29) 
If instead we choose an exponential model for the decaying part of the LC (e.g., Kauffmann & Haehnelt, 2000), , then the relic mass becomes:
(30) 
Therefore LC models with large values of mimic the exponentially decaying model, and we do not consider the exponential LC model further.
If the decaying parameter is independent of mass, then Eqns. (26) and (29) imply . In other words, the scalings between and , and between and should be the same for QSOtriggering merger remnants (e.g., early type galaxies). Some theoretical models and simulations predict (e.g., King, 2003; Di Matteo et al., 2005; Springel et al., 2005a), where momentum is conserved in selfregulated feedback; while others invoking energy conservation predict (e.g., Silk & Rees, 1998; Wyithe & Loeb, 2003). Two recent determinations of the local black hole masshalo mass relation gave (Ferrarese, 2002) and (Baes et al., 2003). Within uncertainties seems to be a reasonable scaling, therefore we do not consider further mass dependence of . Given the definition of (Eqn. 28), our light curve model (24) is then universal in the sense that it scales with in a selfsimilar fashion.
As an example, Fig. 2 shows three light curves with for and , in which case . The time that a QSO spends above of its peak luminosity is typically for sharp decaying curves, but it can be substantially longer for extended decaying curves.
2.3. BH Mass and Eddington Ratio Distributions
Now we have specified the light curve model and the peak luminosityhalo mass relation, and tied the luminosity function to the QSOtriggering merger rate. We continue to derive the instantaneous BH mass and Eddington ratio distributions at fixed instantaneous luminosity and redshift , which can be compared with observationally determined distributions (e.g., Babić et al., 2007; Kollmeier et al., 2006; Shen et al., 2008a; Gavignaud et al., 2008).
Suppose that at redshift we observe quasars with instantaneous luminosity . These quasars consist of objects triggered at different earlier redshift and are at different stages of their evolution when witnessed at (e.g., Eqn. 12). There is a characteristic earlier redshift , determined by
(31) 
where is the cosmic time. Quasars triggered between
are all in the rising part of the LC; while quasars
triggered earlier than are all in the decaying part of the
LC. In order to contribute to , a quasar should have
peak luminosity if it is triggered right at
, and it should have otherwise. Therefore
the triggering redshift distribution (Eqn. 13)
peaks around
(32) 
where is given by Eqn. (13) and is
determined by the decaying half of the LC (Eqn. 27), i.e., the triggering redshift
is a monotonically increasing function of
(note that here refers to the
instantaneous BH mass observed at ). The denominator in the
above equation is to normalize the distribution – we literally
augment the distribution with the pileup of objects with
triggered within . The
fraction of such objects is about for and becomes negligible at lower
luminosities. This procedure removes the spike ( function)
at the constant BH mass in the
distribution, which is an artifact of our model
3. The Reference Model
We are now ready to convolve the LC model with the QSOtriggering rate Eqn. (2.1) to predict the LF using Eqn. (12). Before we continue, let us review the model parameters and available observational/theoretical constraints.

QSO triggering rate [, , , ]: We choose our fiducial value of following the traditional definition of major mergers ( or ). For the minimum halo mass below which efficient accretion onto the BH is hampered we simply choose . We model the maximum halo mass cut above which the gasrich merger fraction drops as a function of redshift:
(34) with . Therefore at lower redshift halos have a smaller upper threshold. Since at rich group to cluster size halos no longer host luminous QSOs we tentatively set . The parameters and control the the shape of the LF at both the faint and bright luminosity ends. The global gasrich merger fraction, , is more challenging to determine. Direct observations of the cold gas fraction as function of redshift and stellar mass (and therefore halo mass) are still limited by large uncertainties. Moreover, how gasrich a merger needs to be in order to trigger efficient BH accretion is not well constrained. For these reasons, we simply model as a twopiece function such that when , the peak of the bright quasar population; and linearly decreases to at . We note again that should be regarded as the fraction of major mergers that trigger QSO activity.

The relation (, , ): some merger event simulations (e.g., Springel et al., 2005a; Hopkins et al., 2005; Lidz et al., 2006) reveal a correlation between the peak luminosity and the mass of the (postmerger) halo: at , with a lognormal scatter dex. Other analytical arguments predict , where energy is conserved during BH feedback (e.g., Silk & Rees, 1998; Wyithe & Loeb, 2003). We choose values between since the above simulations are also consistent with the slope. Furthermore, we allow the normalization parameter to evolve with redshift, with . Thus for the same peak luminosity, the host halos become less massive at higher redshift. This decrease of the characteristic halo mass with redshift is also modeled by several other authors, although we do not restrict to (e.g., Lapi et al., 2006), or more complicated prescriptions (e.g., Wyithe & Loeb, 2003; Croton, 2009), since the form of the (or ) relation is also different in various prescriptions. A smaller characteristic halo mass is easier to account for the QSO abundance at high redshift since there are more mergers of smaller halos, but it also reduces the clustering strength. The intrinsic scatter of the relation, , will have effects on both the luminosity function and clustering. Larger values of lead to higher QSO counts, but will also dilute the clustering strength due to the upscattering of lower mass halos (e.g., White et al., 2008). The relation establishes a baseline for the mapping from halos to SMBHs. Although in our fiducial model we adopt the above rather simple scaling [] for the redshift evolution in , we will discuss alternative parameterizations for in §4.2 and §4.4.

The light curve model (, , , ): our chosen fiducial value for is