AGN evolution from a galaxy evolution viewpoint

Abstract

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 L∗ can be produced by either an evolution in the distribution of Eddington ratios, or in the mbh/m∗ mass ratio, or both. To try to break this degeneracy we look at the distribution of AGN in the SDSS (mbh,L) plane, showing that an evolving ratio mbh/m∗∝(1+z)2 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 mbh−σ relation and could explain why the local mbh−σ relation is tighter than mbh−m∗ even if σ is not directly linked to black hole growth.
Irrespective of the evolution of mbh/m∗, 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.

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 z∼2 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 M∗ of the mass-function ϕ(m) of the star-forming population has been more for less constant back to at least z∼2.5 (and likely to z∼4) 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.

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 L∗ at redshifts below z∼2. 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 ϕ(m) 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 M∗. 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

L=1038.1⋅λ⋅mbh

(2)

where L is given in erg s−1, λ is the Eddington ratio and mbh 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 mbh 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 mbh 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 mbh−σ 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 (mbh−σ) and with the stellar mass of the bulge in the galaxy (mbh−mbulge) (Magorrian et al., 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Marconi & Hunt, 2003). Canonical values of the ratio between mbh and mbulge are between 10−3 and 10−2.7(Häring & Rix, 2004). The recent extensive review by Kormendy & Ho (2013) has argued for larger values, up to 10−2.3, 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 mbh−n where n is Sersic index describing the light profile of the galaxy, mbh−mhalo where mhalo is mass of the dark matter halo and mbh−m∗ where m∗ 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 mbh−m∗ is more fundamental then mbh−mbulge 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 mbh−m∗ 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.

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 ϕ(m∗) 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-LAGN 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 (mbh,L) 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 ϕ(L) can be simply constructed from the galaxy ϕ(m∗) 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 ϕSF(m∗,z) and ϕ(L,z) 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 mbh/m∗ ratio and show that one particular choice can explain the observed distribution of luminous quasars in the (mbh,L) 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 mbh/m∗ implementation, differences between mbh/m∗ and mbh−σ 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 ΩM=0.3, ΩΛ=0.7 and H0= 70 km s−1 Mpc−1. Luminosities will be given in units of erg s−1 and will refer to bolometric luminosities, unless specified otherwise. We use the term ”dex” to denote the antilogarithm, i.e. n dex = 10n. We also define all distribution functions, i.e. the star-forming galaxy mass function ϕSF(m∗), the associated star-forming galaxy black hole mass function ϕBH(mbh), the AGN luminosity function ϕ(L) 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 ϕSF(m∗), ϕBH(mbh) and ϕ(L) are Mpc−3 dex−1.

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 (mbh∝mβ∗), with some scatter, and we will for simplicity set β∼1.

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 β∼1, 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.

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 mbh/m∗, 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 η=0 for initial analytic simplicity, before reintroducing it with η∼0.5 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 ϕBH(mbh,z) with a probability distribution of the Eddington ratio λ,

ϕ(L,z)=∫ϕBH(mbh,z)ξ(λ,z)dlogλ,

(3)

where ϕ(L,z) is the resulting QLF, and ξ(λ,z) 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 ϕBH(mbh) (itself derived from the stellar mass function of star-forming galaxies ϕSF(m∗)) 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 mbh/m∗ scaling to get the mass function of SMBH in star-forming galaxies, with scatter as required (2nd panel). We convolve it with the Eddington ratio distribution (in the 3rd 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 γ2) 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 δ2 we can thus equate

γ2=δ2,

(4)

