AGN evolution from a galaxy evolution viewpoint

AGN evolution from a galaxy evolution viewpoint

Neven Caplar Simon J. Lilly Benny Trakhtenbrot11affiliation: Zwicky Fellow Institute for Astronomy, Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093, Zurich, Switzerland

We explore the connections between the evolving galaxy and AGN populations. We present a simple phenomenological model that links the evolving galaxy mass function and the evolving quasar luminosity function, which makes specific and testable predictions for the distribution of host galaxy masses for AGN of different luminosities. We show that the normalisations of the galaxy mass function and of the AGN luminosity function closely track each other over a wide range of redshifts, implying a constant “duty cycle” of AGN activity. The strong redshift evolution in the AGN can be produced by either an evolution in the distribution of Eddington ratios, or in the mass ratio, or both. To try to break this degeneracy we look at the distribution of AGN in the SDSS () plane, showing that an evolving ratio reproduces the observed data and also reproduces the local relations which connect the black hole population with the host galaxies for both quenched and star-forming populations. We stress that observational studies that compare the masses of black holes in active galaxies at high redshift with those in quiescent galaxies locally will always see much weaker evolution. Evolution of this form would produce, or could be produced by, a redshift-independent relation and could explain why the local relation is tighter than even if is not directly linked to black hole growth. Irrespective of the evolution of , the model reproduces both the appearance of “downsizing” and the so-called “sub-Eddington boundary” without any mass-dependence in the evolution of black hole growth rates.

Subject headings:
galaxies: active - galaxies: evolution - galaxies: mass function - quasars: general - quasars: supermassive black holes

July 12, 2019

1. Introduction

Over the last few years there has been a lot of interest in the cosmic “co-evolution” of supermassive black holes (SMBH) and the stellar populations of the galaxies that they reside in. This has been motivated on the one hand by the tight scaling relations that have been established between the masses of the black holes and various parameters of the host galaxies (e.g. Kormendy & Richstone, 1995; Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002; Marconi & Hunt, 2003; Häring & Rix, 2004; Ferrarese & Ford, 2005; Greene & Ho, 2006; Greene & Ho, 2007; Aller & Richstone, 2007; Greene et al., 2008; Graham et al., 2011; Sani et al., 2011; Vika et al., 2012; McConnell & Ma, 2013) and on the other by the overall similarities between the evolution of the star-formation rate density of the Universe and the luminosity density that is ascribed to SMBH accretion (e.g. Boyle & Terlevich, 1998; Franceschini et al., 1999; Marconi et al., 2004; Heckman et al., 2004; Hasinger et al., 2005; Silverman et al., 2009; Zheng et al., 2009). This paper seeks to link directly the evolution of the galaxy and AGN populations via a simple phenomenological model.

Looking first at galaxies, almost all galaxies can be broadly classified into a few distinct populations. The majority of star-forming (SF) galaxies have a star-formation rate (SFR) that is strongly correlated to their existing stellar mass, producing the so-called “Main Sequence” in which the specific SFR (or sSFR) varies only weakly with stellar mass. This characteristic sSFR of the Main Sequence however increases strongly with look-back time (Daddi et al., 2007; Elbaz et al., 2007; Pannella et al., 2009; Speagle et al., 2014) and is about a factor of twenty higher at compared with the present day value. A small percentage of star-forming galaxies have significantly elevated sSFR, above the Main Sequence. It appears that the fraction of these “outliers” is more or less constant to z 2 (Sargent et al., 2012) and that they represent of order 10% of the integrated star-formation that is occurring at any epoch. There is also a large population of “quenched” galaxies in which the star formation is substantially lower than on the Main Sequence, producing an sSFR that is much lower than the inverse Hubble time. Our understanding of the physical processes that lead to the quenching of star-forming galaxies is still quite limited and a number of plausible physical mechanisms have been proposed (some of which involve AGN directly). However, the main empirical, or phenomenological, features of this quenching process are quite well understood based on the characteristic of the evolving population(s) of galaxies.

As the available data on the galaxy population at substantial look-back times has improved, new analysis techniques have been introduced that take a purely empirical (phenomenological) approach to the data, see e.g Peng et al. (2010) and Behroozi et al. (2013a). These are complementary to the semi-analytic models of galaxy evolution (e.g. Somerville et al., 2008, Conroy & Wechsler, 2009, Henriques et al., 2014).

The approach in Peng et al. (2010) was to identify a few striking simplicities exhibited by the galaxy population(s) and to explore the consequences of these, where possible analytically, in terms of the most basic continuity equations linking the galaxy population(s) at different epochs. The Peng et al. (2010) analysis was based on dividing the galaxy population into two components, the star-forming Main Sequence (including the outliers) and the quenched population of passive galaxies. Much of the Peng et al. (2010) formalism is based on the observation that the characteristic Schechter of the mass-function of the star-forming population has been more for less constant back to at least (and likely to ) despite the substantial increase in stellar mass (by a factor of 10-30) of any galaxy that stays on the star-forming Main Sequence over this same time period (Bell et al., 2003; Bell et al., 2007; Ilbert et al., 2010; Pozzetti et al., 2010; Ilbert et al., 2013). The Peng et al. (2010) continuity formalism is very successful at reproducing the single and double Schechter (Press & Schechter, 1974) shapes of the mass functions of star-forming and passive galaxies in SDSS and also explains the quantitative relations between the Schechter parameters of these mass functions which constitute a test of the approach. The formalism allows easy computation of things like the quenching rate of galaxies, the mass function of galaxies that are undergoing quenching at any epoch, and so on. The alternative phenomenological approach of “abundance matching” (Behroozi et al., 2013a) provides similar results, in terms of mass functions and star formation histories. With these recent developments, we now have a self-consistent, empirical (”phenomenological”) description of the evolving galaxy population at least back to z 4 in terms of the evolving mass-functions of both the star-forming and quenched populations of galaxies.

Turning to the AGN, the most basic description of the evolving population is the bolometric luminosity function (i.e. quasar luminosity function, QLF), i.e. . Large homogeneous samples of AGN have been created, from optical surveys such as Sloan Digital Sky Survey (SDSS) (Schneider et al., 2010; Ross et al., 2013) and deep X-ray surveys, out to redshifts (e.g. Hasinger et al., 2005; Brown et al., 2006; Silverman et al., 2008; Brusa et al., 2010; Civano et al., 2011; Kalfountzou et al., 2014; Ueda et al., 2014). These QLF studies have found that the QLF is best described by a double power law, i.e. two power-law segments broken at a characteristic luminosity ,


These studies tend to agree that the increase in the number density of luminous quasars with increasing redshift is mostly driven by an evolution of at redshifts below . The exact shape of the QLF and need for evolution of other parameters is still debated (Ueda et al., 2003; Barger et al., 2005; Croom et al., 2009; Aird et al., 2010; Assef et al., 2011).

The double power-law shape of the QLF contrasts the galaxy mass function which is clearly better described by one or more Schechter functions, i.e. a power-law at low masses (or luminosities) and an exponential cut-off at masses above a characteristic mass . This difference in the shapes of these two most basic descriptions of the two populations is interesting in that AGN are being harboured in galaxies and one might naively expect similarities between these two functions. This difference will form an important part of the current analysis.

The luminosity of an individual AGN can be expressed as


where is given in erg s, is the Eddington ratio and is the mass of the central black hole in units of solar mass. The SMBH mass is therefore a key quantity in understanding both an individual AGN and the AGN population. Unfortunately, black hole masses can be measured reliably with dynamical modeling only for small number of nearby sources in which the SMBH is generally quiescent (e.g. Gültekin et al., 2009; Schulze & Gebhardt, 2011; Graham & Scott, 2013; McConnell & Ma, 2013; Rusli et al., 2013; Kormendy & Ho, 2013 and references within). Objects that are actively accreting and/or that are further away need to be analyzed with different techniques. For AGNs that show broad lines in their spectra it is possible to determine by conducting reverberation mapping campaigns (see review by Peterson, 2013 and references within). The results of such campaigns makes it possible to construct “single-epoch” mass estimators, which allow the determination of of AGN based on luminosity and line width. Systematic uncertainties in these estimators are calibrated using a sample of local AGNs for which black hole masses are statistically known from the of local (inactive) galaxies.

Early suggestions that the distribution of specific accretion rates onto black holes has no dependence on black hole or galaxy mass (Yu et al., 2005; Kollmeier et al., 2006; Merloni & Heinz, 2008), have gained further observational support from analysis of the PRIMUS and COSMOS survey fields (Aird et al., 2012, Bongiorno et al., 2012). We will adopt a similar assumption in this paper.

Observations in the local Universe have revealed several interesting correlations between the mass of the SMBH in the center of a given galaxy and various quantities describing the surrounding galaxy. Specifically, quite tight relations have been established with the stellar velocity dispersion () and with the stellar mass of the bulge in the galaxy () (Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Marconi & Hunt, 2003). Canonical values of the ratio between and are between and (Häring & Rix, 2004). The recent extensive review by Kormendy & Ho (2013) has argued for larger values, up to , mainly due to differentiation between bulges and pseudo-bulges, by omitting mergers in progress and by the inclusion of dark matter into dynamical modelling. Other scaling relations have been proposed, such as where is Sersic index describing the light profile of the galaxy, where is mass of the dark matter halo and where is integrated stellar mass of the galaxy (Sani et al., 2011). All these relations tend to show more scatter, but the basic difficulty is that most of these galaxy parameters are strongly correlated with each other. It has been claimed that the correlation between is more fundamental then relation (Marleau et al., 2013), especially in star-forming hosts, which could then help to explain the AGN sources in bulgeless, pure disk, galaxies that have been observed (Simmons et al., 2013). In this paper we will base the analysis on simply because the available mass functions of galaxies generally utilise the integrated stellar mass rather than the bulge mass, especially at high redshift.

