Supermassive Black Holes in the Hierarchical Universe: A General Framework and Observational Tests


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 self-consistently: an initial exponential growth at constant Eddington ratio of order unity until it reaches the peak luminosity, followed by a power-law decay. Assuming that the peak luminosity correlates with the post-merger 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 X-ray 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 log-normal, 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 – large-scale structure of universe – quasars: general – surveys

1. 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 phase2 (e.g., Salpeter, 1964; Zel’dovich & Novikov, 1964; Lynden-Bell, 1969; Soltan, 1982; Small & Blandford, 1992; Salucci et al., 1999; Yu & Tremaine, 2002; Shankar et al., 2004; Marconi et al., 2004). The statistics of the local dormant SMBHs and distant active QSOs therefore provide clues to the build-up of the local BH density across cosmic time, and have implications for the hierarchical structure formation paradigm, as well as for the interactions between accreting SMBHs and galaxy bulges.

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. Gas-rich 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 ULIRG-quasar 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 small-scale 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 merger-based 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 merger-based 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 friends-of-friends (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 X-ray 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; Lynden-Bell, 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:


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 equation3 (cf. Small & Blandford, 1992):


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 X-ray band (e.g., Ueda et al., 2003; Hasinger et al., 2005; Silverman et al., 2005; Barger et al., 2005). While wide-field optical surveys provide the best constraints on the bright end of the QSO luminosity functions, deep X-ray surveys can probe the obscured faint end AGN population. Although it has been known for a while that the spatial density of optically-selected bright quasars peaks at redshift , the spatial density of fainter AGNs selected in the X-ray 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 so-called downsizing of the cosmic SMBH assembly.

Hopkins et al. (2007) compiled QSO luminosity function data from surveys in various bands (optical, X-ray, 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 large-scale 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 two-point 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:


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 Press-Schechter (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 look-back 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:


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 friends-of-friends 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):


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 semi-analytical 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; Boylan-Kolchin 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 stripping4 and hence the duration of dynamical friction is greatly extended.

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; Boylan-Kolchin 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).

Figure 1.— Merger timescales and (sub)halo merger rates. Upper: comparison of the subhalo merger timescale from three dynamical friction prescriptions. Bottom: the redshift evolution of the halo merger rate (black solid lines), and the subhalo merger rate using the Boylan-Kolchin et al. (2008) (black dashed lines), and Jiang et al. (2008) (red dashed lines) prescriptions for the merger timescale, for three (post)merger halo masses.

The fitting formula for in Jiang et al. (2008), who used hydro/-body simulations, has the following form:


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 Boylan-Kolchin et al. (2008) is given by:


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 :


where is a function of and is determined by:


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 :


for the fitting formula of in Jiang et al. (2008), or


for the fitting formula of in Boylan-Kolchin 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 Boylan-Kolchin et al. (2008) (solid line), Jiang et al. (2008) (dashed line), and a commonly-used SAM prescription: (dotted line). For the major merger regime the fitting formula in Boylan-Kolchin 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 Boylan-Kolchin 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 Boylan-Kolchin 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 QSO-triggering 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
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
Power-law slope of the decaying phase
The QSO-triggering rate
Minimum mass ratio 0.25
Exponential lower mass cut
Exponential upper mass cut

Note. – Parameter values are for the fiducial model.

Table 1Notation and Model Parameters

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 semi-analytical 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 gas-rich major merger event ().

  • The build-up 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 gas-rich major merger is triggered and self-regulated 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 mass-dependent 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 luminosity5:


where is the QSO number density in the luminosity range at redshift ; is the QSO-triggering 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:


where again is tied to via the light curve model. Obviously only quasars with can contribute to .

The QSO-triggering rate is related to the halo merger rate by:


where we parameterize the fraction as


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 gas-rich 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 power-law with log-normal 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):


where the mean relation is:


with normalization and power-law 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 :


We derive the probability distribution of host halo mass at given instantaneous luminosity and redshift, , as:


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 QSO-triggering 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 :


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 merger-based 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 power-law decay model with . Once a light curve model is chosen, the total accreted mass during the QSO phase is simply6:


The evolution of the luminosity Eddington ratio is (neglecting the seed BH mass):


where is the Eddington luminosity per .

Unfortunately none of the three LC models is self-consistent 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 self-consistently (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 power-law (e.g., Yu & Lu, 2008):


where is determined by:


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:


where we require so that a BH cannot grow infinite mass. The evolution of the luminosity Eddington ratio is therefore:


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


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 Eddington-limited accretion, it assumes the electron scattering cross section, and in practice super-Eddington accretion cannot be ruled out (e.g., Begelman, 2002). So we consider within the range .

The relic BH mass, given eqn. (26), is simply


If instead we choose an exponential model for the decaying part of the LC (e.g., Kauffmann & Haehnelt, 2000), , then the relic mass becomes:


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 QSO-triggering 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 self-regulated feedback; while others invoking energy conservation predict (e.g., Silk & Rees, 1998; Wyithe & Loeb, 2003). Two recent determinations of the local black hole mass-halo 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 self-similar fashion.

Figure 2.— Example light curves for and , for . Larger values of lead to more rapid decay. The time that the QSO spends above half of its peak luminosity is short, Myr for large values of .

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

Figure 3.— Predicted bolometric luminosity functions in our fiducial model (solid lines). Overplotted are the complied bolometric LF data from Hopkins et al. (2007) for optical (black circles), soft X-ray (green squares) and hard X-ray (red triangles) samples. The dashed lines are the predicted LF for a model where the normalization in the mean relation is reduced by an additional amount at , tuned to fit the LF at (see §4.2 for details). The blue dotted lines are the predictions based on an alternative prescription for the redshift evolution in the normalization of the relation as discussed in §4.4.

Now we have specified the light curve model and the peak luminosity-halo mass relation, and tied the luminosity function to the QSO-triggering 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


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 around7 since decreases rapidly when increases. Moreover, at the observing redshift , all the contributing quasars triggered between will have Eddington ratio and BH mass ; while all the contributing quasars triggered before will have Eddington ratio and BH mass . The probability distribution of instantaneous BH mass at the observing redshift and luminosity is therefore:


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 model8. Eqn. (32) gives the equivalent distribution of Eddington ratios at and at instantaneous luminosity .

Given the luminosity function (12) and the distribution of (32) we can determine the active BH mass function (i.e., the mass function of active QSOs above some luminosity cut ) as


3. The Reference Model

Figure 4.— Distributions of triggering redshift (upper panels) and halo mass (lower panels) in our fiducial model, for three instantaneous luminosities and two values of observing redshifts. The distributions of triggering redshift generally peak around (marked by arrows) given by Eqn. (31). For lower luminosity, the distribution of the triggering redshift is more extended. As a consequence, there is a significant population of faded low-luminosity QSOs which have massive hosts.

We are now ready to convolve the LC model with the QSO-triggering 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 gas-rich merger fraction drops as a function of redshift:


    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 gas-rich 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 gas-rich a merger needs to be in order to trigger efficient BH accretion is not well constrained. For these reasons, we simply model as a two-piece 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 up-scattering 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 (