with γ2 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 γ1 will be given by the steeper of the low end slope of ϕBH(mbh) and the low end slope of ξ(λ), which we denote by δ1 (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

γ1=max(−αBH,δ1).

(5)

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 mbh vs. m∗ 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 αBH=αSF/β.

We note at this point that there is good observational evidence that the faint end QLF slope, γ1, is similar to the observed low mass slope of the mass function of star-forming galaxies, αSF, allowing for the reversal of sign in our definition.
The QLF faint end slope, γ1 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 αSF 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 αSF∼−0.4.

The data are therefore consistent with the idea that γ1∼−αSF. This suggests that γ1 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 β∼1. We will henceforth assume that this is the case, i.e. that αBH=αSF≈−γ1. This then means that the low end slope δ1 of ξ(λ) can take any value that is shallower than this, i.e. δ1<−αSF∼0.4.
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, fd, of objects are active above a certain threshold value λ∗ with an Eddington ratio distribution that is a single power law with slope δ2, while (1-fd) of objects are completely inactive. Effectively this sets the low end slope δ1∼−∞.

The exact functional form of ξ(λ) can then be written as

ξ(λ)=dNNdlogλ=ξ∗λ(λλ∗)−δ2,λ>λ∗,

(6)

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

In an Appendix we show that the knee of the QLF, L∗, 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.

L∗=1038.1⋅λ∗⋅M∗BH⋅ΔL(γ2),

(7)

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

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

ϕ∗QLF=ϕ∗SF⋅ξ∗λ⋅Δϕ(γ2).

(8)

where the Δϕ denotes once again a small multiplicative factor, weakly dependent on γ2 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 L∗ is connected with the normalization of the contributing functions at M∗BH and λ∗. Simply put, if there are more objects at mass M∗ and more of them are active at λ∗, then we expect also more objects with corresponding luminosity L∗.

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

ξ(λ)=dNNdlogλ=ξ∗λ(λ/λ∗)δ1+(λ/λ∗)δ2,

(9)

where we set δ1=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 δ1<−αBH.

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 γ2) and thereafter decay exponentially with a constant decay time τ. In this case, the concept of ”duty cycle” fd 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

ξ∗λ=η⋅τ,

(10)

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

fd=ξ∗λδ2ln(10)

(11)

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

If δ1<−αBH, 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 δ1 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 mbh/m∗ relation, we can also connect the observed L∗ of the QLF to the M∗ of the star-forming galaxy mass function with

L∗=1038.1⋅λ∗⋅M∗⋅(mbhm∗)⋅ΔL(δ1,γ2).

(13)

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 M∗ of the galaxy population (M∗∼1010.8M\sun) and Φ relative to the ϕ∗SF of the galaxy population. All luminosities are relative to the L∗ 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 M∗ and ϕ∗SF of the input Schechter mass function of the star-forming galaxies. These distributions are plotted for different AGN luminosities relative to the knee luminosity L∗ 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 γ2+αSF, i.e. with αSF∼−0.4 and γ2∼2 (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 δ1, 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 L∗, 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 M∗. Above L∗, the peak host galaxy mass is always close to M∗ because there are so few more massive galaxies. The numbers of M∗ 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 (αSF+γ2, with γ2∼2) that of the passive galaxies which is given by αSF+1 (see Peng et al. (2010)). In mathematical form, the mass function of galaxies at any AGN luminosity above L∗ will have a Schechter shape

ϕ(m∗;L>L∗)∝(m∗M∗)αSF+γ2exp(−m∗M∗).

(14)

At lower luminosities, well below L∗, the behaviour changes and a region of intermediate slope appears. For black hole masses between (1038.1)−1(mbh/m∗)−1(λ∗)−1L and (1038.1)−1(mbh/m∗)−1(λ∗)−1L∗, the mass function will have the slope given by α+δ1. At very low masses we will again always recover the slope of αSF+γ2.

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 α+δ1. 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

mpeak∼M∗1038.1mbhλ∗L|αSF|>|δ1|mpeak∼M∗1038.1mbhλ∗L∗|αSF|<|δ1|.

(15)

This difference in behaviour can be understood as follows: In the |αSF|>|δ1| 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 mbh. 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 M∗) 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 δ1 discussed above, which can shift the peak of ϕ(m).

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, ϕBH(mbh), 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 ϕ∗SF and M∗ parameters in the Schecter mass function of star-forming galaxies. Data shows minimal evolution of M∗ until at least z≈ 2.5, but significant evolution in ϕ∗SF.

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 (ϕ∗SF≡ϕ∗SF(z), M∗≡M∗(z)) 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

ϕSF(m∗,z)=dNdlogm∗=ϕ∗SF(z)(m∗M∗(z))αSFexp(−m∗M∗(z)),

(16)

where we model the redshift dependence as

logϕ∗SF(z)=a0+a1κ+a2κ2+a3κ3,logM∗(z)=a0+b1κ+b2κ2+b3κ3,

(17)

with κ≡log(1+z).
The slope of the low mass end, αSF, was kept constant at the local value of αSF=−0.4 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).

The evolution of parameters of mass function in Figure 3 confirms that M∗ 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 M∗ that does not change with redshift. It is important to stress that our results will not depend critically on the exact functional evolution of M∗ above z≳2.5, but will use the by now well established observational fact that M∗ does not change at z≲2.5 and that the evolution in the galaxy population since that time is associated with a “vertical” evolution in ϕ∗SF. 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 0<z<6. 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. ϕ∗QLF, L∗, γ1 and γ2) 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 ϕ∗QLF at a constant value. We now know that ϕ∗SF of the galaxy population changes significantly over the redshift range of interest and, as developed in the previous section, the relationship of ϕ∗QLF relative to ϕ∗SF is of great interest in the context of the duty cycle. In contrast, γ1 is set, in our convolution model, by the low mass slope αSF 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 γ1. Therefore, we choose to have a redshift-independent faint end slope of the luminosity function, γ1 and to allow ϕ∗QLF 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 γ2 to vary. We do this because, as we have seen, γ2 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 ϕ∗QLF is set to be exactly the redshift dependence of ϕ∗SF.

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 ϕ∗QLF and L∗ 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 L∗ and ϕ∗QLF 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 ϕ∗QLF with redshift. Bottom left: Projection showing explicitly change of L∗ with redshift. Figure 6.— Same as Figure 5 for QLF fit at redshifts z≤3, in which we demanded that ϕ∗QLF∝ϕ∗SF, 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 ϕ∗SF which were determined by fitting the data from Ilbert et al. (2013) at a single redshift. Black squared show the values of ϕ∗QLF 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 ϕ∗SF/ϕ∗QLF. 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 ϕ∗QLF

The first result is that we notice a strong similarity between the observed redshift dependence of ϕ∗QLF in Figure 5 and the observed ϕ∗SF 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 ϕ∗QLF and ϕ∗SF 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 ϕ∗SF.

Referring back to Section 3.3, a constant ratio between ϕ∗QLF and ϕ∗SF 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 ϕ∗QLF to be exactly that of ϕ∗SF, 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 ϕ∗SF available (z<3.5), while the constant multiplicative offset between these parameters remains a free parameter to be determined by the fitting procedure, i.e.

ϕ∗QLF∝ϕ∗SF.