In fact, almost all of the properties of galaxies that host AGN are still widely debated. Although there are undoubtedly some active AGN found in quenched galaxies, the bulk of radiatively efficient AGN seem to reside in star-forming galaxies (Netzer, 2009; Silverman et al., 2009; Schawinski et al., 2010; Koss et al., 2011; Cimatti et al., 2013; Rosario et al., 2013; Matsuoka et al., 2014), especially in those relatively low luminosity systems where the host galaxy can be most easily discerned. Nevertheless, it is still not entirely clear what fraction of AGN are situated in quenched galaxies and whether this fraction changes with AGN luminosity.

Many studies have suggested that there is some redshift evolution in mass scaling relations (Peng et al., 2006; Decarli et al., 2010; Merloni et al., 2010; Trakhtenbrot & Netzer, 2010; Bennert et al., 2011; Sijacki et al., 2014), but others have found no significant evolution (Jahnke et al., 2009; Cisternas et al., 2011; Mullaney et al., 2012; Schramm & Silverman, 2013) or have argued that an observed evolution can be fully explained with selection effects (Schulze & Gebhardt, 2011; Schulze & Wisotzki, 2014). The redshift evolution of the scaling relations is still debated and we will explore different possibilities in the current paper within our model framework. We will also emphsize the methodological issues associated with the choice of samples.

Several studies have also attempted to constrain different aspects of black hole evolution with empirical or semi-empirical approaches which model the evolving black hole mass function, accretion rate, QLF, duty cycle, accretion efficiency or some subset of these quantities (Merloni, 2004; Shankar et al., 2004; Yu & Lu, 2004; Merloni & Heinz, 2008; Hopkins & Hernquist, 2009; Shankar et al., 2009; Shen, 2009; Cao, 2010; Li et al., 2011; Steinhardt et al., 2011; Shankar et al., 2013b; Conroy & White, 2013; Novak, 2013; Goulding et al., 2014).

In this work, we aim to apply a similar phenomenological approach that has proved so successful in describing the evolving galaxy population evolution to the evolving AGN population. We will take the simplest observables of the population and try to infer the underlying simplicities of the situation and use straightforward prescriptions to construct a simple model that can both explain the salient observational data and which can be used to make simple testable predictions for other quantities. A focus of this work is to try to use our improved understanding of the evolving galaxy population, and especially of the evolving mass function of galaxies to high redshifts, to interpret the evolution of the AGN population. We aim thereby to create a simple global model to interpret the evolving AGN population and to enable us to evaluate biasses that may arise in observational work, e.g. through the use of luminosity-selected samples. The model is based on the construction of the AGN QLF via a double convolution of the underlying host galaxy mass function with a scatter function representing the black-hole to stellar mass ratio, and with an Eddington ratio distribution.

Hopkins et al. (2008b) and Hopkins & Quataert (2010) have also used the observed galaxy mass functions as a basis for modelling the AGN population using AGN lightcurves obtained from simulations. A more similar approach to ours has been used by Aird et al. (2013) with an observationally driven model that connects AGN and galaxy evolution out to z 1, based on observational data from the PRIMUS survey (Aird et al., 2012). They have shown that is possible to reproduce the luminosity function measurements with a mass-independent distribution of specific accretion rates, i.e. linking the AGN QLF to the galaxy mass function. Hickox et al. (2014) has also used a simple phenomenological approach and shown that the main observational trends (such as QLF and SFR- correlations) can be recovered by assuming that all star-forming galaxies host an AGN and that star formation and black hole accretion are correlated on 100 Myr time-scales, i.e. in this case they were linking the growth of the black holes and stellar populations. A very recent paper of Veale et al. (2014) has, like the current work, explored the links between the evolving galaxy (and halo) mass functions and the observed AGN QLF via a convolution approach. Veale et al. (2014) have emphasised the degeneracies between the “mass function” and “mass growth” approaches when only the AGN QLF is considered. In the current work, we incorporate other observational data, most notably the distribution of AGN in the plane, to move beyond the information in the QLF alone.

The layout of the paper is as follows. In Section 2 we first state the three simple Ansätze that are used to construct the model. We then show in Section 3 how the quasar can be simply constructed from the galaxy and Eddington ratio distribution via a convolution, show how the parameters of these distribution functions are connected, and make testable predictions of the mass distribution of the host galaxies of quasars selected at different luminosities. In Section 4 we review the observed epoch-dependent and functions and then in Section 5 we compare these within the framework of the convolution model to derive interesting conclusions about the duty cycle of AGN.