(19)

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) ( χ2thiswork/d.o.f.=2.1; χ2Hopkins/d.o.f.=2.0), 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 ”ϕ∗SF-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 χ2/d.o.f.=2.1. 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 ϕ∗SF(z) and ϕ∗QLF(z) track each other throughout cosmic time (at least since z∼3, 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”, fd, 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 L∗ 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 L∗qlf

The other striking feature of the quasar luminosity function is the strong redshift evolution of L∗, which increases by almost two orders of magnitude back to z∼2. At higher redshifts, the L∗(z) 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 ϕ∗SF(z) and ϕ∗QLF(z) tracking each other.

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

As discussed in Section 4.1, we know that the characteristic M∗ of the galaxy population does not change, especially in the redshift range where increase of L∗ is most prominent (z≤2), 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.1. The quasar mass-luminosity plane

Figure 8.— Mass-luminosity plane shown in 3 representative redshift bins. Our model with assumed non-evolving relation mbh/m∗=10−2.8 is in the top row and our model with assumed evolving relation mbh/m∗=10−3(1+z)2 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 mbh/m∗ scaling relations and an assumed scatter. We will then assign Eddington ratios from the evolving ξ(λ,z) 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 mbh/m∗ and the redshift evolution in λ∗ must combine to produce the observed evolution in the QLF L∗.

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 mbh/m∗ and ξ(λ,z) 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 z=2, 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 (mbh,L) plane in 3 representative redshift bins if we assume a redshift independent mbh/m∗ with the standard value of mbh/m∗≈10−2.8. We introduce a log-normal scatter in this relationship of 0.5 dex to account for scatter in mbh/m∗ relationship. The solid diagonal line indicates the Eddington limit (λ=1).

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 mbh/m∗ 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 L∗ is much higher then locally and is a reflection of the high λ∗ implied for an unevolving mbh/m∗ relation shown in Figure 17.

It should be noted that the problem gets worse if we reduce the scatter in the mbh/m∗ 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 mbh/m∗ 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 mbh/m∗ relation. We adopt as a heuristic example the form mbh/m∗∝(1+z)n with n=2. 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 ∼(1+z)4 increase in L∗ back to z∼2 would be due to an equal split between a (1+z)2 change in mbh/m∗ and a (1+z)2 change in characteristic Eddington ratio λ∗, remembering that these two changes are degenerate in our convolution model without the mbh/m∗ data of Figure 8.

Finally we note that, quite independently of any assumptions about mbh/m∗, 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 L 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 M∗, 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 (mbh,L) 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 mbh/m∗ correlation in quenched systems

In this section we show how it is possible to reproduce the observed tight correlation between mbh/mbulge 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 mbulge 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 mbulge in what follows with mquench.

Figure 10.— Top left: Resulting mbh/m∗ relation in quenched galaxy systems today, if we only take into account systems that have quenched after z=2. 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 mbh/m∗ relation in quenched systems today, for full range of quenching history (z<5). Orange and blue line are generated in same way as in previous panel. Merging, uncertain evolution in mbh/m∗ relation and errors in estimated rate of quenching at high redshift could bias this result, as discussed in text. Bottom left: Measured mbh/mbulge by Kormendy & Ho (2013). Orange line shows again best fit quoted in Kormendy & Ho (2013). Bottom right: Measured mbh/mbulge 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 mbh/m∗ 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 z=2 and then include all galaxies that have quenched since z=5.
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 z<2, may break down at higher redshifts. When we ignore quenching at 2<z<5, 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

100(mbhm∗,quench)=a⋅(m∗,quench1011Msun)b

(20)

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 (a,b)=(0.40,0.09) and a scatter of 0.53 dex for the case of quenching only from z=2, and values of (a,b)=(0.63,0.18) with scatter of 0.6 dex for the case of quenching from z=5. These values should be compared with the observed values of (a,b)=(0.49+0.06−0.05,0.14±0.08) 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 mbh/mbulge relation is the underlying mbh/m∗ 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 mbh/mbulgerelation is positioned roughly at the value of mbh/m∗ at the mean redshift of quenching at that mass, which is z∼1−1.5
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 mbh/m∗ value, because of the evolution in mbh/m∗. 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 mbh/mbulge 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 mbh/mbulge relation (b≠0) even though the input mbh/m∗ 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 (1+z) so that the usual Faber-Jackson type scaling relation would be expected to evolve as

mstar∝σ4(1+z)2.

(21)

If mbh/m∗ scales as (1+z)2, as we have been exploring in this part of the paper, then this naturally produces an mbh−σ relation that is independent of redshift.

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

We stress that this tighter mbh−σ 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 mbh/m∗ is in fact caused by the constancy of an underlying mbh−σ causal connection. On the other hand, there could well be other causes for mbh/m∗ to evolve as (1+z)2, in which case the tight mbh−σ relation would simply be a coincidence.

6.3. mbh/m∗ in AGNs in the local Universe

Finally, we turn to estimates of the mbh/m∗ 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 z<0.6 quasars in the SDSS Stripe 82. They found that the quasars are predominately hosted in massive star-forming galaxies, with relatively large mbh/m∗ ratios of around 10−2.5. 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 η=0.4 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 mbh/m∗ 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 mbh/m∗ relation in Stripe 82, mimicking analysis of Matsuoka et al. (2014). Blue line shows mean mbh/m∗ relation in full sample (before luminosity cut was applied), while the shaded area shows 1−σ 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 10−2.8, 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 mbh/m∗∝(1+z)2 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.

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 mbh/m∗ and Eddington ratio distribution, ξ(λ), successfully reproduces the observed quasar luminosity function (by construction) and also reproduces observations of the distribution of quasars in the (mbh,L) plane, the black hole to bulge mass relation of quenched galaxies and measurements of mbh/m∗ for low redshift AGN. In this Discussion section, we develop further some of the astrophysical implications of the preceding results.

It is important to appreciate that a quite rapid evolution in mbh/m∗ in active (star-forming) AGN systems does not imply an equally rapid evolution of mbh/m∗ 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 mbh/m∗ 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 z∼2 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 mbh/m∗ 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 mbh nor m∗ will have changed once a given object becomes inactive (i.e. passive/quenched). It is therefore clear that the mean mbh/m∗ 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 mbh/m∗ of star forming systems at z ∼ 1-1.5 to the mbh/m∗ ratio of passive galaxies seen today that quenched at z ∼ 1-1.5, then we would expect to see no change in mbh/m∗, 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 mbh/m∗ 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 mbh/m∗ scaling as (1+z)2 for active systems, then we would expect to see (1+z)0.8 if we compare comparisons of active systems at z∼2 with quenched systems at z∼0, 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 mbh−σ 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 (1+z)2 dependence in mbh/m∗.

Figure 12.— Evolution of the mean mbh/m∗ 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 mbh/m∗ 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 mbh/m∗ upwards in luminosity selected active samples, which goes in the opposite direction. This is indicated by the thin black lines which show the observed mbh/m∗ for active systems for luminosity-selection (labelled relative to L∗) 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. mbh/m∗∼(1+z)n. However, this is an arbitrary redshift dependence, and there are reasons why it cannot hold (with n∼2) 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),

rsSFR(z)=˙m∗m∗=−1(1+z)√ΩM(1+z)3+ΩΛ1m∗dm∗dz,

(22)

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

m∗(z)=m∗(0)(exp[2⋅0.07√ΩM+ΩΛ3H0ΩM−2⋅0.07√ΩM(1+z)3+ΩΛ3H0ΩM]).

(23)

This modest increase in stellar mass for a star-forming galaxy sets a maximal change in mbh/m∗ ratio even in the most extreme case that mbh does not increase at all. We show this maximal evolution in Figure 13 and compare it with the (1+z)2 evolution used elsewhere in the paper. It can be seen that the maximal evolution is actually slower then (1+z)2 at redshifts z≲0.7. Galaxies are growing so slowly at low redshifts that it is not possible to create a strong evolution in mbh/m∗ 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 mbh/m∗. The maximal change in mbh/m∗ 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 mbh/m∗∝(1+z)2 should best be interpreted as implying a change of a factor of about ten to z∼2, rather than a precise particular dependence on redshift.

Figure 13.— Comparison of the maximal evolution of the mbh/m∗ relation in a star-forming galaxy that stays on the Main Sequence, with our assumed mbh/m∗∝(1+z)2. At low redshift galaxies are growing so slowly in the stellar mass that strong evolution in the mbh/m∗ 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 ξ(λ,z). 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 L∗ 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 mbh/m∗ mass ratio and characteristic Eddington ratio, λ∗. The second is the redshift evolution of ϕ∗QLF that changes the number densities uniformly with luminosity. We have argued in this work that the evolution of ϕ∗QLF is a direct consequence of the evolution of ϕ∗SF, 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 ϕ∗SF and ϕ∗QLF

One of the most striking results of this paper is that the observed evolution of ϕ∗QLF of the QLF tracks the observed evolution of ϕ∗SF 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, ηm, which may be written

ηm∼sSFR(z)⋅(m∗/M∗),

(24)

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)

sSFR∼{(1+z)3,when z<2(1+z)5/3,when z>2.

(25)

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 ηAGN (from Equation 10) to be equal to ηm. These can only be equated for our inferred constant ξ∗λ for a particular redshift and mass dependence of the decay time τ.

τ(m,z)∼sSFR(z)−1(m/M∗)−1ξ∗λ.

(26)

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 ηm and ηAGN.

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 ξ(λ,z), while the former is itself a convolution of the galaxy mass-function with the mbh/m∗ 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:

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, L∗, is proportional to the product of the M∗ of the galaxy mass function, the ratio mbh/m∗ and the position λ∗ of the break in the Eddington ratio distribution while the ϕ∗QLF normalization of the QLF is proportional to the product of the ϕ∗SF 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 M∗) for different AGN luminosities (relative to L∗). At high luminosities (above the AGN L∗) this is a Schechter function with the star-forming M* but a faint end slope given by αSF+γ2∼1.5.

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

In contrast, the QLF L∗ evolves strongly with redshift, with evolution being at least L∗∝(1+z)3 up to z∼2. Given that there is strong evidence for no change in the galaxy M∗ over this redshift range, this evolution in L∗ 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, mbh/m∗, 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 mbh/m∗∝(1+z)2 in star-forming systems. This implies that the observed ∼(1+z)4 increase in L∗ back to z∼2 would be due to an equal split between a (1+z)2 change in mbh/m∗ and a (1+z)2 change in characteristic Eddington ratio λ∗.

We show that this is compatible with the observed mbh/m∗ 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