Up until this point, the model is completely general. In Section 6, we then consider different possibilities for the evolution of the ratio and show that one particular choice can explain the observed distribution of luminous quasars in the ( plane and local scaling relations of black holes in both star-forming and passive galaxies and the possible differences in evolution of active and passive systems. Section 7 presents a discussion of some further implications of both the general model and the specific implementation, differences between and the appearance of downsizing within the AGN population, and the pervasive biasses that enter into luminosity-selected samples. The paper then concludes with a summary section.

Throughout the Paper, we will assume a CDM cosmology, with parameters =0.3, =0.7 and = 70 km s Mpc. Luminosities will be given in units of erg s and will refer to bolometric luminosities, unless specified otherwise. We use the term ”dex” to denote the antilogarithm, i.e. dex = 10. We also define all distribution functions, i.e. the star-forming galaxy mass function , the associated star-forming galaxy black hole mass function , the AGN luminosity function and the probability distribution of Eddington ratio , in log space. This leads to power-law exponents that differ by unity relative to distribution functions defined in linear space. Therefore, the units of , and are Mpc dex.

2. Ansätze

The essence of this phenomenological approach is to make a limited number of simple Ansätze that allow us to construct a model using analytic techniques, or very elementary numerical modeling, based on straightforward representations of the most important features of the observational data. These Ansätze are in a sense “assumptions” and a decisive observational disproof of any of them would largely invalidate the model. Clearly they are very unlikely to be exactly true. Their value, as the basis for a simple “toy model”, is that we believe they are likely to be broadly true.

The three Ansätze in the current work are as follows:

  • radiatively efficient AGNs are found in star-forming galaxies,

  • the probability distribution of the Eddington ratio does not depend on the black hole mass of the system,

  • the mass of the central black hole is linked to the stellar mass (), with some scatter, and we will for simplicity set .

The justification for the first and the third of these has been reviewed in the Introduction and the second is closely related to the assumption of Aird et al. (2013). We will justify these, and the choice of , further below. Of course, there are a number of other implicit assumptions that are being made: observational data is not wildly wrong, individual black hole and stellar mass measurements are not systematically biased, the cosmological model is more or less correct and so on. We will not consider these further.

3. The origin of the quasar luminosity function

With our Ansätze, we can use our knowledge of the mass function of galaxies to construct the black hole mass function. To do this, we start with the mass function of star-forming galaxies and then impose a black hole to stellar mass ratio , with an additional log-normal scatter of , to create a black hole mass function in star-forming hosts. This will serve as the basis for the radiatively efficient AGN population. In first part of the paper, we will set the scatter for initial analytic simplicity, before reintroducing it with when we evaluate the model numerically.

If radiatively efficient accretion onto central black holes is only occurring in star-forming hosts, the quasar luminosity function can be created from a convolution of the black hole mass function in star-forming galaxies with a probability distribution of the Eddington ratio ,


where is the resulting QLF, and is the probability distribution of the Eddington ratio, , as defined in Equation (2).

3.1. Convolution of Schechter functions with power law Eddington ratios

As noted above in Equation (3), the QLF will, with our Ansätze, be a convolution of the black hole mass function (itself derived from the stellar mass function of star-forming galaxies ) and the distribution of Eddington ratios , which gives the distribution of luminosities for black holes of a given mass. In this section we look at the general features of the distribution that are needed to produce a double power-law QLF from a Schechter-like mass function and derive the connections that will exist between the parameters of these different functions.

Figure 1.— Schematic representation of our procedure to create the QLF. Starting from the star-forming mass function (leftmost) we use an scaling to get the mass function of SMBH in star-forming galaxies, with scatter as required (2 panel). We convolve it with the Eddington ratio distribution (in the 3 panel) to create QLF in the rightmost panel. Blue and red Eddington ratio distributions differ in the choice of parametrization of low end slope. The faint end slope of QLF will be same as low mass slope of galaxy mass function or low end slope of Eddington ratio distribution, depending on relative steepness of these two slopes. A short summary of the connections between parameters of the QLF and parameters of the contributing functions is given in Equation (12).

First, we immediately note that in a model with a mass-independent and an input Schechter mass function, the only way to produce a power law at the bright end of the QLF (with slope ) is to have a that is also a power law of the same slope at high values of the Eddington ratio. This power-law will then ensure that there is a high-end power-law in the QLF even as the mass function drops off exponentially. This is shown in the Appendix. Given the limited dynamic range of the data, other representations of the QLF instead of a double power law are possible, e.g. a log-normal bright end or a modified Schechter function (e.g. Hopkins et al., 2007, Aird et al., 2012, Veale et al., 2014). We choose the conventional double power law for analytic simplicity and because it certainly provides a reasonable representation of the available data. Denoting the slope of the high end of Eddington ratio distribution by we can thus equate


with the bright end slope of the QLF.

At the faint end, we may also expect a power-law QLF but now the QLF faint end slope will be given by the steeper of the low end slope of and the low end slope of , which we denote by (see also discussion in Veale et al. (2014)). Figure 1 illustrates what is happening with faint end of QLF (with slope) for different Eddington ratio assumptions.

We can express this as


The natural conclusion of a mass-independent Eddington rate distribution is that the faint end of QLF is set up by either the low end slope of underlying black hole mass function or by the low end slope of the Eddington ratio function. If the logarithmic slope of the vs. relation is as above, then the faint end slope of the black hole mass function will be related to that of the star-forming galaxies by .

We note at this point that there is good observational evidence that the faint end QLF slope, , is similar to the observed low mass slope of the mass function of star-forming galaxies, , allowing for the reversal of sign in our definition. The QLF faint end slope, is usually observed to be between 0.3 and 0.9 (Hasinger et al., 2005; Hopkins et al., 2007; Aird et al., 2010; Masters et al., 2012; Ueda et al., 2014). On the other hand the observed values of range from to -0.6 to -0.2 (Baldry et al., 2008; Pérez-González et al., 2008; Peng et al., 2010; González et al., 2011; Kajisawa et al., 2011; Baldry et al., 2012; Lee et al., 2012) with most newer studies converging around .

The data are therefore consistent with the idea that . This suggests that is indeed being set by the low end slope of the black hole mass function (and not the low end of the Eddington ration distribution) and further that . We will henceforth assume that this is the case, i.e. that . This then means that the low end slope of can take any value that is shallower than this, i.e. . The high and low end power-laws of the Eddington ratio distribution will meet in a knee at a characteristic Eddington ratio which we denote by at which point the value of has a characteristic value .

3.2. The simplest Eddington ratio distribution that reproduces the shape of the QLF

To further demonstrate our approach we use the simplest function possible for the Eddington ratio distribution that reproduces, analytically, the required broken power law shape of QLF. We use a ”triangular” distribution in which some fraction, , of objects are active above a certain threshold value with an Eddington ratio distribution that is a single power law with slope , while (1-) of objects are completely inactive. Effectively this sets the low end slope .

The exact functional form of can then be written as


where in constrained by requirement that all of the objects have to be either active or inactive, so in this case has to integrate to a duty cycle .

In an Appendix we show that the knee of the QLF, , will be related to the parameters of the original distribution functions through a formula describing the population which is analogous to Equation (2) of individual black holes, i.e.


where the factor denotes a small multiplicative factor that is weakly dependent on and varying by less than 0.15 dex.

In the same Appendix, we also show that the normalization of the QLF, i.e. the value of at will be linked to the normalization of the galaxy mass function and the normalization of the Eddington ratio distribution, ,


where the denotes once again a small multiplicative factor, weakly dependent on and varying by less then 0.15 dex. This is not surprising and is stating the fact that the normalization of objects at given luminosity is connected with the normalization of the contributing functions at and . Simply put, if there are more objects at mass and more of them are active at , then we expect also more objects with corresponding luminosity .

3.3. Broken power law Eddington ratio distribution and generalized duty cycle

Even though the triangular distribution of is a simple and analytically tractable one, it is unlikely to describe the real AGN population. In the remainder of this paper we will adopt a more realistic broken power law distribution of Eddington ratio given by


where we set =0 for further analysis. This function fits both observations and hydrodynamical simulations better, which show no sudden cutoff in the Eddington ratio distribution at some single value (e.g. Novak et al., 2011, Kelly & Shen, 2013). This distribution diverges logarithmically at the low end. Since the integral of will therefore reach unity at some low value of , all black holes are “active” at some very low level and we need no longer consider “inactive” ones.

As discussed in Section 3.1, this distribution will also reproduce the same shape of the QLF as the ”triangular” distribution just discussed, provided that .

We note that such a distribution, (with sharp break at , instead of ”smooth” version above), would naturally arise if individual AGN are boosted to some initial Eddington ratio above (the distribution of which is given by ) and thereafter decay exponentially with a constant decay time . In this case, the concept of ”duty cycle” in the previous section is replaced by the value of . If the “boost plus decay” scenario is relevant, then it is easy to see that, if AGN are activated at some fractional rate per star-forming galaxy per unit time (e.g. Hopkins et al., 2005), then


so that is the fraction of black holes that are activated in the time interval corresponding to their subsequent decay time. This is quite a useful definition of duty cycle.

Of course, other scenarios could also produce this distribution and a more general “duty cycle” could be defined as the fraction of black holes accreting above (as in the triangular distribution above). In this particular case, as shown in the Appendix, this definition of duty cycle would be given by


The point is that the value of is a good indicator of a generalized duty cycle.

If , then the exact choice below is of no great importance to the convolution model and the connections between the parameters, which we can derive analytically in the simplified model as in Equations (4), (5), (7) and (8) will still hold. Specifically, we can write


where we have now also explicitly shown a dependence in the factors, because this corrective factor will also depend on our exact choice of low slope of the Eddington ratio distribution. Now using the relation, we can also connect the observed of the QLF to the of the star-forming galaxy mass function with


3.4. Predictions of the convolution model

Figure 2.— Expected shape of the black hole and the black hole host galaxy mass function, selected in different AGN luminosity bins. Mass are ploted relative to the Schechter of the galaxy population () and relative to the of the galaxy population. All luminosities are relative to the of the AGN luminosity function. The black dashed line is connecting the masses where the mass functions reach maximum value.

If the QLF is indeed produced by the convolution described in the previous two sections, then we can immediately see in general terms how the mass distribution of AGN, or host galaxies of AGN, selected at a given luminosity will vary with that selection-luminosity. This is shown in Figure 2 where we show the mass functions of galaxies and black holes. The masses and number densities are normalized respectively by the and of the input Schechter mass function of the star-forming galaxies. These distributions are plotted for different AGN luminosities relative to the knee luminosity in the quasar luminosity functions (see also Equation (1)). To generate this plot we used the model for considered in the previous section, i.e. with a flat below and with 0.5 dex scatter in the black hole to stellar mass relation. With this normalization the figures applies at all redshifts. Several points should be noted on this Figure.

First, all of the mass functions (of both black holes and host galaxies) always show a peak at some mass. This is quite unlike the input mass function of star-forming galaxies, and thus also of black holes, which increase monotonically to lower masses. The peak arises because the low mass end of the input mass function is modulated by the high part of , simply because it is the highest Eddington ratios that pull the lowest mass black holes (and hosts) into a given luminosity bin. The low mass end of these mass distributions should therefore have a slope of , i.e. with and (see below), we predict a low mass slope of +1.6. This is a quite generic and robust prediction of the convolution model and is independent of the choice of , the low behavior of discussed above. This is because for most luminosities, the high mass end of the mass functions in Figure 2 is dominated by the exponential fall-off of the input mass function.

Second, at high luminosities, above , the effect of the AGN luminosity on the host galaxy mass function is to change the numbers of hosts, but not to change the peak mass or, to first order, the shape of the mass function. This is because of the steep exponential fall off in the input mass function of galaxies above . Above , the peak host galaxy mass is always close to because there are so few more massive galaxies. The numbers of galaxies at the peak is determined by how high Eddington ratios are required to bring these objects into the luminosity range in question. For these high AGN luminosities, the shape of the mass function of (star-forming) hosts is very similar to that of passive galaxies except that the faint end slope is significantly steeper (, with ) that of the passive galaxies which is given by (see Peng et al. (2010)). In mathematical form, the mass function of galaxies at any AGN luminosity above will have a Schechter shape


At lower luminosities, well below , the behaviour changes and a region of intermediate slope appears. For black hole masses between and , the mass function will have the slope given by . At very low masses we will again always recover the slope of .

The location of the peak in the black hole mass function therefore depends on the slopes of both underlying distributions (Eddington ratio and galaxy mass function), i.e. on the sign of . Not surprisingly, this intermediate zone also appears in the host galaxy mass function with the same slope. This produces a peak host galaxy mass of


This difference in behaviour can be understood as follows: In the case the dominant contribution to the luminosity function will come from low mass objects accretion at high Eddington ratio which are more numerous than high mass objects accreting at low Eddington ratios. In the second case, roles are simply reversed and main contributor to the QLF will be high mass AGNs with low Eddington ratios.

If black hole masses were very tightly correlated with host galaxy masses, then clearly the same behavior would be seen in the black hole mass function(s) since these would always be a simple (mass-)scaling of the galaxy mass function(s). However, the effect of scatter in the black hole to host galaxy mass relation is quite marked. This is visible in Figure 2. The reason is clear: the mass function of black holes will now be much smoother than that of the host galaxies since the log normal scatter smooths out the sharp exponential drop of the galaxies at high masses (recall that Figure 2 used a 0.5 dex scatter). Of course, it might be argued that the black hole mass is somehow more fundamental, and that the galaxy mass function should be obtained by smoothing it. However, we know the galaxy mass function quite well, and it has an exponential cutoff (see Peng et al., 2010, Baldry et al., 2012). We do not, at this stage, know the underlying mass function of black holes in star-forming galaxies. With this scatter, there is a much smoother variation of the peak mass with AGN luminosity.

The dashed lines show the variation of the peak mass with luminosity. The fact that the difference between these dashed curves on the two panels decreases with luminosity illustrates the well known “Lauer effect” (Lauer et al., 2007) whereby the typical objects at high AGN luminosities will generally have been scattered above the mean black hole host mass relation. It should be noted, however, that there is no change in the low mass slope of the mass function(s) due to scatter and this remains a very robust prediction of our model.

It is worth noting in Figure 2 that the peak of the black hole mass distribution varies quite weakly with luminosity, i.e. a change in luminosity of 1 dex is associated with a smaller increase in the peak . This means that we will generally see higher Eddington ratios in higher luminosity quasars even though the Eddington ratio distribution is taken to be strictly independent of black hole mass. We return to this point in Section 6.1 below.

We stress that the (solid) curves in Figure 2 that show the mass functions of AGN host galaxies (with masses normalized to the Schechter ) at different AGN luminosities (computed relative to the knee of the luminosity function) are an easily testable prediction of the convolution approach, modulo the effects linked to the choice of and discussed above, which can shift the peak of .

4. Data

We next turn to demonstrate how our model relates to the observations of the mass function of star-forming galaxies and QLF. In this section, we will estimate the redshift evolution of parameters that describe these two populations.

4.1. Mass function of star-forming galaxies

As explained previously, we need to know the mass function of star-forming galaxies to deduce the SMBH mass function in star-forming galaxies, , which is an integral part of predicting the QLF in Equation (3). We also want to verify the predictions of Peng et al. (2010) and Peng et al. (2012) for the time evolution of the parameters of the galaxy mass function (e.g. equation B3 from Peng et al., 2012).

Figure 3.— Evolution of and parameters in the Schecter mass function of star-forming galaxies. Data shows minimal evolution of until at least 2.5, but significant evolution in .

In order to do this, we fit the data for star-forming galaxies from Ilbert et al. (2013) at all redshifts with a single Schechter function in which the parameters (, ) are smoothly varying functions of redshift. Fitting all the data in this fashion, instead of fitting at single redshifts, enables us to follow better the evolution of the parameters. It also reduces the sensitivity to small variations which may compromise fits at a single redshift and add to the degeneracies between parameters. This procedure of simultaneously fitting all of the data of the galaxy mass function at different redshifts is analogous to the standard procedure often used in determining the evolving QLF (e.g. Hopkins et al., 2007, Ueda et al., 2014). The functional form that we use for the redshift evolution of the galaxy mass function is


where we model the redshift dependence as


with . The slope of the low mass end, , was kept constant at the local value of as there is considerable evidence that there is little or no change in the low mass slope for the redshift range considered (Peng et al., 2014).

& -2.55 0.04
& -0.26 0.06
& -1.6 0.12
& -0.88 0.16
& 10.9 0.04
& -0.53 0.11
& 3.36 0.11
& -3.75 0.13

Table 1Evolution of star-forming galaxy mass function, data from Ilbert et al. (2013)

The evolution of parameters of mass function in Figure 3 confirms that is more or less constant up until at least redshift 2.5, i.e. it verifies the Ansatz of Peng et al. (2010) which establishes a single quenching mass scale that does not change with redshift. It is important to stress that our results will not depend critically on the exact functional evolution of above , but will use the by now well established observational fact that does not change at and that the evolution in the galaxy population since that time is associated with a “vertical” evolution in . We have also performed this analysis with the dataset from Muzzin et al. (2013) and the compilation of galaxy mass functions from Behroozi et al. (2013b) and our conclusions presented in the rest of this paper do not change.

4.2. Quasar luminosity function

Hopkins et al. (2007) combined measurements of quasar luminosity functions in different bands, fields and redshifts in order to characterise the bolometric QLF at epochs . In their work, the best fit luminosity-dependant bolometric correction and luminosity and redshift-dependent column density distributions are used in order to construct an estimate of the bolometric QLF which should be consistent with all of the various individual surveys. Although now several years old, we believe that this synthesized QLF compilation remains the most comprehensive and is used as the basis of the current paper.

We proceed with fitting their tabulated data with a double power law QLF, as given by Equation (1). We fit both at each redshift individually, and carry out a “full fit” where the parameters (i.e. , , and ) are all constrained to be smoothly varying functions of redshift, i.e. adopting a similar approach as we used earlier in our fitting of the galaxy mass function in Section 4.1. This “full” fitting of a redshift-dependent double power-law is essentially the same as the approach used in Hopkins et al. (2007), with one important difference. In order to avoid degeneracies in the fits it is necessary, both here and in Hopkins et al. (2007), to fix one or more of the parameters. Hopkins et al. (2007) chose to require a redshift-independent at a constant value. We now know that of the galaxy population changes significantly over the redshift range of interest and, as developed in the previous section, the relationship of relative to is of great interest in the context of the duty cycle. In contrast, is set, in our convolution model, by the low mass slope of the mass function of star-forming galaxies. As mentioned before there is not much evidence (see Peng et al., 2014) that this changes significantly, if at all, over the redshift range of interest, nor compelling evidence for a change in . Therefore, we choose to have a redshift-independent faint end slope of the luminosity function, and to allow to vary with redshift. It should be noted that our fitting procedure is the same as a “luminosity and density evolution” model (e.g. Aird et al., 2010) except that we are allowing the bright end slope to vary. We do this because, as we have seen, is set by the high behavior of the Eddington ratio distribution .

Figure 4.— Fits to the bolometric dataset compiled in Hopkins et al. (2007), with the broken power law fit of Equation (1). The orange curve shows our ”full fit” done with parametrization in Equation (18). Red curve shows fit in which redshift dependence of is set to be exactly the redshift dependence of .

The parameters in the luminosity function are allowed to vary as


with again . We follow Hopkins et al. (2007) in fitting and as a polynomial in (rather then ), but this choice is not of great importance.

& -4.35 0.06 & -1.79 0.04
& 0.059 0.07 & -
& -3.2 0.3 & -
& 1.73 0.2 & -
& 44.67 0.06 & 44.67 0.04
& 4.02 0.07 & 4.09 0.09
& 3.78 0.08& 3.5 0.07
& -4.68 0.1 & -3.85 0.07
& -5.7 0.2 & -5.98 0.15
& 0.5 0.03 & 0.52 0.04
& 1.64 0.05 & 1.63 0.06
& 0.54 0.05 & 0.56 0.06
& -0.12 0.01 & -0.106 0.03

Table 2Best fits to QLF data from Hopkins et al. (2007)

Additional degrees of freedom were added to the fit until we can find no appreciable quantitative difference in the quality of the fit. The values of the double power-law parameters in the smoothly varying “full-fit” are given in Table 2 and the fits are shown with orange curves in Figure 4.

The redshift dependences of and are shown in the panels of Figure 5 and we discuss these results in the next Section.

Figure 5.— Top left: Redshift evolution of QLF parameters and in our ”full fit”. The contours are showing 1- allowed regions of parameter space for fits which were made with data at each individual redshift. The filled circles are the result of a global fit, for which the resulting QLF is shown in Figure 4. Uncertainty contours for this fit are not shown here for clarity. Top right panel: Projection showing explicitly the change of normalization with redshift. Bottom left: Projection showing explicitly change of with redshift.
Figure 6.— Same as Figure 5 for QLF fit at redshifts , in which we demanded that , i.e. that redshift dependence of theses two variables is the same.
Figure 7.— Top: Comparison of the redshift evolution of the normalization of star-forming galaxy mass function and normalization of quasar luminosity function. Black circles show the values of which were determined by fitting the data from Ilbert et al. (2013) at a single redshift. Black squared show the values of which were determined by fitting the QLF data at a single redshift and are also shown in Figures 5 and 6. Bottom: Redshift evolution of the . This ratio is remarkably constant over redshift range for which data is available.

5. Results from comparing the evolution of the QLF and the galaxy mass function

5.1. The evolution of

The first result is that we notice a strong similarity between the observed redshift dependence of in Figure 5 and the observed in Figure 3. Both drop by 0.5 dex by redshift 2. Their relative evolution is explicitly compared in the bottom right panel of Figure 7, which shows that a constant ratio between and is perfectly consistent with the data. We note that the fact that the density normalization of the QLF decreases slowly with redshift has been seen in several previous analyses, for instance in the ”luminosity and density evolution” model of Aird et al. (2010), which has a 0.4 dex decrease in normalization between redshifts 0 and 2, Croom et al. (2009) shows a very similar decrease of normalization in his ”luminosity and density evolution” model, while Hasinger et al. (2005) show a decrease in normalization of 0.5 dex between redshift 0.6 and 2.4. Here we highlight the striking similarity of this behavior to the observed decline in .

Referring back to Section 3.3, a constant ratio between and implies a constant duty cycle of the black holes in the context of our convolution model. To explore this further, we now repeat the fitting procedure but set the evolution of to be exactly that of , i.e. we impose that the evolution of the QLF normalization is the same as the evolution of the normalization of star-forming galaxies in the redshift range where we have available (z3.5), while the constant multiplicative offset between these parameters remains a free parameter to be determined by the fitting procedure, i.e.


The results of this fitting are given in Table 2 and the resulting QLF is shown with red curve in Figure 4. It can be seen that the fit is extremely good, being only marginally worse then the full fit of Hopkins et al. (2007) ( ; ), both of which have the comparable number of free parameters. The main driver for the slightly worse fits in our work is the deviation of the fit from data for low luminosities at low redshifts. The data in this regime may be contaminated by contributions from the stellar populations of the hosts, as discussed in Hopkins et al. (2007).

The parameter evolution in this ”-matched” fitting are shown as linked filled circles in Figure 6. As one can see, the points are typically situated on the edges of the individual 1- contours, which is to be expected given that . We conclude that the change of normalization of the QLF is perfectly consistent with the change of normalization of the star-forming galaxy mass function over the entire redshift range for which we have the measurements of normalization of the star-forming galaxy mass function.

The fact that as far as we can tell the and track each other throughout cosmic time (at least since , and quite possibly since earlier epochs also, is an interesting result. It means that the factor which is connecting these two quantities, which in the convolution model is (slightly modified by ), remains constant.

As we have discussed above, this suggests that the general ”duty cycle”, , stays constant over cosmic time (e.g. Equations (8) and (11))(see also Conroy & White, 2013). Clearly, this would not be the case for other definitions of “duty cycle” that are based, for instance, on the fraction of black holes radiating above some luminosity or accreting above some Eddington ratio, if or evolves with time, as it does (see below), but we believe that our definition of duty cycle is the most natural one (as discussed above).

5.2. The evolution of

The other striking feature of the quasar luminosity function is the strong redshift evolution of , which increases by almost two orders of magnitude back to . At higher redshifts, the certainly flattens out and probably declines, although this cannot be established at high significance. The strong initial rise with redshift is seen independent of whether we use the full fit, or the fit constrained to have and tracking each other.

In our convolution model, the steep rise in (see Equation (13)) could have been caused in principle by one or more of (a) an evolution of , (b) an evolution of of the galaxy mass function or (c) an evolution of the mass ratio .

As discussed in Section 4.1, we know that the characteristic of the galaxy population does not change, especially in the redshift range where increase of is most prominent (), so case (b) will not apply. There is however complete degeneracy between cases (a) and (c), i.e. the distribution of specific accretion rates of AGN and the black hole to stellar mass ratio. This is clear in Equation (13) and as has been pointed out by e.g. Veale et al. (2014). This degeneracy can only be broken if we have information on the black hole masses of the AGN. We will explore this in the next section of the paper.

6. Testing the model with black hole mass data

6.1. The quasar mass-luminosity plane

Figure 8.— Mass-luminosity plane shown in 3 representative redshift bins. Our model with assumed non-evolving relation is in the top row and our model with assumed evolving relation is in the row below. Observational data from Shen et al. (2011) and Trakhtenbrot & Netzer (2012) are shown in bottom two rows. The thick black line is the Eddington limit, dashed orange lines show the calculated luminosity selection limit for lowest and highest redshift in the bin. The dotted and dashed black lines represent FWHM=1000 km/s and 1500 km/s, respectively, and are shown here to indicate which objects could be missed in observations because of FWHM limit in quasar selection. Contours are set at 10, 20 etc. values of estimed probability distribution of objects. Outermost objects are represented as individual dots.

In this part of the paper we will compare predictions of how the black hole mass-luminosity plane of broad-line AGNs should be populated in our model, and compare these with SDSS data. After creating a mock sample of star-forming galaxies in an SDSS volume, we will populate them with black holes assuming different redshift-dependent scaling relations and an assumed scatter. We will then assign Eddington ratios from the evolving distribution and apply an obscuration prescription from Hopkins et al. (2007). These two functions are chosen so that the QLF is reproduced as in the previous section: in other words the redshift evolution in and the redshift evolution in must combine to produce the observed evolution in the QLF .

For the observational distribution, we use data from two observational studies (Shen et al., 2011; Trakhtenbrot & Netzer, 2012) to show that our results do not depend critically on the data choice and to give the reader a graphical impression of the uncertainties involved in this kind of measurement. By comparing our mock data with the observed distribution we will be able to see which combination of redshift-dependent and best reproduces the observational data. We take into account the obscuration factor and the differences in bolometric correction between different works and apply the same bolometric correction to the data and model (namely one used in Shen et al., 2011) to make them directly comparable.

By using this mock sample approach we can fully account for biasses introduced by the luminosity-selection of the quasar samples. We recreate data only up to , as black hole mass estimates for higher redshifts are based on the broad C IV emission line, which was shown to be far less reliable for these purposes (Shen et al., 2008; Trakhtenbrot & Netzer, 2012 and references therein).

We first show in the topmost panels of the Figure 8 the modelled distribution of quasars in the black hole mass-luminosity plane in 3 representative redshift bins if we assume a redshift independent with the standard value of . We introduce a log-normal scatter in this relationship of 0.5 dex to account for scatter in relationship. The solid diagonal line indicates the Eddington limit ().

The data from (Shen et al., 2011; Trakhtenbrot & Netzer, 2012) are plotted in the two bottom rows of panels in Figure 8. In these panels, the diagonal dashed and dotted lines indicate the locus of black hole masses for two constant FWHM (of 1000 km/s and 1500 km/s respectively) of the emission lines (H and MgII) that were used to infer the black hole masses. Systems of lower FWHM will not appear in the samples, e.g. the limit used in Shen et al. (2011) is set at 1200 km/s.

We see that, although there is good agreement at low redshifts, a non-evolving relation produces far too many objects with masses that are too small or, equivalently, which have very high Eddington ratios, with around 50 being super-Eddington in the final redshift bin, while only around 2 of objects are super-Eddington in the data. This is a simple consequence of the fact that is much higher then locally and is a reflection of the high implied for an unevolving relation shown in Figure 17.

It should be noted that the problem gets worse if we reduce the scatter in the relation as we then have even fewer high mass black holes. We note that the comparison could not be expected to be perfect because our data is constrained to reproduce the QLF from Hopkins et al. (2007); although SDSS data and the optical QLF is the main contributor to the Hopkins QLF in this luminosity range, there are small contributions from other surveys as well as slightly different bolometric and obscuration corrections which will induce small differences. Nevertheless, the disagreement with a non-evolving ratio is too large to be due to this.

A much better agreement is obtained (second row of Figure 8) if we adopt an evolving relation. We adopt as a heuristic example the form with . The agreement with the observed distribution is considerably better and there are now far fewer objects crossing the Eddington limit ( 3 ) at high redshifts. This means that the observed increase in back to would be due to an equal split between a change in and a change in characteristic Eddington ratio , remembering that these two changes are degenerate in our convolution model without the data of Figure 8.

Finally we note that, quite independently of any assumptions about , our convolution model naturally recreates the apparent “sub-Eddington boundary” that has been emphasized by Steinhardt & Elvis (2010b), by which we mean the flat upper envelope . This refers to the fact that at all redshifts there seems to be a lack of objects at high masses close to the Eddington limit, which can also be seen in the Figure 16 of Trakhtenbrot & Netzer (2012). This can be observed on the Figure 8 where the upper contours of the red regions have slopes that are noticeably shallower than the 45 degree slope that corresponds to a constant Eddington ratio, thereby giving the impression of an absence of high luminosity high sources.

This behaviour is quite counter-intuitive, and was interpreted by Steinhardt & Elvis (2010b) as being caused by some new physical effect that somehow limits accretion onto more luminous quasars. Various authors have proposed alterations to the measurement methods in mass-luminosity plane which could reduce or eliminate this effect (e.g. Rafiee & Hall, 2011a; Rafiee & Hall, 2011b). We show instead that this behaviour is actually expected from the convolution model presented in this paper!

We commented earlier in Section 3.4 that while the typical black hole mass increases with luminosity in our convolution model, it does so sub-linearly, so that the “typical” Eddington ratio must also increase with luminosity. This can be seen also in these plots: the ridge line that is defined by the peak in the mass distribution at a given is indeed steeper than the 45 degree line of constant Eddington ratio, which is why the shallower slope of the boundary defined by the contours is so counter-intuitive. However, because the number of more massive black holes falls very steeply with mass, because of the decline of the galactic mass function above , the contours of constant surface density, given by the contours in the plots in Figure 8 and in the Steinhardt & Elvis (2010b) analysis, are actually shallower than the 45 degree line.

This is more clearly seen in Figure 9 were we show the appearance of our full sample, before the application of the SDSS luminosity cut. At lower luminosities the effect gets more and more pronounced and we expect that deeper surveys of a given area would find more super-Eddington objects at small masses.

We emphasize that the reproduction of the Steinhardt & Elvis (2010b) effect in Figure 8 is achieved within our convolution model with an Eddington ratio distribution that is completely independent of black hole mass. The effect noted by Steinhardt & Elvis (2010b) would not be seen (in our model) if the distribution of points in the plane was normalized to the total number of black holes (in star-forming galaxies) at a given mass, which is, of course, unfortunately not observable. Our explanation for the sub-Eddington boundary in Steinhardt & Elvis (2010b) is therefore a kind of “plotting bias” arising from what is plotted, rather than a “selection effect” per se, coming through the construction of the sample via, e.g., emission line widths.

Figure 9.— Full sample from our model, before applying the luminosity cuts to produce the second row in Figure 8. For clarity, only one out of every 25 points is shown.

6.2. Establishing the correlation in quenched systems

In this section we show how it is possible to reproduce the observed tight correlation between in the local Universe even if the black hole to stellar mass relation in the star-forming galaxies is strongly evolving (see also Croton, 2006; Hopkins et al., 2006a). In what follows, we loosely equate with the stellar mass of a star-forming galaxy when star-formation quenches, which is also the point at which we assume the black hole ceases growing. Readers uncomfortable with making this association may simply substitute in what follows with .

Figure 10.— Top left: Resulting relation in quenched galaxy systems today, if we only take into account systems that have quenched after . Orange line shows best fit quoted in Kormendy & Ho (2013). Blue line shows best fit to the data generated from 40 realizations of the model. Top right: Resulting relation in quenched systems today, for full range of quenching history (). Orange and blue line are generated in same way as in previous panel. Merging, uncertain evolution in relation and errors in estimated rate of quenching at high redshift could bias this result, as discussed in text. Bottom left: Measured by Kormendy & Ho (2013). Orange line shows again best fit quoted in Kormendy & Ho (2013). Bottom right: Measured by Häring & Rix (2004).

For this analysis, we use results from the simple galaxy evolution model of Birrer et al. (2014) to construct a set of evolving star-forming galaxies and their quenched descendants. The details of the Birrer et al. (2014) model are not important for the present purpose since it reproduces well the overall evolution of the galaxy population, which is all that we require here.

For each quenched galaxy seen in the model at the present epoch, we know the mass and redshift at which it quenched and can therefore compute the black hole mass from the adopted redshift-dependent relation for (star-forming) galaxies, adding also the adopted scatter. We assume that there is no stellar mass growth after quenching and that there is no central SMBH mass growth after quenching

After this procedure we are left with a mock sample of quenched galaxies, with their black hole masses known, after which we can account for sample selection effects. This is not trivial, as measurements of masses of black holes are from heterogeneous sources and no single luminosity or other cut is possible.

We therefore decided to create an empirical selection in galaxy mass so that the distribution of galaxy masses in the mock sample broadly matches that of the passive early type galaxies which have had their SMBH mass measured. We show results for two situations, where we first only consider galaxies that quenched after and then include all galaxies that have quenched since . We differentiate between these two cases since we expect merging, which is not explicitly accounted for here, to have a larger impact for galaxies that quench earlier, and because the assumed mass ratio scaling, which we think is reasonable approximation (see above) at redshift , may break down at higher redshifts. When we ignore quenching at , we are losing only about 10% of today’s quenched population, but the model is on a much firmer footing.

We have fitted the simulated data, derived from 40 random realizations of the model, with a relation of the form


by regressing the black hole mass on the stellar mass, and compute the scatter of the simulated galaxies around this relation. We derive values of and a scatter of 0.53 dex for the case of quenching only from z=2, and values of with scatter of 0.6 dex for the case of quenching from z=5. These values should be compared with the observed values of and a intrinsic scatter of 0.29 dex derived in Kormendy & Ho (2013). While the predicted scatter appears to be slightly larger than observed, subsequent merging of galaxies, that has not been modelled here, will act to reduce the scatter (Hirschmann et al., 2010, Jahnke & Macciò, 2011).

The origin of relation is the underlying relation in star-forming galaxies assumed in the model, with added scatter due to the fact that galaxies of a given stellar mass today have quenched at various redshifts and, for a given stellar mass, the spread in black hole masses is amplified by spread in the quenching redshifts. The mean of the relation is positioned roughly at the value of at the mean redshift of quenching at that mass, which is over a wide range of galaxy masses.

It is interesting to note that objects that have quenched more recently would be expected on average to have a lower value, because of the evolution in . It is possible that these recently quenched galaxies could be associated with pseudo-bulges instead of bulges. Kormendy & Ho (2013) has indeed argued that pseudo-bulges do indeed have a lower ratio.

We also note in passing that if there is a trend for the typical quenching redshift to increase with increasing stellar mass (e.g. because of interplay of mass and environment quenching, Peng et al., 2010) then this would introduce a tilt in the local relation () even though the input in star-forming galaxies was perfectly linear.

The link with quenching redshift then prompts another interesting point. There is very good evidence that the sizes of galaxies, at a given mass, are smaller at high redshift (Daddi et al., 2005; Newman et al., 2012; Carollo et al., 2013; Shankar et al., 2013a). We would therefore expect, through virial arguments, that the velocity dispersions, at a given mass, would also be higher. This has been directly observed in a few cases (e.g. van Dokkum et al., 2014). For analytic simplicity, we assume that the scale radius of galaxies scales as so that the usual Faber-Jackson type scaling relation would be expected to evolve as


If scales as , as we have been exploring in this part of the paper, then this naturally produces an relation that is independent of redshift.

This has two implications: first, we would not expect to see any significant evolution in the observed relation (see e.g. recently Shen et al., 2015). Second, we would expect the present-day relation for passive galaxies to be tighter than the relation, because of the aforementioned broadening of the latter from the range of .

We stress that this tighter relation would be present even if the velocity dispersion is not playing a direct role in the growth of black holes. Of course, it is also possible that the evolution in the is in fact caused by the constancy of an underlying causal connection. On the other hand, there could well be other causes for to evolve as , in which case the tight relation would simply be a coincidence.

6.3. in AGNs in the local Universe

Finally, we turn to estimates of the relations for AGNs in the local Universe, where it is possible to estimate independently the mass of the galaxy that is hosting an optically active AGN. Matsuoka et al. (2014) made a careful decomposition into nuclear and host contributions of the images of quasars in the SDSS Stripe 82. They found that the quasars are predominately hosted in massive star-forming galaxies, with relatively large ratios of around . Matsuoka et al. (2014) were aware of the possibility that selection effects could bias this value upwards but could not estimate their magnitude. We can now use our model to examine the expected size of this bias.

To do this we simulate sources within the same sky area as Stripe 82 and recreate objects above the AGN luminosity cut that would have been selected for quasar spectroscopy in the SDSS sample, with scatter of dex. The results are shown in Figure 11, represented in the same way as in the original paper. We see that the observed quasars will have ratios that are indeed much higher then the mean value in the underlying sample, which is indicated by the shaded region in the figure.

Figure 11.— Top panel: Simulated observation of relation in Stripe 82, mimicking analysis of Matsuoka et al. (2014). Blue line shows mean relation in full sample (before luminosity cut was applied), while the shaded area shows spread around mean relation in our model. The black line shows mean and spread of distribution of points in 0.1 redshift slices. Red line is set to , which is approximately local relation from Häring & Rix (2004). Bottom panel: original data from Matsuoka et al. (2014).

Finally, we notice that even though we inserted redshift evolution in our model, this evolution would be quite hard to detect in the sample of Matsuoka et al. (2014), owing to the large spread of points and the selection biasses connected with such a study. Nevertheless, we do still see a slight change in the mean values of the simulated sample that is not seen in the actual data. We discuss this further in Section 7.2 below.

7. Discussion

In the preceding Sections of this paper we have developed a simple generic model for obtaining the evolving AGN luminosity function from a convolution of the evolving galaxy mass function. This generic model makes testable predictions for quantities such as the shape of the mass distribution of host galaxies as a function of AGN luminosity and allows us to derive quantitative connections between the parameters describing the galaxy mass function and the AGN QLF. On the basis of these, we can derive powerful statements about the duty cycle of AGN. We then showed, in the framework of this general model, that a redshift-dependent and Eddington ratio distribution, , successfully reproduces the observed quasar luminosity function (by construction) and also reproduces observations of the distribution of quasars in the () plane, the black hole to bulge mass relation of quenched galaxies and measurements of for low redshift AGN. In this Discussion section, we develop further some of the astrophysical implications of the preceding results.

7.1. Comparing redshift evolution in active and quenched systems

It is important to appreciate that a quite rapid evolution in in active (star-forming) AGN systems does not imply an equally rapid evolution of in the quenched systems. This is because, when we are observing quenched systems, we are effectively observing the integrated history of previous activity. To illustrate this, we perform a simple heuristic exercise in which we determine in quenched systems at given epoch, assuming as above that quenching started at redshift 2. We make the same assumptions as before, namely that there is no stellar or black hole mass growth after quenching.

Our results are shown in Figure 12. Even though star-forming galaxies have changed their black-hole to stellar mass ratio by almost a factor of 10 from redshift to today, the evolution in this ratio for quenched systems is much milder, more like a factor of 3 or even less (0.4 dex), simply because the quenched population includes galaxies that have quenched much earlier. This means that the redshift evolution in for quenched systems will always be much milder than the evolution in this quantity in active systems.

The distinction between active and passive populations becomes even more important if observational studies compare actively accreting systems at high redshift with data on the quiescent population at low redshift. This is unfortunately quite common practice (e.g. Jahnke et al., 2009; Decarli et al., 2010; Merloni et al., 2010; Bennert et al., 2011; Schramm & Silverman, 2013; Schulze et al., 2015), because of the current practicalities of observations. We stress that this approach automatically produces weaker evolution since neither nor will have changed once a given object becomes inactive (i.e. passive/quenched). It is therefore clear that the mean relation in the quenched systems will reflect the mass ratio in the star-forming galaxies at the much earlier epochs when those galaxies actually quenched (as already discussed above). For the local population, this is typically around z 1-1.5.

Clearly, if we compare the of star forming systems at z 1-1.5 to the ratio of passive galaxies seen today that quenched at z 1-1.5, then we would expect to see no change in , even if this ratio had changed a lot within the actively accreting population! Of course, if the high redshift active systems are luminosity-selected then their ratio will likely have been biassed to higher values than the underlying population through the “Lauer effect” discussed in Section 3.4, which goes in the opposite direction (as shown by the parallel black lines in Figure 12). This figure shows that if Nature has an underlying scaling as for active systems, then we would expect to see if we compare comparisons of active systems at with quenched systems at , once we correct for these observational biases. This is quite similar to the evolution seen in at least some of the observational studies cited above (e.g. Jahnke et al., 2009; Decarli et al., 2010; Merloni et al., 2010; Schramm & Silverman, 2013; Schulze et al., 2015).

We have already commented in Section 6.2 that we would also expect to see little or no evolution in for any population of galaxies, however selected. This is due to the higher associated with a given stellar mass at high redshift which cancels out the dependence in .

Figure 12.— Evolution of the mean relation in star-forming (heavy blue line) and quenched systems (heavy red line). The dashed blue line shows the redshift range where this mean relation may not accurately represent star-forming galaxies, as discussed in Section 7.2. The evolution in for the quenched (inactive nuclei) galaxies is always shallower than for the population of actively accreting systems, potentially leading to an underestimate of the evolution if the populations are mixed. On the other hand, the “Lauer effect” will bias the upwards in luminosity selected active samples, which goes in the opposite direction. This is indicated by the thin black lines which show the observed for active systems for luminosity-selection (labelled relative to ) for our standard assumption of = 0.5 dex.

7.2. Mass growth of star-forming galaxies at low redshift

In this paper we have made the conventional assumption that the black hole to stellar mass relation increases as some power of (1+z), i.e. . However, this is an arbitrary redshift dependence, and there are reasons why it cannot hold (with ) at low redshifts (see also Hopkins et al., 2006b). Assuming that a star-forming galaxy is on the Main Sequence we can track the increase of its stellar mass using the Main Sequence sSFR(z),


where is the “reduced specific star-formation rate” (see Lilly et al., 2013), , where R is the fraction of stellar mass that is returned during star formation , which is therefore the inverse mass doubling timescale of the stellar population. Using Gyr from Lilly et al. (2013) and references therein, this can be integrated to give


This modest increase in stellar mass for a star-forming galaxy sets a maximal change in ratio even in the most extreme case that does not increase at all. We show this maximal evolution in Figure 13 and compare it with the evolution used elsewhere in the paper. It can be seen that the maximal evolution is actually slower then at redshifts . Galaxies are growing so slowly at low redshifts that it is not possible to create a strong evolution in ratio in a given galaxy, even if their black holes are not growing at all. This may well be why there is no evidence in Figure 11 for the ”expected” increase in the observed . The maximal change in over this redshift range would have been gentle enough that it would be very hard to observe. For this reason, we emphasize that our statements elsewhere in the paper concerning the evolution in should best be interpreted as implying a change of a factor of about ten to , rather than a precise particular dependence on redshift.

Figure 13.— Comparison of the maximal evolution of the relation in a star-forming galaxy that stays on the Main Sequence, with our assumed . At low redshift galaxies are growing so slowly in the stellar mass that strong evolution in the ratio is not possible, even if the black hole does not grow at all.

7.3. Downsizing

A number of authors (e.g. Hasinger et al., 2005; Barger et al., 2005; Labita et al., 2009; Li et al., 2011) have noted or discussed a “downsizing” of the quasar population. Although different authors often use this term to mean different things, it is most often used to denote the observational fact that lower-luminosity AGNs peak in comoving density at lower redshifts then higher-luminosity AGNs.

It is worth stressing that this may not have much physical significance. We have shown that it is possible to reproduce the strong observed redshift evolution in the QLF with a model based on the observed mass function of star-forming galaxies coupled with a mass-independent (but redshift-dependent) Eddington ratio distribution . This is shown more explicitly in Figure 14 where we show the comoving number density of quasars of different luminosity in the QLF which is reproduced by our model. This emphasizes that the apparent down-sizing signature in the AGN population can appear even though the distribution of Eddington ratios (and thus of specific black hole growth rates) is strictly independent of black hole mass at all redshifts in our model.

It is clear that the apparent “downsizing” in our model arises as a natural consequence of two competing effects which are independent of mass. The first is the redshift evolution of the which shifts the luminosity function uniformly in luminosity, but which therefore produces a differential change in number density with luminosity. This shift is produced by the degenerate combination of evolution of the mass ratio and characteristic Eddington ratio, . The second is the redshift evolution of that changes the number densities uniformly with luminosity. We have argued in this work that the evolution of is a direct consequence of the evolution of , coupled with a constant duty cycle. The combination of these two produces variation in the number density (at fixed luminosity) that changes with luminosity, producing a variation in the redshift of the peak in the number density as well as different rates of evolution at different luminosities.

Figure 14.— Redshift evolution of comoving density of quasars for several different luminosities. The density of lower luminosity AGN peaks at lower redshift then for high luminosity AGNs. This behaviour is naturally produced in our model even though the distribution of Eddington ratios is completly independent of black hole mass.

7.4. Coherent evolution of and

One of the most striking results of this paper is that the observed evolution of of the QLF tracks the observed evolution of of the star-forming galaxy mass function. This in turn implies that the knee of the Eddington ratio distribution has a more or less constant value, even though the change of the characteristic ”knee” is dramatic. We have argued that is a good measure of a generalized “duty cycle” of quasars. Indeed, if the luminosity of individual quasars decays with a timescale then (Equation 10) will be direct measure of the birth-rate of quasars in star-forming galaxies.

By applying simple continuity equations to the observed evolution of the star-forming galaxy mass function, Peng et al. (2010) derived an expression for the rate of the mass-quenching process, , which may be written


where sSFR(z) is the redshift dependent specific star-formation rate of the star-forming Main Sequence. It has been well established that sSFR was much higher in the past and a useful representation is (Lilly et al., 2013 and references therein)


AGN activity in ”quasar-mode” is one of the many processes that have been proposed (e.g. Granato et al., 2004; Hopkins et al., 2008a; King, 2010) to drive the mass-quenching of galaxies. If a single quasar event is responsible for quenching, then this would require that (from Equation 10) to be equal to . These can only be equated for our inferred constant for a particular redshift and mass dependence of the decay time .


We would require quasars to fade faster at high redshift and at high galaxy masses. We could well imagine ways in which this would occur, e.g. because of the shorter dynamical timescales of galaxies at high redshift.

However, the above picture is probably over-simplistic. We could well expect the physics of quenching to be more complex. It is quite plausible that the energy source for quenching is AGN activity but that the effectiveness of this energy injection depends on the stellar or halo mass of the system (the Peng et al., 2010 quenching law can be written in terms of a redshift-independent survival probability that depends only on mass). This would then break the simple link between and .

8. Summary and conclusions

We have presented a simple, phenomenological model that aims to link the evolving galaxy population with the evolving AGN population. We use our observational knowledge of the evolving galaxy mass function and of the evolving quasar luminosity function (QLF) to connect these two populations and to create a global model to interpret the AGN population, including biases associated with sample selection.

Our model is based on three observationally motivated Ansätze, namely that

  • radiatively efficient AGNs are found in star-forming galaxies,

  • the probability distribution of the Eddington ratio does not depend on the black hole mass of the system,

  • the mass of the central black hole is linked to the stellar mass.

The QLF is then a straightforward convolution of the black hole mass function with the distribution of Eddington ratios , while the former is itself a convolution of the galaxy mass-function with the relation. These heuristic assumption ensure that our model is simple enough to be analytically tractable, while still capturing the main characteristics of the galaxy and AGN population.

The main conclusions of this model analysis can be summarized as follows:

  1. The “broken” or “double” power law form of the quasar luminosity function is a consequence of the underlying Schechter mass function and a “double” power law, mass independent, Eddington ratio distribution. We show how the parameters of the QLF are straightforwardly connected with the input functions. Most importantly, the knee of the QLF, , is proportional to the product of the of the galaxy mass function, the ratio and the position of the break in the Eddington ratio distribution while the normalization of the QLF is proportional to the product of the normalization of the star-forming galaxy mass function and the normalisation of the Eddington ratio distribution, which can be loosely interpreted as a ”duty cycle”.

    Our simple convolution model makes clear and testable predictions for the distribution of host galaxy masses (relative to the star-forming galaxy Schechter ) for different AGN luminosities (relative to ). At high luminosities (above the AGN ) this is a Schechter function with the star-forming M* but a faint end slope given by .

  2. There is a remarkable consistency in the redshift evolution of normalization of SF mass function and the normalization of QLF. These two characteristic densities track each other closely out to redshifts of , and possibly to even higher redshifts. This implies that the generalised “duty cycle” of AGN is surprisingly constant with redshift.

  3. In contrast, the QLF evolves strongly with redshift, with evolution being at least up to . Given that there is strong evidence for no change in the galaxy over this redshift range, this evolution in is driven by an evolution in the characteristic “knee” in the Eddington ratio distribution or in the mass scaling between stellar mass and black hole mass, , or some combination of the two. The QLF evolution is degenerate in changes of these two quantities.

We then explore this degeneracy by comparing predictions of our model, incorporating the relevant selection cuts, for the distribution of systems in the SDSS AGN mass-luminosity plane(s) using black hole mass estimates for individual AGN. We find a good match with an evolving in star-forming systems. This implies that the observed increase in back to would be due to an equal split between a change in and a change in characteristic Eddington ratio .

We show that this is compatible with the observed relations in both quenched and in star-forming galaxies in the local Universe, both in terms of the mean relations and the scatter.

We also make the important point that much weaker evolution, more like , would be deduced by observers, after correcting for the “Lauer bias”, if (as they usually of necessity do) they compare black hole masses in active (star-forming) systems at high redshift () with those in quiescent systems at low redshift. Similarly, much weaker evolution would also be seen if black hole masses were to be compared solely within the passive population at different redshifts.

The inferred evolution in star-forming systems is likely simplistic, and unlikely to hold all the way down to the lowest redshifts simply because there is not enough star-formation to change fast enough, even if black holes are not growing at all. Our use of a evolutionary model should be interpreted more loosely as a way of getting a factor of 2.5 change between in active systems at when compared with the same relation in passive systems today, rather than as a precise relation that absolutely holds at all redshifts.

We however make the interesting point that an evolution of the form would, when coupled with the observed changes in the sizes of galaxies, which typically scale as roughly , have the feature of producing a relation that would be more or less independent of redshift. This could then explain why the present-day relation would have lower scatter than the relation, even if plays no direct role in black hole growth. Alternatively, a fundamental redshift-independent relation could be the physical origin of an apparent evolution.

We stress that the most basic features of the model do not depend on the redshift dependence of , which is driven largely by the possibly uncertain observational estimates of black hole masses in high redshift AGN.

Not least, quite independent of the form of any evolution in , the generic model naturally reproduces the counter-intuitive “sub-Eddington boundary” in the () plane that has been noted by Steinhardt & Elvis (2010b) without the need to invoke any new physical effects. The generic model also produces the apparent “down-sizing” of the AGN population (e.g. Hasinger et al., 2005). We stress that both of these apparently mass-dependent effects are achieved in the model with an Eddington ratio distribution that is completely independent of black hole mass at all redshifts.

Acknowledgements We thank Simon Birrer for providing us with the full output of the model presented in Birrer et al. (2014). We also thank Peter Behroozi for providing us with the galaxy mass functions that were used as input in Behroozi et al. (2013b). Sebastian Seehars has helped in creating fitting procedures. In developing this work, we have benefited from stimulating discussions with Charles Steinhardt, Andreas Schulze, Kevin Schawinski and David Rosario. We are also very grateful for the very perceptive comments provided by the anonymous referee. This work has been supported by the Swiss National Science Foundation.


  • Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  • Aird et al. (2010) Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS, 401, 2531
  • Aird et al. (2012) Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746, 90
  • Aird et al. (2013) —. 2013, ApJ, 775, 41
  • Alexander & Natarajan (2014) Alexander, T., & Natarajan, P. 2014, ArXiv e-prints, arXiv:1408.1718
  • Aller & Richstone (2007) Aller, M. C., & Richstone, D. O. 2007, ApJ, 665, 120
  • Assef et al. (2011) Assef, R. J., Kochanek, C. S., Ashby, M. L. N., et al. 2011, ApJ, 728, 56
  • Azadi et al. (2014) Azadi, M., Aird, J., Coil, A., et al. 2014, ArXiv e-prints, arXiv:1407.1975
  • Baldry et al. (2008) Baldry, I. K., Glazebrook, K., & Driver, S. P. 2008, MNRAS, 388, 945
  • Baldry et al. (2012) Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621
  • Barger et al. (2005) Barger, A. J., Cowie, L. L., Mushotzky, R. F., et al. 2005, AJ, 129, 578
  • Behroozi et al. (2013a) Behroozi, P. S., Marchesini, D., Wechsler, R. H., et al. 2013a, ApJ, 777, L10
  • Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013b, ApJ, 770, 57
  • Bell et al. (2003) Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289
  • Bell et al. (2007) Bell, E. F., Zheng, X. Z., Papovich, C., et al. 2007, ApJ, 663, 834
  • Bennert et al. (2011) Bennert, V. N., Auger, M. W., Treu, T., Woo, J.-H., & Malkan, M. A. 2011, ApJ, 742, 107
  • Birrer et al. (2014) Birrer, S., Lilly, S., Amara, A., Paranjape, A., & Refregier, A. 2014, ApJ, 793, 12
  • Bongiorno et al. (2012) Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, 3103
  • Boyle & Terlevich (1998) Boyle, B. J., & Terlevich, R. J. 1998, MNRAS, 293, L49
  • Brown et al. (2006) Brown, M. J. I., Brand, K., Dey, A., et al. 2006, ApJ, 638, 88
  • Brusa et al. (2010) Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348
  • Cao (2010) Cao, X. 2010, ApJ, 725, 388
  • Caputi et al. (2011) Caputi, K. I., Cirasuolo, M., Dunlop, J. S., et al. 2011, MNRAS, 413, 162
  • Carollo et al. (2013) Carollo, C. M., Bschorr, T. J., Renzini, A., et al. 2013, ApJ, 773, 112
  • Cimatti et al. (2013) Cimatti, A., Brusa, M., Talia, M., et al. 2013, ApJ, 779, L13
  • Cisternas et al. (2011) Cisternas, M., Jahnke, K., Bongiorno, A., et al. 2011, ApJ, 741, L11
  • Civano et al. (2011) Civano, F., Brusa, M., Comastri, A., et al. 2011, ApJ, 741, 91
  • Conroy & Wechsler (2009) Conroy, C., & Wechsler, R. H. 2009, ApJ, 696, 620
  • Conroy & White (2013) Conroy, C., & White, M. 2013, ApJ, 762, 70
  • Croom et al. (2009) Croom, S. M., Richards, G. T., Shanks, T., et al. 2009, MNRAS, 399, 1755
  • Croton (2006) Croton, D. J. 2006, MNRAS, 369, 1808
  • Daddi et al. (2005) Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680
  • Daddi et al. (2007) Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156
  • Decarli et al. (2010) Decarli, R., Falomo, R., Treves, A., et al. 2010, MNRAS, 402, 2453
  • Eddington (1922) Eddington, A. S. 1922, MNRAS, 83, 98
  • Elbaz et al. (2007) Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33
  • Ferrarese & Ford (2005) Ferrarese, L., & Ford, H. 2005, Space Sci. Rev., 116, 523
  • Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9
  • Fine et al. (2008) Fine, S., Croom, S. M., Hopkins, P. F., et al. 2008, MNRAS, 390, 1413
  • Franceschini et al. (1999) Franceschini, A., Hasinger, G., Miyaji, T., & Malquori, D. 1999, MNRAS, 310, L5
  • Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13
  • González et al. (2011) González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, 735, L34
  • Goulding et al. (2014) Goulding, A. D., Forman, W. R., Hickox, R. C., et al. 2014, ApJ, 783, 40
  • Graham et al. (2011) Graham, A. W., Onken, C. A., Athanassoula, E., & Combes, F. 2011, MNRAS, 412, 2211
  • Graham & Scott (2013) Graham, A. W., & Scott, N. 2013, ApJ, 764, 151
  • Granato et al. (2004) Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580
  • Greene & Ho (2006) Greene, J. E., & Ho, L. C. 2006, ApJ, 641, L21
  • Greene & Ho (2007) —. 2007, ApJ, 670, 92
  • Greene et al. (2008) Greene, J. E., Ho, L. C., & Barth, A. J. 2008, ApJ, 688, 159
  • Gültekin et al. (2009) Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 695, 1577
  • Häring & Rix (2004) Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89
  • Hasinger et al. (2005) Hasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441, 417
  • Heckman et al. (2004) Heckman, T. M., Kauffmann, G., Brinchmann, J., et al. 2004, ApJ, 613, 109
  • Henriques et al. (2014) Henriques, B., White, S., Thomas, P., et al. 2014, ArXiv e-prints, arXiv:1410.0365
  • Hickox et al. (2014) Hickox, R. C., Mullaney, J. R., Alexander, D. M., et al. 2014, ApJ, 782, 9
  • Hickox et al. (2009) Hickox, R. C., Jones, C., Forman, W. R., et al. 2009, ApJ, 696, 891
  • Hirschmann et al. (2010) Hirschmann, M., Khochfar, S., Burkert, A., et al. 2010, MNRAS, 407, 1016
  • Hopkins et al. (2008a) Hopkins, P. F., Cox, T. J., Kereš, D., & Hernquist, L. 2008a, ApJS, 175, 390
  • Hopkins & Hernquist (2009) Hopkins, P. F., & Hernquist, L. 2009, ApJ, 698, 1550
  • Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2005, ApJ, 630, 716
  • Hopkins et al. (2006a) —. 2006a, ApJS, 163, 1
  • Hopkins et al. (2008b) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008b, ApJS, 175, 356
  • Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS