Mon. Not. R. Astron. Soc. 000, 1–31 (2014)

Printed 3 December 2014

(MN LATEX style file v2.2)

arXiv:1412.0754v1 [astro-ph.GA] 2 Dec 2014

The cosmic growth of the active black hole population at
1 < z < 2 in zCOSMOS, VVDS and SDSS
A. Schulze1,2⋆, A. Bongiorno3, I. Gavignaud4, M. Schramm1, J. Silverman1, A. Merloni5,
G. Zamorani6, M. Hirschmann7, V. Mainieri8, L. Wisotzki9, F. Shankar10, F. Fiore3,
A. M. Koekemoer11, G. Temporin12
1 Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
2 Kavli Institute for Astronomy and Astrophysics, Peking University, 100871 Beijing, China
3 INAF-Osservatorio Astronomico di Roma, Via di Frascati 33, 00040, Monteporzio Catone, Rome, Italy
4 Departamento de Ciencias Fisicas, Universidad Andres Bello, Av. Republica 220, Santiago, Chile
5 Max-Planck-Institut für extraterrestrische Physik (MPE), Giessenbachstrasse 1, Garching bei München, D-85748, Germany
6 INAF-Osservatorio Astronomico di Bologna, Via Ranzani 1, 40127, Bologna, Italy
7 UPMC-CNRS, UMR7095, Institut d’Astrophysique de Paris, Boulevard Arago, F-75014 Paris, France
8 ESO, Karl-Schwarschild-Strasse 2, Garching bei München, D-85748, Germany
9 Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
10 School of Physics & Astronomy, University of Southampton, Southampton SO17 1BJ, UK
11 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218
12 Institute of Astro- and Particle Physics, University of Innsbruck, Technikerstrasse 25, A-6020 Innsbruck, Austria

Accepted 2014 November 28. Received 2014 November 19; in original form 2014 October 6

ABSTRACT

We present a census of the active black hole population at 1 < z < 2, by constructing
the bivariate distribution function of black hole mass and Eddington ratio, employing
a maximum likelihood fitting technique. The study of the active black hole massfunction (BHMF) and the Eddington ratio distribution function (ERDF) allows toclearly disentangle the AGN downsizing phenomenon, present in the AGN luminosity
function (AGN LF), into its physical processes of black hole mass downsizing andaccretion rate evolution. We are utilizing type 1 AGN samples from 3 optical surveys(VVDS, zCOSMOS and SDSS), that cover a wide range of 3 dex in luminosity over ourredshift interval of interest. We investigate the cosmic evolution of the AGN population
as a function of AGN luminosity, black hole mass and accretion rate. Compared toz = 0 we find a distinct change in the shape of the BHMF and the ERDF, consistent
with downsizing in black hole mass. The active fraction or duty cycle of type 1 AGN atz ∼ 1.5 is almost flat as a function of black hole mass, while it shows a strong decrease
with increasing mass at z = 0. We are witnessing a phase of intense black hole growth,which is largely driven by the onset of AGN activity in massive SMBHs towards z = 2.
We finally compare our results to numerical simulations and semi-empirical models
and while we find reasonable agreement over certain parameter ranges, we highlightthe need to refine these models in order to match our observations.
Key words: Galaxies: active - Galaxies: nuclei - quasars: general - quasars: supermassive black holes

1

INTRODUCTION

Supermassive black holes (SMBHs) constitute a fundamental component of galaxies, as almost every massive galaxy
habours a SMBH in its centre (Kormendy & Richstone
1995). A complete picture of galaxy evolution requires⋆

E-mail: andreas.schulze@ipmu.jp

c 2014 RAS

an understanding of the growth history of supermassive black holes. This is demanded by observational
evidence for co-evolution between the SMBH and its
host galaxy, suggested by the observed tight correlation between SMBH mass and the properties of thegalaxy’s spheroidal component, e.g. with stellar velocity
dispersion (e.g. Ferrarese & Merritt 2000; Gebhardt et al.
2000; Tremaine et al. 2002; McConnell & Ma 2013),

2

Schulze et al.

bulge luminosity and bulge mass (Magorrian et al. 1998;
Marconi & Hunt 2003; Häring & Rix 2004). Furthermore,
the cosmic black hole accretion rate density traces well thestrong evolution in the cosmic star formation rate density
since z ∼ 3 (e.g. Boyle & Terlevich 1998; Marconi et al.
2004; Silverman et al. 2008), implying a link between starformation and black hole accretion in a global sense. Sucha connection between star formation and AGN activity
is also directly seen for luminous AGN (e.g. Netzer et al.
2007; Lutz et al. 2008; Silverman et al. 2009; Rosario et al.
2012).
The Soltan argument suggests that most black hole
growth takes place in luminous AGN phases (Soltan 1982;
Salucci et al. 1999; Yu & Tremaine 2002). Therefore, thedemographics of the AGN population enables an assessmentof the cosmic SMBH growth history. The main demographic
quantity is the AGN luminosity function (AGN LF), whichis established over a wide range in redshift, luminosity
and wavelength (e.g. Boyle et al. 2000; Wolf et al. 2003;
Ueda et al. 2003; Hasinger et al. 2005; La Franca et al.
2005; Richards et al. 2006a; Bongiorno et al. 2007;
Silverman et al. 2008; Croom et al. 2009; Aird et al.
2010; Assef et al. 2011; Fiore et al. 2012; Ueda et al. 2014).
The evolution of the AGN LF allows to constrain the blackhole growth history using a continuity equation to evolve theblack hole mass function through redshift (Yu & Tremaine
2002; Yu & Lu 2004; Marconi et al. 2004; Merloni 2004;
Merloni & Heinz 2008; Yu & Lu 2008; Shankar et al. 2009;
Draper & Ballantyne 2012). These studies found that the
local total black hole mass function is consistent with beinga relic of previous AGN activity. They also link the AGN”downsizing”, observed in the AGN LF (Ueda et al. 2003;
Hasinger et al. 2005; La Franca et al. 2005; Croom et al.
2009), to the growth of SMBHs, i.e. the most massive
black holes grow at earlier epochs, while lower mass black
holes are still actively growing at later times, a behavior
which is also called anti-hierachical black hole growth(Marconi et al. 2004; Merloni 2004).
However, having only information on the AGN luminosity does not allow us to distinguish the two quantities that govern SMBH growth i.e. SMBH mass (M• )and accretion rate. The latter can be expressed as the
Eddington ratio λ = LBol /LEdd , with the Eddington luminosity LEdd = 1.3 × 1038 (M• /M⊙ ) erg s−1 . This isa critical limitation, since AGN show a wide distributionof accretion rates (Heckman et al. 2004; Babić et al. 2007;
Kauffmann & Heckman 2009; Schulze & Wisotzki 2010;
Aird et al. 2012; Bongiorno et al. 2012). Additional assumptions on the accretion rate distribution or alternatively themean accretion rate are required in the above studies to link
the bolometric AGN luminosity to black hole growth. Thus,
while the AGN LF is relatively easy to determine observationally, its information content is limited.
This situation is similar to galaxy evolution studies,
where the galaxy LF already provides important information, but the study of the stellar mass function in connection
to the specific star formation rate (sSFR) largely enhances
our understanding (e.g. Elbaz et al. 2007; Peng et al. 2010).
Galaxies show a downsizing trend, equivalent to black holes,in the sense that for more massive galaxies the bulk of theirmass buildup happens at earlier cosmic times compared toless massive galaxies (e.g. Cowie et al. 1996; Bundy et al.

2006; Pérez-González et al. 2008). While for galaxy evolution, stellar mass and sSFR are key parameters, the same istrue for black hole growth with black hole mass and specificaccretion rate (i.e. Eddington ratio).
To break the degeneracy in AGN luminosity and disentangle the AGN LF into their underlying distribution functions, additional information on the SMBH mass has to be
added. The BHMF and ERDF provide additional observational constraints on models of galaxy evolution, black hole
growth and black hole - galaxy co-evolution using empirical approaches (Shankar et al. 2013), semi-analytic models(Fanidakis et al. 2012; Hirschmann et al. 2012) or numerical simulations (Di Matteo et al. 2008; Hirschmann et al.
2014; Sijacki et al. 2014). The ERDF also allows further constraints on black hole accretion models and the triggeringhistory of black holes (Yu et al. 2005; Hopkins & Hernquist
2009).
Additional constraints on black hole - galaxy
co-evolution are obtained from the cosmic evolution of the relation between M• and bulge properties(Croton 2006; Booth & Schaye 2011; Dubois et al. 2011;
Anglés-Alcázar et al. 2013). However, its observational
determination is affected by sample selection effects(Lauer et al. 2007; Schulze & Wisotzki 2011). To modeland account for selection effects that will bias observations,precise knowledge of the underlying distribution functions,
including the BHMF, is required (Schulze & Wisotzki 2011;
Schulze & Wisotzki 2014).
The observational determination of BHMF and ERDF
requires a well-defined AGN sample, as it is the case forthe luminosity function determination, but in addition requires a black hole mass estimate for every object in thesample. The most direct black hole mass estimates for large
statistical AGN samples use the ”virial” method to estimate M• from a single-epoch spectrum of a type 1 AGN
with broad emission lines in their optical spectra. We here
restrict our investigation to these broad line AGN and emphasize that our definition of an active black hole only refers
to this unobscured AGN population (we will try to correctfor obscured AGN in section 5.2).
With calibrations available for the broad Hα, Hβ, Mg ii
and C iv line (e.g. Vestergaard 2002; McLure & Dunlop
2004; Greene & Ho 2005; Vestergaard & Peterson 2006;
Trakhtenbrot & Netzer 2012; Park et al. 2013), this enables the estimate of M• for large AGN samples out tohigh-z (e.g. McLure & Dunlop 2004; Netzer & Trakhtenbrot
2007; Shen et al. 2008; Trump et al. 2009; Shen et al. 2011).
However, these ”virial” black hole masses are still subject to significant uncertainties and possible systematics(e.g. Vestergaard & Peterson 2006; Shen 2013; Denney et al.
2013). While these are particularly important for individual
objects, they are less critical for demographic studies of the
AGN population, as will be carried out here.
The establishment of this ”virial” method for blackhole mass estimation allowed the empirical determinationof the active BHMF. There have been two main approaches
used for the determination of the BHMF and ERDF. The
first is directly borrowed from the determination of the
AGN LF, using the classical 1/Vmax estimator (Schmidt
1968) with the same volume weights as for the AGN LF(Wang et al. 2006; Greene & Ho 2007; Vestergaard et al.
2008; Vestergaard & Osmer 2009; Schulze & Wisotzki 2010;
c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2Shen & Kelly 2012; Nobuta et al. 2012). We refer to this approach as luminosity weighted Vmax method. While it provides a non-parametric and model independent estimate ofthe distribution function, in general it suffers from severe
incompleteness due to the improper volume weights which
do not correct for active black holes below the flux limit of
the survey (Kelly et al. 2009; Schulze & Wisotzki 2010, seealso Appendix B). Furthermore, it does not account for theuncertainty in the virial M• estimates. Therefore, in general
it does provide a biased census of the active SMBH population.
For a proper treatment of the survey selection function
the BHMF and ERDF need to be determined jointly by fitting an analytic model for the bivariate distribution functionof M• and λ to the observations. This forward modeling of
the distribution functions can be achieved via a Maximum
likelihood method, as developed by Schulze & Wisotzki(2010), or within a Bayesian framework (Kelly et al. 2009).
Schulze & Wisotzki (2010) used the Maximum likelihood method to determine the BHMF and ERDF at z < 0.3
from the Hamburg/ESO Survey (Wisotzki et al. 2000), establishing the local zero point of these distribution functions.
Their results support the notion of downsizing in black holemass. However, a full exploration of the growth history of
black holes requires the determination of BHMF and ERDF
at higher redshifts. The SDSS QSO sample has been used
to study the evolution of black hole activity by determining
the BHMF for 0.4 < z < 5 (Kelly et al. 2010; Shen & Kelly
2012; Kelly & Shen 2013), employing a Bayesian framework(Kelly et al. 2009). These authors also found evidence forblack hole mass downsizing in the BHMF and for z > 2 also
in the ERDF (Kelly & Shen 2013). However, since the SDSS
only covers the bright end of the AGN LF, their results arelargely limited to the high mass end of the BHMF and tohigh Eddington ratios and are affected by large incompleteness at lower values, which introduces large uncertainties
in this regime. To resolve these limitations and probe the
BHMF and ERDF to lower masses and λ, deep AGN surveys have to be used. Nobuta et al. (2012) presented the
BHMF and ERDF at z ∼ 1.4 based on an X-ray selected
sample from the Subaru XMM-Newton Deep Survey (SXDS,
Ueda et al. 2008), utilizing the Maximum likelihood methodof Schulze & Wisotzki (2010). Since deep surveys with spectroscopic follow-up are currently limited to small areas, they
did not probe the bright end of the luminosity function.
Therefore, for a wide luminosity coverage and thus a wide
coverage of black hole mass and Eddington ratio it is important to combine deep, small area surveys with shallower,large area surveys.
In this paper we follow this strategy. We are using two
well-defined AGN samples from the deep, small area surveysof the VVDS and zCOSMOS, and combine our results withthe large area, shallow SDSS to constrain the BHMF and
ERDF within the redshift range 1.1 < z < 2.1, covering awide range in M• and λ. Our redshift range is chosen suchthat it covers the broad Mg ii line in the spectra for all threesamples.
In Section 2 we present the individual samples we areusing. The black hole masses and Eddington ratios of thethree samples are presented in Section 3. In Section 4 we
present our results for the bivariate distribution function.
We discuss these results in Section 5 and conclude in Secc 2014 RAS, MNRAS 000, 1–31

3

tion 6. The Appendix gives more details on the selection
functions of the employed samples and present results on the
BHMF, ERDF and AGN LF based on the Vmax method.
Throughout this paper we assume a Hubble constant of
H0 = 70 km s−1 Mpc−1 and cosmological density parameters Ωm = 0.3 and ΩΛ = 0.7.

2

THE SAMPLE

2.1

The VVDSThe VVDS is a pure magnitude limited spectroscopic survey(Le Fevre et al. 2005, 2013; Garilli et al. 2008). The target
selection is purely based on the I band magnitude without
any further selection criteria. It consists of a wide survey(Garilli et al. 2008), spread over 3 fields, with an effective
area of 4.5 deg2 for 17.5 < IAB < 22.5 and a deep survey(Le Fevre et al. 2005), spread over 2 fields, with an effective area of 0.62 deg2 for 17.5 < IAB < 24.0. The I band
imaging has been designed to be complete to much faintermagnitudes than these limits, with limiting magnitudes of
IAB = 24.8 and IAB = 25.3 at 5σ in the wide and the deep
fields, respectively (McCracken et al. 2003; Le Fèvre et al.
2004). Targets have been nearly randomly selected for spectroscopy with the VIMOS multi object spectrograph on the
VLT. The survey used the LRRED grism covering the spectral range 5500 − 9500 Å at a resolution of R ∼ 230, corresponding to a spectral resolution element of ∼ 33 Å. Type 1
AGN have been manually identified based on the presence ofbroad emission lines (FWHM> 1000 km s−1 ) in their spectra. These simple selection criteria provide a well-defined
AGN sample, free of a potential bias by pre-selection via
color and/or morphology. Details on the type 1 AGN sample
are given in Gavignaud et al. (2006) and Gavignaud et al.
(2008), while the type 1 AGN luminosity function, basedon the ”epoch 1” sample is presented in Bongiorno et al.
(2007).
We here use the AGN sample from VVDS ”epoch 2”,
presented in Gavignaud et al. (2008). It consists of 298
AGN, 235 with a secure redshift and 63 with a degenerate redshift, because only a single broad emission line without any other feature is present in the spectrum. Black holemasses based on Mg ii have been measured for 86 out of 95
within 1.0 < z < 1.9, for which the data quality was sufficient to estimate M• , and are presented in Gavignaud et al.
(2008). These 86 AGN form our VVDS-sample. We additionally included the redshift degenerate objects that could lie
in our redshift range (61 AGN), measured their M• , assuming the broad line is Mg ii, and statistically account for theircontribution to the full sample (see Appendix A).

2.2

zCOSMOS

Like the VVDS, zCOSMOS is a magnitude limited spectroscopic survey (Lilly et al. 2007) using VIMOS. It consistsof zCOSMOS-bright (Lilly et al. 2009), targeting ∼ 20000
galaxies over the full COSMOS field and zCOSMOS-deep,targeting ∼ 10000 galaxies over the central 1 deg2 .
We here use the dataset of the zCOSMOS-bright
20k spectroscopic catalog. We do not include the type-1
AGN detected in the deep sample (Bongiorno et al. 2012;

Schulze et al.

log Lbol [erg s−1 ]

4ple (Richards et al. 2006a; Shen & Kelly 2012). The sample
is magnitude limited to 15.0 < i < 19.1 below z = 2.9 andthe target selection is based on multi-colour information.
The effective area of the survey is 6248 deg2 (Shen & Kelly2012).
We further restrict the uniform QSO sample to the redshift range 1.1 < z < 2.1, containing 28322 QSOs. Black holemasses, based on the broad Mg ii line have been estimatedby Shen et al. (2011). We reject M• estimates with measurement error larger than 0.5 dex, leading to a sample of 27257
QSOs with available M• , forming our SDSS-sample.

47

46

45

44
1.2

1.4

1.6

z

1.8

2.0

Figure 1. Redshift versus bolometric AGN luminosity for the 3
AGN samples within the redshift range considered in this work,
1.1 < z < 2.1. The orange circles, blue squares and the red contours are for the VVDS, zCOSMOS and SDSS sample respectively(the VVDS sample is restricted to 1.1 < z < 1.9 and redshift degenerate objects are shown as open circles).

In Figure 1 we show the redshift distribution against luminosity for the three samples within the redshift range westudy in this work (1.1 < z < 2.1).
The study of the demographics of the AGN populationfrom these well-defined samples requires knowledge on therespective selection functions. In Appendix A we provide details on the applied selection functions for the three surveys.

3
Rosario et al. 2013), since they do not constitute a welldefined sample, mainly due to the BzK colour preselection, targeting z > 1.5 galaxies. Targets have been chosenfrom HST F814W imaging (Koekemoer et al. 2007) with15 < IAB < 22.5 across the 1.7 deg2 COSMOS field. The
target selection is based on a combination of nearly random
selection and compulsory targets. The latter were observedwith a higher chance and were largely based on X-ray detection (Hasinger et al. 2007; Brusa et al. 2010), which strongly
favors AGN. For the statistical analysis this has to be taken
into account by proper weighting (see Appendix A). Spectra
have been obtained with the VIMOS multi object spectrograph within the range 5500 − 9500 Å at a resolution ofR ∼ 580.
Type 1 AGN are again identified by the presence ofbroad emission lines with FWHM> 1000 km s−1 in the
spectra. Results on the zCOSMOS AGN sample from the
10k catalog are presented in Merloni et al. (2010). Blackhole masses from broad Mg ii can be estimated within theredshift range 1.1 < z < 2.1. zCOSMOS-bright contains
160 broad line AGN within this redshift range, for 145 ofthem we could measure a reliable M• . These AGN form ourzCOSMOS-sample.
We note that the COSMOS field has a rich multi wavelength coverage from X-rays to radio. This gives us for example a more precise handle on the bolometric luminosities
and the host galaxy contribution of the AGN sample.
2.3

SDSS

For the SDSS sample we largely follow the work ofShen & Kelly (2012). We start from the SDSS DR7 QSO
sample (Schneider et al. 2010), consisting of 105,783 QSOs.
We restrict this sample to the subset that has been selectedby a uniform target selection algorithm (Richards et al.
2002) and thus form a well-defined broad line AGN sam-

3.1

BLACK HOLE MASSES AND EDDINGTONRATIOSBlack hole masses

Black hole mass estimates from single epoch broad line
AGN spectra can be obtained via the ”virial method”. Thismethod employs the virial relation M• = f RBLR ∆V 2 /G,
where the broad line cloud velocity ∆V is derived fromthe broad emission line width, f is a scale factor that accounts for the geometry and kinematics of the BLR (usually calibrated by scaling to the local M• −σ∗ relationship(e.g. Onken et al. 2004; Woo et al. 2010; Woo et al. 2013)),and RBLR is the broad line region size. While the lattercan be directly measured via reverberation mapping (e.g.
Blandford & McKee 1982; Peterson et al. 2004), for single
epoch spectra, commonly the thereby established scaling
relation to the continuum luminosity (Kaspi et al. 2000;
Bentz et al. 2009, 2013) is used.
We (re-)computed black hole masses consistently for all
3 samples using the formula from Shen et al. (2011):
0.62 
2
L3000
FWHM
6.74
M• = 10M⊙ (1)
1044 erg s−1
1000 km s−1
The Mg ii broad line widths and the continuum luminosities at 3000 Å have been measured for the SDSS sample
by Shen et al. (2011). For the VVDS sample they are provided by Gavignaud et al. (2008), but only for AGN with secure redshift. For the redshift degenerate AGN in the VVDS
sample, we measured the Mg ii FWHM and L3000 using exactly the same procedure as for the rest of the VVDS sample,
assuming the broad line in the spectrum is indeed Mg ii.
For the zCOSMOS sample we used the VIMOS spectra
to fit the broad lines. A few of the objects also have SDSS
spectra and we find consistent results for these. The fitting
procedure is similar to Gavignaud et al. (2008) and outlinedin Schramm & Silverman (2013).
For both, the VVDS and the zCOSMOS sample, thec 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 251.0

log Lbol (IAB)

46.5

LAGN/LTot

0.9

0.8

0.7

45.5

45.0

44.5
44.5

0.6

44.5

45.0

45.5

46.0

46.5

L3000,Tot
Figure 2. Average host galaxy contamination correction for thecontinuum luminosity at 3000 Å, derived from the SED decomposition of COSMOS type 1 AGN. The black squares show theaverage correction for the sample in 4 bins of total L3000 , andthe dashed line represents the interpolation for the average correction. We applied this correction to L3000 for the VVDS andzCOSMOS samples.

Mg ii line region is modeled by a pseudo continuum, consisting of a power law and the broadened Fe ii template from
Vestergaard & Wilkes (2001) and one or two Gaussian components for the broad line (always two Gaussians for the
VVDS AGN). The Mg ii broad line FWHM was measured
from the best fit line profile and the monochromatic continuum luminosity at 3000 Å L3000 is taken from the power law
fit. Details on the black hole mass measurement for COSMOS AGN will be given in Schramm et al. (in prep.; seealso Merloni et al. 2010).
For the SDSS AGN sample the host galaxy contamination to L3000 is negligible at these redshifts (Shen et al.
2011). However, this is not the case for the deeper VVDS andzCOSMOS samples, thus L3000 needs to be corrected for thiscontamination. For these samples we applied an average hostgalaxy correction to L3000 . We used the SED decompositionof 428 type 1 AGN in the COSMOS field (Bongiorno et al.
2012) to compute the average host correction in several binsof total L3000 and interpolated between the bins. The resulting average host correction as a function of L3000 is shownin Figure 2 and was applied to the continuum luminosity at
3000 Å for the VVDS and zCOSMOS AGN sample. While
we may over- or underestimate the host contribution for
individual objects, this will statistically account for the systematic effect of the host contribution. While the zCOSMOSAGN in our sample also have individual correction values,
these are consistent with this average host correction and
we decided to use the latter also for zCOSMOS, to ensure
consistency between the two samples.

Bolometric Luminosities

For the SDSS sample we compute the bolometric luminosity
from the i band magnitude, following Shen & Kelly (2012).
c 2014 RAS, MNRAS 000, 1–31

log Lbol (L3000)

44.0

45.0

45.5

46.0

46.5

log Lbol (L3000)

46.5

0.5

3.2

46.0

46.0

45.5

45.0

44.5
44.5

45.0

45.5

46.0

46.5

log Lbol (SED)
Figure 3. Comparison of bolometric luminosity estimates. The
left panels compare bolometric luminosity computed from IAB(upper panel) and from L3000 (lower panel) to the bolometricluminosity computed from the rest-frame SED for the 107 AGN
common to our zCOSMOS sample and the XMM-COSMOS sample studied by Lusso et al. (2012). The right panel compares thebolometric luminosities obtained from IAB and from L3000 witheach other. In this panel we show in addition to the zCOSMOS
sample (blue squares) the VVDS sample (orange circles). The
solid line shows the one-to-one relation.

Since the SDSS selection function is defined by its i band
flux limit, this enables a direct mapping of the selection
function from magnitude and redshift dependence to M• , λand z dependence. This largely simplifies the determinationof the BHMF and ERDF. Alternatively, either a distributionof Lbol at a given magnitude and redshift has to be assumed
to determine Ωe (Lbol , z), or an individual selection function
for every object has to be defined.
We apply the K-correction from Richards et al. (2006a)
to compute the absolute i-band magnitude normalized atz = 2 from the i band flux and translate this to the continuum luminosity at 2500 Å following Richards et al. (2006a).
We compute the bolometric luminosity from this L2500 estimate, assuming a constant correction factor fBC (2500 Å) =5 (Richards et al. 2006b; Shen & Kelly 2012). The bolometric luminosity computed this way is on average fully consistent with the bolometric luminosity obtained directly from
L3000 , by applying a constant bolometric correction factorof 5.15 to L3000 (Richards et al. 2006b; Runnoe et al. 2012).
For the VVDS and the zCOSMOS sample we follow
the same strategy of using the apparent magnitude tocompute L2500 , and derive the bolometric luminosity via
Lbol = 5L2500 . This ensures consistency between the samples and minimizes possible systematics between the samples when combining the datasets. To compute L2500 , weadopt the K-correction from Richards et al. (2006a) to compute Mi (z = 2), but transformed this K-correction to the
CFHT/CFH12K I-band (VVDS) and the HST/ACS F814W

Schulze et al.

filter (zCOSMOS), respectively, using the QSO templatefrom Vanden Berk et al. (2001). For these two filters thistransformation is small over our redshift range. We statistically correct L2500 for host galaxy contribution, using the
average host correction presented above (see Figure 2). Alternatively, we also compute Lbol from the L3000 , as measured from the spectral fitting, applying a constant bolometric correction factor fBC (3000 Å) = 5.15 (Richards et al.
2006b). We show in the upper right panel of Figure 3 thatboth estimates of Lbol give consistent results.
In both cases we assume a constant bolometric correction factor. It has been suggested that the bolometric correction might depend on luminosity (Marconi et al. 2004;
Hopkins et al. 2007; Trakhtenbrot & Netzer 2012). Whilethe bolometric correction in X-rays shows a clear dependence on luminosity or Eddington ratio (e.g. Marconi et al.
2004; Vasudevan & Fabian 2007; Jin et al. 2012; Lusso et al.
2012), this dependence is weaker for the optical bolometric correction (e.g. Marconi et al. 2004; Lusso et al. 2012;
Runnoe et al. 2012), with large uncertainties in the low luminosity regime. Therefore a constant bolometric correction
is a good approximation over the luminosity regime we areprobing.
A large fraction of our zCOSMOS sample is included in
the X-ray selected XMM-COSMOS catalogue (Brusa et al.
2010). Lusso et al. (2012) measured the bolometric luminosity directly for XMM-COSMOS AGN by integration of therest-frame SED from X-rays to 1µm. In the left panel of Figure 3 we compare our bolometric luminosities, using boththe IAB and L3000 prior to the host galaxy correction, with
their results for the 107 AGN in common. We find a good
agreement of these Lbol estimates to the direct integration
results, which still include host galaxy light, presented in
Lusso et al. (2012).

3.3

The Mass-Luminosity plane

In Figure 4 we first show the bivariate distribution of thethree AGN samples in the L3000 -FWHM plane. This gives
the distribution of the underlying observables for M• , λ andLbol . We are using the term ’distribution’, to refer to theobserved location of the sample properties in these kind ofdiagrams, while we use the term ’distribution function’ toexplicitly refer to the underlying demographic quantity, i.e.
accounting for the selection effects that affect the observed
distributions.
We here and in most other figures use the following
colour and symbol convention: (1) The VVDS sample isshown by orange circles. Filled circles are objects with secure redshift and open circles indicate degenerate redshiftobjects. (2) The zCOSMOS sample is shown by blue squares.
(3) The SDSS sample is represented by red contours.
We see the distinct luminosity regimes occupied by thebright SDSS sample and the deep VVDS and zCOSMOS
samples. For the VVDS sample there is little overlap in luminosity with the SDSS sample, while the zCOSMOS sample extends a bit more into the high luminosity regime. Thisis also illustrated in the bolometric luminosity histogram,shown in the right panel of Figure 5. It demonstrates theneed to combine deep and bright surveys for a wide luminos-

log FWHMMgII

6

4.0

3.5

3.0
43

44

45

46

47

log L3000
Figure 4. Distribution of the observable quantities 3000 Å continuum luminosity and Mg ii FWHM for the three samples. Theorange circles show the VVDS sample, where filled circles are
secure redshift objects and open circles are redshift degenerate
objects. The blue squares represent the zCOSMOS sample andthe red contours show the SDSS sample.

ity coverage. On the other hand, the FWHM distributions
for the three samples are largely consistent.1
In Figure 5 we compare the histograms of M• , λ andLbol for the VVDS, zCOSMOS and SDSS samples. The
VVDS and zCOSMOS samples have almost consistent blackhole mass distributions, covering the range 107 . M• .
109.4 and Eddington ratio distributions. Both samples extend about 1 dex lower in M• , compared to the SDSS sample, which approximately covers the range 108 . M• . 1010
and also probe lower Eddington ratios than the SDSS sample, down to 0.01 of the Eddington rate. Compared to
the SDSS sample they show a wider dispersion in theirobserved mass, luminosity and Eddington ratio distributions, consistent with previous studies on deep AGN samples(Gavignaud et al. 2008; Trump et al. 2009; Merloni et al.
2010). We emphazise that these observed distributions do
not constitute the underlying distribution functions, but are
affected by the specific survey selection criteria.
This aids the interpretation of the observed massluminosity and mass-Eddington ratio bivariate distribution,
shown in Figure 6. The red contours show the apparent bivariate distribution in M• -Lbol or M• -λ of the SDSS sample. For the studied redshift range, the SDSS flux limit corresponds to an absolute luminosity limit of ∼ 46.3 erg s−1 ,i.e. SDSS is not sensitive to AGN below this luminosity.
This becomes evident when adding the VVDS and zCOSMOS samples, which smoothly fill in the lower luminosity
parameter range missed by the SDSS. At the lowest masses

1

The zCOSMOS sample tends to show a slightly higher mean
FWHM compared to SDSS and VVDS (with hlog FWHMi = 3.70
vs. 3.61). We will investigate if this offset is intrinsic or caused bysystematics in the spectral analysis in future work.
c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2

7

60
50

N

40
30
20
10
0

7

8

9

10 -2.0-1.5

log M •

-1.0

-0.5

0.0

0.5

44

45

log λ

46

47

48

log LbolFigure 5. Histograms of black hole masses (left panel), Eddington ratios (middle panel) and bolometric luminosities (right panel) for
the VVDS (orange dashed line), zCOSMOS (blue solid line) and SDSS (red dashed dotted line). The VVDS sample includes redshift
regenerate objects with their respective weight. The SDSS sample is scaled down by a factor of 200.

0

46

log λ

log Lbol [erg s−1 ]

47

45

44

43

-1

-2

7

8

9

10

log M • [M ⊙ ]

7

8

9

10

log M • [M ⊙ ]Figure 6. Observed bivariate distribution for the three samples in the black hole mass-luminosity plane (left panel) and black holemass-Eddington ratio plane (right panel). The red contours show the apparent distribution for the SDSS sample, the blue squares are forthe zCOSMOS sample and the orange circles show the VVDS sample. For the latter, open circles indicate redshift degenerate objects.
The black solid, dashed and dotted lines indicate Eddington ratios of 1, 0.1 and 0.01, respectively.

and Eddington ratios the flux limits of VVDS and zCOSMOS also lead to an apparent absence of objects even inthese deep samples.
On the other hand, the lack of AGN at the highest
masses and highest Eddington ratios even in the SDSS (upper right corner in the right panel) is physical and a directconsequence of the decrease of the space density both in the
BHMF and the ERDF in this regime, as further shown below. The same holds for VVDS and zCOSMOS, while theirmuch smaller volumes already lack moderately high M• andhigh λ AGN, which necessitates the addition of the SDSS
data. Combined with the flux limit, this causes the apparent
anti-correlation between M• and λ in the individual surveys,c 2014 RAS, MNRAS 000, 1–31

i.e. they cannot detect low mass, low Eddington ratio AGN
below the flux limit, while high mass, high Eddington ratio
AGN do not exist in these volumes.
For the combined sample the observed distribution
is approximately confined within the Eddington limitand λ = 0.01, consistent with previous results on
type 1 AGN (e.g. Kollmeier et al. 2006; Trump et al. 2009;
Schulze & Wisotzki 2010; Shen et al. 2011). The Eddington
limit serves as an upper boundary for the accretion rate ofa black hole. The low observed number of λ < 0.01 broadline AGN could be a physical effect, caused by the change in
the accretion mode to a radiatively inefficient accretion flow
(RIAF), which leads to the disappearance of the BLR (e.g.

8Schulze et al.

Nicastro 2000; Yuan 2007; Trump et al. 2011). This changein accretion mode is thought to happen around 10−2 − 10−3 .
Alternatively the apparent limit could be driven by observational limitations to detect the corresponding very broad
and low luminosity lines beneath the spectral noise and thebroad Fe ii emission.
The study of the active black hole mass function, requires a clear definition of an active black hole. In general AGN show a broad range of levels of activity (e.g.
Heckman et al. 2004; Hickox et al. 2009; Aird et al. 2012;
Bongiorno et al. 2012), down to accretion rates well below 0.01(e.g. Soria et al. 2006; Ho 2008; Gallo et al. 2010;
Miller et al. 2012). A different definition of an active blackhole will naturally result in a different active BHMF
(Goulding et al. 2010).
In this study we define an active black hole as a type 1
AGN with an Eddington ratio λ > 0.01, analogous to
Schulze & Wisotzki (2010). The lower Eddington ratio limit
is mainly a practical choice, corresponding to our observed
range. While we observe AGN below this limit in our sample, their low number and restricted black hole mass range
does not allow us to put strong constraints on the ERDF
below λ = 0.01. This differs somewhat from the more standard definition via an AGN luminosity limit but we think ismore physically motivated, in particular when studying thebivariate distribution of M• and λ. As shown in the rightpanel of Figure 6, for this definition of an active black hole
the SDSS sample suffers from severe incompleteness, while
the deeper VVDS and zCOSMOS samples are much lessaffected. However, even these deep surveys become increasingly incomplete towards the lowest M• and lowest Eddington ratios.
Therefore, we restrict the three samples to M• >107 M⊙ and λ > 0.01 for the determination of the BHMF
and ERDF, leaving 147,137 and 27238 AGN in the samples
from the VVDS, zCOSMOS and SDSS, respectively.
We further note that besides low accretion AGN our
census does not include obscured (type 2) AGN. Their issimply due to the fact that the virial method to estimate
M• for large statistical samples is limited to unobscured(broad line) AGN. In section 5 we will make an attempt
to account for the obscured population in the analysis.

4

THE ACTIVE BLACK HOLE MASSFUNCTION AND EDDINGTON RATIODISTRIBUTION FUNCTION

Early studies on the active BHMF focused on the 1/Vmax
method, directly applying the volume weights from the
AGN LF (Greene & Ho 2007; Vestergaard et al. 2008;
Vestergaard & Osmer 2009). However, given our definitionof an active black hole, this approach introduces unaccounted sample incompleteness at the low mass end of
the BHMF, since in this mass range not the full Eddington ratio range is sampled, due to the flux limit of thesurvey and the missing objects are not statistically accounted for (Kelly et al. 2009; Schulze & Wisotzki 2010;
Shen & Kelly 2012). Furthermore, this approach does notcorrect for the uncertainty in the virial black hole mass estimates (Kelly et al. 2009; Shen & Kelly 2012). These lim-

itations can be resolved by employing e.g. the maximum
likelihood approach described in the next section.
On the other hand, if interpreted with the necessary
caution, the determination of the binned BHMF and ERDF
via the 1/Vmax method can serve as a useful tool, to guidethe more refined parametric analysis presented in the following. We present and discuss our results for the BHMF
and ERDF as well as for the AGN LF, using the 1/Vmax
method, in Appendix B.
4.1

The maximum likelihood method

A proper determination of the intrinsic BHMF and the
ERDF requires a joint modeling of these two distributionfunctions to take the selection function fully into account.
In Schulze & Wisotzki (2010) we developed a maximum
likelihood approach to estimate the intrinsic BHMF and
ERDF simultaneously. We here present a modified framework, based on our previous work, but generalized to more
flexible models for the bivariate distribution function, morediverse selection functions and adapted to the combinationof several different surveys. A further improvement, compared to Schulze & Wisotzki (2010) is the inclusion of a correction for virial black hole mass uncertainties.
We aim at estimating the bivariate distribution function of black hole mass end Eddington ratio Ψ(M• , λ, z),i.e. Ψ(M• , λ, z)d log M• d log λ gives the space density of active black holes with black hole mass between log M• andlog M• + d log M• and Eddington ratio between log λ andlog λ+d log λ at the redshift z. The BHMF, ERDF and AGN
LF are all different projections of this bivariate distributionfunction.
The bivariate distribution function can be also expressed as a function of M• and Lbol , Ψ(M• , Lbol , z).
Marginalising this over M• gives the bolometric AGN LF,
i.e.
Z
Φ(Lbol , z) =Ψ(M• , Lbol , z)d log M• .
(2)
The survey selection modifies the intrinsic distributionfunction Ψ(M• , λ, z) to an observed distributionΨo (M• , λ, z) = Ω(M• , λ, z)Ψ(M• , λ, z) ,(3)

where Ω(M• , λ, z) is the selection function of the survey.
The maximum likelihood technique aims at minimizingQthe likelihood function S = −2 ln L, where L = N
i=1 pi isthe product of the individual likelihoods for the observed
objects (Marshall et al. 1983). This multivariate probability distribution p(M• , λ, z), i.e. the probability to observe
an AGN with black hole mass M• , Eddington ratio λ andredshift z, is given by the normalized observed bivariate distribution:pi (M• , λ, z) =1dV
Ωi (M• , λ, z) Ψ(M• , λ, z)
,
Ni
dz

(4)

where Ωi (M• , λ, z) is the selection function for the i-th object and
ZZZ
dV
Ni =
Ωi (M• , λ, z)Ψ(M• , λ, z)d log M• d log λ dz
dz
(5)is the normalization for the i-th object . If the selection
function is the same for all objects, also Ni will be the samec 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2
for them and corresponds to the expected number of objects
from the model. However, the method is also applicable for aselection function which varies from object to object. In thiscase Ni corresponds to the total number of objects, assuming
all objects would have the selection function Ωi (M• , λ, z).
Using the probability distribution given by Equation 4,
we minimize the function:
S = −2NX
i=1[ln Ψ(M• i , λi , zi ) − ln Ni ] ,(6)where N is the number of objects in the respective survey. This maximum likelihood method represents a forward
modeling approach. We have to assume a specific parametric model for the bivariate distribution function and we fit
the probability distribution predicted by this model in theM• − λ plane to the observations.
A drawback in this approach is that the normalizationΨ∗ of the bivariate distribution function cannot be determined directly in the fit. We therefore determined it in asecond step by integrating over the best fit model and scaling the predicted number of objects to the observed number,
Ψ∗ =

Nobs
.
Nmod

(7)

We will first determine the bivariate distribution function for each survey individually, but our main analysis
is based on the combined dataset from the three surveys.
We combine multiple surveys with their respective effective area and survey depth using the approach proposed by
Avni & Bahcall (1980). The overlap in area and luminosity
between SDSS, zCOSMOS and VVDS is marginal and there
are no objects in common to several samples. Therefore thesethree surveys can be treated as being independent.
The combined Pselection function follows asNjjΩe (M• , Lbol , z) =j=1 Ωe (M• , Lbol , z), where Nj isthe number of combined surveys and is adopted for the
combined dataset. Thus, it is assumed that each object
with its given luminosity, redshift and black hole mass
could potentially have been detected in each of the surveys,
if it otherwise satisfies the respective selection function ofthat survey.
Since the selection function, specifically the spectroscopic success rate (see Appendix A), is initially defined
as a function of magnitude m and redshift z, we need tomap Ω(m, z) → Ω(M• , λ, z) = Ω(M• , Lbol , z). We here follow the approach by Shen & Kelly (2012) and Kelly & Shen
(2013), by computing Lbol from the optical magnitude, as
discussed in Section 3.2. This ensures a direct and consistent mapping from magnitude to bolometric luminosity andtherefore from Ω(m, z) to Ω(M• , Lbol , z).
For the determination of the underlying distribution function it has to be noted that virial black holemasses are thought to have an error of ∼ 0.2 − 0.4 dex(Vestergaard & Peterson 2006; Park et al. 2012), as derived
from the scatter around reverberation mapping masses. We
here account for this uncertainty, an improvement comparedto Schulze & Wisotzki (2010)2 , by convolving Ψ(M• , λ, z)

9with the measurement uncertainty before comparing themodel distribution to the data:ZΨe (M•,e , λ e , z) = g(M•,e , λ e | M• , λ ) Ψ(M• , λ , z) d log M• d log λ .

(8)We assume a log-normal scatter distribution in M• and Lbol ,(µe − µ)2
(le − l)2
1
exp −
−
,
g(M•,e , λ e | M• , λ ) =222πσVM σbol
2σVM
2σbol
(9)
where µ = log M• is the black hole mass and l =log Lbol is the bolometric luminosity. For the uncertainty
in the bolometric correction σbol we assumed a scatter of
0.05 dex (Marconi et al. 2004), and for the uncertainty inthe virial mass σVM = 0.2 dex, consistent with previous
work (Kelly et al. 2010; Kelly & Shen 2013). The commonly
assumed measurement uncertainty associated to an individual viral mass estimate is ∼ 0.3 − 0.4 dex. However, for our
purposes we are only interested in the effect of uncertainties
on the apparent distribution of M• and λ of our statisticalsample, i.e. by which amount these distributions are actuallybroadened by using viral masses instead of true masses. Several studies suggest this value is actually smaller, ∼ 0.2 dex(Kollmeier et al. 2006; Fine et al. 2008; Steinhardt & Elvis
2010; Kelly et al. 2010). For a detailed discussion of thisdifference, its possible origin and consequences see e.g. Shen
(2013). Adopting a large value for σVM corresponds to a significant correction to the distribution functions, in particularwhen this value becomes comparable to the actual width ofeither the observed M• or λ distribution. In this sense, ourchoice of σVM = 0.2 dex for the scatter in the virial method
can be considered as a conservative value, since it corresponds to the smallest correction suggested by our current
knowledge of the virial mass uncertainties.

4.2

Parametric model

For the maximum likelihood fitting we have to assume a specific parametric model for the bivariate distribution function
Ψ(M• , λ , z). We assume a BHMF independent of Eddington
ratio, while we allow for a mass dependence in the ERDF.
This is motivated by the results of Kelly & Shen (2013), whofind a mass dependence in the ERDF when using virial blackhole mass estimates based on Mg ii, as we do in this study.
Ψ(M• , λ , z) = ρλ (λ , M• , z) ρ• (M• , z)ρz (z) ,(10)where ρλ (λ , M• , z) is an ERDF term, ρ• (M• , z) is the
BHMF term and ρ(z) is a redshift evolution term. The
BHMF is given by integration of the bivariate distributionfunction over λ,ZΦ• (M• , z) = Ψ(M• , λ , z) d log λ .
(11)Equivalently, the ERDF is given by integration of the bivariate distribution function over M• ,Z
Φλ (λ , z) = Ψ(M• , λ , z) d log M• .
(12)

2

The results for the local BHMF and ERDF presented in
Schulze & Wisotzki (2010) only marginally change if we correct for the uncertainty in virial black hole masses, assumingc 2014 RAS, MNRAS 000, 1–31

σVM = 0.2. The main effect is a slight decrease at the high accretion end of the ERDF (see Schulze 2011).

10

Schulze et al.
-3

log Φ •

-4
-5
-6
-7

VVDS

-8
7

zCOSMOS
8

9

10

7

log M • [M ⊙ ]

SDSS

8

9

10

7

log M • [M ⊙ ]

8

9

10

log M • [M ⊙ ]-3log Φλ

-4
-5
-6
-7
-8
-2

-1

0

-1

-2

log λ

0

log λ

-2

-1

0

log λ

Figure 7. Best fit maximum likelihood results for the BHMF (upper panels) and ERDF (lower panels) using different parametric models.
We show the different distribution functions for our 3 surveys independently, in the left panels for the VVDS, in the middle panels forzCOSMOS and in the right panels for SDSS. We used 4 different models for the bivariate distribution function: a modified Schechter
function for the BHMF with either a Schechter function (solid line) or a log normal function (dashed line) for the ERDF, or a double
power law for the BHMF, again either with a Schechter function (dashed dotted line) or a log normal function (dotted line) for the
ERDF. We also show the luminosity weighted Vmax results with the open circles (VVDS), squares (zCOSMOS) and triangles (SDSS).
Note that the maximum likelihood results are not a fit to these binned data, but correct for incompleteness in the binned estimate aswell as for uncertainty in the virial mass estimates. The former causes the deviation between our best fits and the bins at low M• / low
λ, while the latter is responsible for the difference at high M• / high λ.

We will start with the simplified assumption of a redshift independent BHMF and ERDF within the given redshift bin. In particular for a narrow redshift bin this is avalid assumption, while we would expect evolution over thefull range 1.1 < z < 2.1. We will investigate redshift evolution in the bivariate distribution function within the full
redshift range we are probing in more detail further below.
For the BHMF we use two different parametric models, that have been shown to be adequate to describe thelow-z BHMF (Schulze & Wisotzki 2010). First, we use amodified Schechter function (e.g. Aller & Richstone 2002;
Schulze & Wisotzki 2010), given by

α+1
β !

Ψ∗
M•
M•
ρ• (M• ) =
.
(13)
exp −
log10 e M•∗
M•∗
Second, we also utilize a double power law fit for the BHMF,
given byΨ∗ / log 10 e.
ρ• (M• ) =∗ −(α+1)(M• /M• )+ (M• /M•∗ )−(β+1)(14)

For the ERDF we use a Schechter function (Schechter 1976)
as our reference model, where we allow for a mass dependence in the break Eddington ratio:αλ +1
1λλρλ (λ, M• ) =.
exp −log10 e λ∗ (M• )λ∗ (M• )(15)
We parameterize the mass dependence in λ∗ by:log λ∗ (M• ) = log λ∗,0 + kλ (log M• − log M•,c ) ,(16)

where we fixed log M•,c = 8.0. We also tested a second orderpolynomial for the mass dependence in λ∗ , but found thebest fit second order term parameter to be consistent withzero.
Additionally, we tested a log-normal distribution, restricted at our lower Eddington ratio limit log λ = −2:1(log λ − log λ∗ (M• ))2√ρλ (λ, M• ) =
exp −
,
2σλ2
log10 e 2πσλ(17)where we also use equation (16) for the mass dependence inthe peak value λ∗ (M• ).
We parameterize the redshift evolution in the bivariate
distribution function of BHMF and ERDF by a density evolution term, where the normalization evolves in a log-linear
manner:ρz (z) = 10γ(z−zc ) .

(18)

We fixed zc = 1.6 to the central value of our redshift range.
The absolute normalisation of the BHMF and
ERDF is obtained by integrating Ψ(M• , λ , z) over λ orM• respectively. We chose integration intervals of −2 <log λ < 1 and 7 < log M• < 11. Our reference model has
7 free parameters determined in the fitting (M∗ , α, β, λ∗,0 ,αλ , kλ , γ) plus the normalisation Ψ∗ determined by Equation (7).
4.3

Results

We start by fitting the bivariate distribution function of
BHMF and ERDF to the three surveys independently, before combining the data sets. We fit the data with the maxc 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2
1.1 <z <1.3

1.3 <z <1.5

1.7 <z <1.9

1.9 <z <2.1

11

1.5 <z <1.7log Φ •

-4

-6

-8
7

-4

log Φ •

8910log M • [M ⊙ ]

-6

-8
7

8

9

10

7

log M • [M ⊙ ]

8

9

10

log M • [M ⊙ ]Figure 8. Active black hole mass function for the combined data set from VVDS, zCOSMOS and SDSS in narrow redshift bins over
the range 1.1 < z < 2.1, based on our maximum likelihood fitting approach. The magenta lines show the best fit results within the
narrow redshift bin only, while the black lines give the best fit model to the full range 1.1 < z < 2.1 evaluated at the central redshiftof the individual z-bin. The solid lines are for a model parameterization with a modified Schechter function BHMF and the dashedlines for a double power law BHMF. The magenta dotted line indicates the best fit double power law BHMF model for the redshift bin
1.5 < z < 1.7, shown for reference. The vertical magenta dotted line and magenta dashed line indicate the position of the break of the
BHMF at 1.5 < z < 1.7 and in the respective z bin.

1.1 <z <1.3

1.3 <z <1.5

1.7 <z <1.9

1.9 <z <2.1

1.5 <z <1.7log Φλ-4-6-8log Φλ

-4

-2

-1

0

log λ

-6

-8
-2

-1

log λ

0

-2

-1

0

log λ

Figure 9. Eddington ratio distribution function for the combined data set from VVDS, zCOSMOS and SDSS in narrow redshift bins
over the range 1.1 < z < 2.1, based on our maximum likelihood fitting approach. The different lines are the same as in Figure 8.

c 2014 RAS, MNRAS 000, 1–31

12Schulze et al.

z =1.3z =1.9

log λ

0.0

-0.5

-1.0

-1.5

-2.0
7

8

9

10

log M • [M ⊙ ]

7

8

9

10

log M • [M ⊙ ]Figure 10. Bivariate distribution function of black hole mass and Eddington ratio Ψ(M• , λ) for our best fit modified Schechter function
BHMF model at two redshifts. The contours show lines of constant space density, from 10−10 to 10−5 , separated by a factor of 10 each.

imum likelihood approach outlined in section 4.1, using theparametric models presented in section 4.2 over the full redshift range 1.1 < z < 2.1. At this stage, we ignore any
possible redshift evolution in Ψ(M• , λ , z), apart from thedensity evolution term (equation 18). We will postpone amore detailed investigation of evolution within our redshift
range to the analysis of the combined sample. We will comment on the results for the SDSS in narrower redshift bins
in section 4.4, when we directly compare our results to thework by Kelly & Shen (2013).
The results for the BHMF and ERDF for 4 different
parametric models for each survey are shown in Figure 7.
We use a modified Schechter function or a double power law
for the BHMF and a Schechter function or a (truncated) lognormal distribution for the ERDF. We also show the binned
luminosity weighted Vmax results for reference as open symbols. We emphasize that the maximum likelihood results are
not a fit to the binned distribution functions, but accountfor their inherent limitations, as discussed above. The binned
distribution functions are only shown for reference, keeping
their limitations in mind.
As expected, for the SDSS sample we find a large completeness correction compared to the luminosity weighted
1/Vmax result, due to the bright flux limit of the sample.
Contrary, for both VVDS and zCOSMOS the completenesscorrections are comparatively small. Since the VVDS sample
consists of the wide I < 22.5 mag and the deep I < 24 magfields, its completeness correction is the smallest, also compared to zCOSMOS, which only has a I < 22.5 mag field.
Note that at the high M• / high λ end the best fit falls off
steeper than the binned estimate, since the latter does notaccount for the virial SMBH mass uncertainties.
The 4 parametric model combinations give largely consistent results. Differences occur mainly in regions with poor
statistics and/or where large completeness corrections haveto be applied. This verifies that our results are robust to thespecific assumed parametric model as long as the completeness corrections are only moderate. This applies in particular over the range for which the luminosity weighted Vmax
results are consistent with the maximum likelihood model.
The SDSS dataset shows a larger dependence on the para-

metric model in the range with significant completeness corrections, i.e. M• < 109 M⊙ and λ < 0.1, demonstrating thelimitation of the SDSS alone to securely constrain the low
mass end of the BHMF and the low Eddington ratio end of
the ERDF at these redshifts.
In the following we will restrict our analysis to the
Schechter function ERDF and will not further consider the
log-normal ERDF model. The Schechter function model ismore flexible, since it allows for a turnover at low λ, but does
not enforce it, in contrast to the log-normal model. Furthermore it ensures an exponential decrease of the space density
towards the Eddington limit, while in the log-normal model
this decrease is in general flatter than exponential. We will
utilize two parametric models: always a Schechter function
for the ERDF and either a modified Schechter function or a
double power law for the BHMF.
These results on the individual surveys demonstrate
their respective strengths and weaknesses. The SDSS is powerful in constraining the high mass end of the BHMF andthe high λ end of the ERDF, while the uncertainties arelarge at lower M• and λ, due to the associated large completeness corrections required. On the other hand VVDS andzCOSMOS can reliably constrain the low M• and λ regime
BHMF and ERDF, while they fail at the high M• and highλ end, due to the poor statistics in this regime caused by the
relative small area of these surveys. Therefore, in the next
step we combine all three surveys to constrain the BHMF
and ERDF over a wide range in M• and λ. For the combined fit, we restrict the SDSS sample to the BH mass range
log M• > 8.5, to reduce the number of objects associatedwith the largest completeness corrections.
To test for redshift evolution in the bivariate distribution function Ψ(M• , λ) within our redshift range, we first
determined the best fit maximum likelihood model in narrowbins of redshift. The results for our two reference parametricmodels are shown in Figure 8 and Figure 9 by the magenta
solid line (modified Schechter function BHMF) and the magenta dashed line (double power law BHMF), respectively.
Both parametric models agree reasonably well for the ERDF
and for the BHMF at log M• > 8. At lower masses the twomodels show a stronger deviation, indicating the larger unc 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2log pλ (λ|M • )

0.2

0.4

0.6

-2

7.0
7.5
8.0
8.5
9.0
9.5
10.0
-1

0

log λ
Figure 11. Dependence of the ERDF on black hole mass. The
different lines show the conditional Eddington ratio distributionfunction at different black hole masses, from 107 to 1010 M⊙ .

certainty of the BHMF in this regime in the narrow z bins,
due to the relative small number of objects at log M• < 8in zCOSMOS and VVDS. At this point we are mainly interested in any detectable trends of evolution between thez bins. With the magenta dotted line we indicate the distribution function for the 1.5 < z < 1.7 bin for reference in all
panels. In addition, we mark the position of the break in the
distribution functions at 1.5 < z < 1.7 and the respective
redshift bin by the vertical dotted and dashed lines in all
panels. We identify two main trends of evolution between
the narrow z bins. First, the break of the BHMF seems toshift to higher M• with increasing redshift, and second, thebreak of the ERDF seems to move to higher λ with increasing redshift. To allow our parametric model to accommodate
such trends, we implement a log-linear redshift dependence
for the break of the BHMF and for the break in the ERDF:
log M∗ (z) = log M∗,0 + c• (z − zc ) ,(19)

log λ∗ (M• , z) = log λ∗,0 +kλ (log M• −log M•,c )+cλ (z −zc ) ,(20)
Furthermore, we allow for redshift evolution in the faint end
slope of the ERDF:
log αλ (z) = αλ,0 + cαλ (z − zc ) .

(21)

We fixed zc = 1.6 for all three equations. With thesethree additional free parameters, we fit our 10 parametermodels for the bivariate distribution function of M• and λ to
the data. The best fit results for the modified Schechter function BHMF (black solid line) and double power law BHMF
(black dashed line) models are shown in Figure 8 and Figure 9 and the best fit parameters are given in Table 1. We
find good agreement between the two parametric models. We
will use the modified Schechter function model as our reference model for further discussion, without implying any
clear preference for either of these parameterizations. We
comment if any differences arise from the two different model
implementations. In Figure 10 we directly show the bivariate distribution function of black hole mass and Eddington
ratio Ψ(M• , λ) for our best fit modified Schechter function
model at two redshifts.
Our
model
parameterization
allowsforac 2014 RAS, MNRAS 000, 1–31

13

M• dependence of the ERDF. Our best fit results for bothmodels supports such a dependence with kλ ≈ 0.10 − 0.15.
In Figure 11 we show the conditional probability distribution of λ at a given M• at several masses. At higher masses
the conditional Eddington ratio distribution is shifted
towards a higher break value. This is also directly visible
in the bivariate distribution function shown in Figure 10.
A similar behavior has been found in Kelly & Shen (2013)
only using the SDSS data, however only for black hole massestimates based on Mg ii, while they see no significant massdependence using Hβ or C iv. At this point it is therefore
unclear if this represents a real trend or is due to systematiceffects in the viral relations.
We also determined the bivariate distribution function
with this parameterization for the three individual samples.
The best fit models are given in Table 1. In Figure 12 wecompare our best fit model to the combined data set to theseresults for the three individual surveys. Since the BHMF andthe ERDF are only a marginalization over the bivariate distribution function Ψ(M• , λ), their respective normalization
depends on the shape of the other distribution function, i.e.
the normalization of the BHMF is determined by the integration over the ERDF and vice versa. This should bealways kept in mind when comparing different BHMFs or
ERDFs independently, since a difference in the ERDF between two fits will manifest in a different normalization of
the BHMF. To reduce the differences caused by this normalization condition, we normalize the distribution function over a range that is well restricted for the 3 surveys,
i.e. we use log λ > −1.5 as lower limit for the normalization of the BHMF and log M• > 8 for the normalization of
the ERDF. The BHMF for the VVDS and zCOSMOS are
in good agreement for M• < 109 M⊙ . On the other hand,the ERDFs for the two deep surveys differ. This might becaused by low number statistics and/or cosmic variance between the small area fields. Combining the two data sets,
as performed here, will alleviate both issues. Furthermore,
we again see the large uncertainties at the high mass end
and high Eddington ratio regime from the deep, small area
surveys. This regime is mainly constrained by the data fromthe SDSS QSO sample. Overall, the combined best fit bivariate distribution function Ψ(M• , λ) provides a good fit
to the three individual data sets and combines their individual strengths.
This can be furthermore demonstrated by a direct comparison to the data. Our maximum likelihood method fits amodel probability distribution in the M• − λ plane to theobserved data in that plane. In Figure 13 we compare theobserved distribution in the M• − λ plane for the 3 surveys,
VVDS (left panel, orange circles), zCOSMOS (middle panel,
blue squares) and SDSS (right panel, red contours), withthis distribution predicted from our best fit bivariate distribution function (black contours). We find a good agreement
between observations and the best fit model. Our bivariate
distribution function provides an excellent representation of
the observed mass and luminosity/ Eddington ratio distributions.
Furthermore, the bolometric AGN LF follows directly from the bivariate distribution function Ψ(M• , λ) bymarginalizing over M• via Equation 2. Thus, our model bivariate distribution function should match by design theobserved bolometric luminosity function. In Figure 14 we14Schulze et al.

Table 1. Best fit model parameters for the bivariate distribution function of black hole mass and Eddington ratio, i.e. the BHMF and
ERDF.
samplemodel

log(Φ∗ )

M•∗

α

β

λ∗,0

αλ

kλ

γ

c•

cλ

cαλ

SDSS

mod.S+S

-4.68

7.69

-1.30

0.50

-1.01

-1.42

0.146

0.00

−

−

−

SDSS

mod.S+S

-4.88

8.06

-1.19

0.57

-1.02

-1.09

0.074

0.00

−

−

−

SDSS

mod.S+S

-4.53

8.16

-1.21

0.59

-1.21

-0.87

0.155

0.00

−

−

−

SDSS

mod.S+S

-8.73

5.83

1.50

0.33

-1.17

-0.23

0.068

0.03

0.275

0.086

0.230

zCOSMOSmod.S+S

-4.73

7.28

-0.67

0.37

-0.84

-1.44

-0.107

0.12

0.292

0.148

0.048

VVDS

mod.S+S

-4.38

7.41

-0.46

0.47

-1.55

0.72

0.032

-0.22

-0.014-0.3250.521

combined

mod.S+S

-5.32

9.09

-1.50

0.96

-1.19

-0.29

0.099

0.05

0.270

0.100

0.094

combined

DPL+S

-5.64

9.23

-1.68

-4.56

-1.24

-0.36

0.148

0.31

0.198

0.114

-0.423

1.1 <z <2.1
-4log Φλlog Φ •

-4

-6

-8

-6

-8

7

8

9

10

-1

-2

0

log λ

log M • [M ⊙ ]Figure 12. Comparison of our best fit BHMF (left panel) and ERDF (right panel) between the combined dataset (black solid line) andthe individual surveys of VVDS (orange dashed line), zCOSMOS (blue solid line) and SDSS (red dashed dotted line). the normalizationof the BHMF is based on the ERDF with log λ > −1.5 and the normalization of the ERDF uses log M• > 8.

0.0

log λ

-0.5

-1.0

-1.5

-2.0
7

8

9

log M • [M ⊙ ]

10 7

8

9

log M • [M ⊙ ]

10 7

8

9

10

log M • [M ⊙ ]Figure 13. Comparison of the distribution of black hole masses and Eddington ratios in the M• − λ plane between the observational
data and the prediction derived from the best fit bivariate distribution function Ψ(M• , λ). The latter is always shown by the black
contours. Left panel: results for the VVDS (orange circles), Middle panel: results for zCOSMOS (blue squares), Right panel: results for
SDSS (red contours). The horizontal solid, dashed and dotted lines indicate Eddington ratios of 1, 0.1 and 0.01, respectively.

c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2Shen & Kelly (2012) for completeness and to give a direct
impression of the spread of possible solutions in the Bayesianframework used in both studies.

log ΦL [Mpc−3 ]

-4

-5

-6

-7

-8
45

46

47

log LBol
Figure 14. Comparison of the bolometric AGN luminosity function derived by marginalizing the best fit bivariate distributionfunction Ψ(M• , λ) to the one directly determined from the opticalmagnitudes. We show the model AGN LF based on the modified
Schechter BHMF model (black solid line) and the double powerlaw BHMF model (black dashed line). These are compared to thebinned estimates from the three surveys, VVDS (orange circles),zCOSMOS (blue squares) and SDSS (red triangles) and the best
fit PLE model to the combined data set (red dashed dotted line).

compare the best fit model bolometric luminosity function
with the observed bolometric AGN LF. Here we use theoptical AGN LF in the Mi (z = 2) band, presented in Appendix B1, and transform it to the bolometric AGN LF for
type 1 AGN, using the conversion outlined in section 3.2.
We find excellent agreement between the AGN LF predictedfrom Ψ(M• , λ) and the direct determination of the AGN LF,
again confirming the consistency of our maximum likelihood
fitting results.

4.4

15

Comparison with previous work

We here compare our results with previous work on the active BHMF and ERDF at z ∼ 1.5. Shen & Kelly (2012) and
Kelly & Shen (2013) both used the same SDSS DR7 dataset,
which we are also incorporating in this study, to determine
the active BHMF over the redshift range 0.3 < z < 5.
Both studies used a Bayesian framework to properly account for the selection function and determine the BHMF
and ERDF. Their work mainly differs in the detailed model
assumptions for their BHMF and ERDF. In particular,
Shen & Kelly (2012) focused more on the active BHMF andused a more simple and restrictive model for the ERDF,
while Kelly & Shen (2013) approached to consistently determine the joint bivariate distribution function of blackhole mass and Eddington ratio by employing a more flexible
parameterisation. Their results for the active BHMF andthe ERDF are shown in Figure 15 and Figure 16. The active BHMFs of both studies agree reasonably well, while the
ERDFs show larger deviations, due to the more restrictive
ansatz in the Shen & Kelly (2012) study, which leads to arather restrictive ERDF model. The work by Kelly & Shen
(2013) is closer to the approach we are following here in
determining the bivariate distribution function of M• andλ. Therefore, we will compare our results mainly to thework by Kelly & Shen (2013). We show the BHMF fromc 2014 RAS, MNRAS 000, 1–31

Our work uses the same data set and the same implementation of the selection function as the two studies above.
But while they used a Bayesian framework, we here employ
the maximum likelihood approach presented in this paperand first outlined in Schulze & Wisotzki (2010). Thus, weare able to directly compare the consistency of these two
methods. Our work also differs in the parametric model used
to fit the data. Kelly & Shen (2013) used a superpositionof multiple Gaussian functions for the bivariate distribution function, while our model parameterization is generally more restrictive, but includes less free parameters. We
already demonstrated above that this parameterization is
fully able to properly describe the data.
In addition to our main redshift range (1.1 < z < 2.1)
we also compute the BHMF and ERDF from the SDSS data
at lower redshifts, 0.5 < z < 1.1, where black hole masses
based on Mg ii are available, in three redshift bins. For thissubset we restricted the SDSS sample to log M• > 7.5 in
the fitting. This allows a comparison with previous studies
over a larger redshift range and will be further used in thediscussion.
In Figure 15 and Figure 16 we show our best fit results for the SDSS data only in each redshift bin by thetwo red lines for two different parametric models (modified
Schechter function or double power law BHMF). We generally find good agreement with the work by Kelly & Shen
(2013), while the difference with the Shen & Kelly (2012)
BHMF is larger. The main deviation occurs in the range
where the completeness in the SDSS sample in below 10 %,as indicated by the vertical dashed line. Thus, in this range
the SDSS suffers from significant uncertainty due to thelarge completeness correction, as already pointed out by
Kelly & Shen (2013). Overall the agreement between theBayesian method and the maximum likelihood approach isremarkable, demonstrating that both methods are equallyable to reliably determine the BHMF and ERDF. This is further shown by comparing to our best fit results for the combined data, including the deep VVDS and zCOSMOS surveys, shown by the black solid line. The BHMF and ERDF
by Kelly & Shen (2013) is generally consistent with the deep
data over the full range they are probing, log M• > 8 andlog λ > −1.5. At 1.3 < z < 1.9 their mass functions appear
to turn over at the lowest masses, which is probably an artifact of their parametric model of Gaussian superpositionsin the range of low completeness. Our deeper data do notshow such a turn over.
In the lowest redshift bin and in the highest redshift
bin we actually compare to results obtained from black holemasses estimated from other lines. While we always use
Mg ii, the BHMF and ERDF at z ≈ 0.6 from Kelly & Shen
(2013) is based on Hβ masses, while the ones at z ≈ 2.0 are
based on C iv masses. For these we see a larger disagreement, most prominent in the z ∼ 2.0 ERDF. This indicatesthat the currently employed virial mass calibrations do notprovide consistent BHMFs and ERDFs. Thus care has to betaken when comparing distribution functions obtained fromdifferent broad lines. In this study we therefore largely focus on using only masses from a single line, namely Mg ii. A

16

Schulze et al.
-4log Φ •

-5

0.5 <z <0.7Hβ0.7 <z <0.9

0.9 <z <1.1

1.1 <z <1.3

1.3 <z <1.5

1.5 <z <1.7

1.7 <z <1.9

1.9 <z <2.1
CIV

-6
-7
-8
-4

log Φ •

-5
-6
-7
-8
8.0

8.5

9.0

9.5

10.0 8.0

log M • [M ⊙ ]

8.5

9.0

9.5

10.0 8.0

log M • [M ⊙ ]

8.5

9.0

9.5

10.0 8.0

8.5

log M • [M ⊙ ]

9.0

9.5

10.0

log M • [M ⊙ ]Figure 15. Comparison of the active BHMF determined in this paper to previous studies based on the SDSS DR7 QSO sample. The
cyan solid line and gray shaded area shows the BHMF by Kelly & Shen (2013) and the yellow solid line gives the BHMF by Shen & Kelly(2012). The red lines show our best fit results using only the SDSS data in the given redshift bin (solid line for a modified Schechter
function BHMF and dashed line for double power law BHMF). The solid black line is our best fit result for the combined data from all
3 surveys and over the full redshift range (1.1 < z < 2.1). The vertical black dashed line indicates the mass at which the completeness in
the SDSS sample drops below 10 %, as given in Kelly & Shen (2013). Note that the BHMF by Kelly & Shen (2013) and Shen & Kelly(2012) in the lowest redshift bin is based on Hβ masses and in the highest redshift bin on C iv masses, while we use Mg ii masses in all
bins.

-4

log Φλ-50.5 <z <0.7Hβ0.7 <z <0.9

0.9 <z <1.1

1.1 <z <1.3

1.3 <z <1.5

1.5 <z <1.7

1.7 <z <1.9

1.9 <z <2.1
CIV

-6
-7
-8
-4

log Φλ

-5
-6
-7
-8
-1.5

-1.0

-0.5

log λ

0.0

0.5 -1.5

-1.0

-0.5

0.0

log λ

0.5 -1.5

-1.0

-0.5

log λ

0.0

0.5 -1.5

-1.0

-0.5

0.0

0.5

log λ

Figure 16. Comparison of the ERDF determined in this paper to the previous study by Kelly & Shen (2013), based on the SDSS
DR7 QSO sample. The lines are the same as in Figure 15. The vertical black dashed line indicates the Eddington ratio at which thecompleteness in the SDSS sample drops below 10 %.

more detailed investigation of this issue is beyond the scopeof this paper.
Recently, Nobuta et al. (2012) used the Subaru XMMNewton Deep Survey (SXDS) to determine the BHMF and
ERDF at z ∼ 1.4, employing the luminosity weighted
1/Vmax method and the maximum likelihood method from
Schulze & Wisotzki (2010). Their sample is X-ray selected
and extends significantly deeper than SDSS (4 × 10−15 erg
s−1 cm−2 in the soft band), over an area of ∼ 1.0 deg2 . Theirresults are shown in Figure 17 by the cyan circles for the
luminosity weighted 1/Vmax method and by the cyan lines
for the maximum likelihood method, where the dashed line

shows their double power law BHMF model and the dasheddotted line shows the Schechter function BHMF model. We
find their binned results in general in good agreement withour binned estimates from VVDS and zCOSMOS. At the
low mass end, the BHMF is in better agreement with zCOSMOS than with VVDS, probably due to a similar effective
flux limit of these two samples.
The maximum likelihood results from Nobuta et al.
(2012) show larger differences. The BHMF agrees well within
8.5 < log M• < 9.5. Similar to VVDS and zCOSMOS,SXDS is not able to reliably constrain the high mass endof the BHMF. At the low mass end their work predicts a
c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2

17

1.1 <z <2.1
-4log Φλlog Φ •

-4

-6

-8

-6

-8

7

8

9

10

log M • [M ⊙ ]

-2

-1

0

log λ

Figure 17. Comparison of our results for the 1 < z < 2 BHMF (left panel) and ERDF (right panel) with the work by Nobuta et al. (2012).
The cyan circles give the luminosity weighted 1/Vmax results from Nobuta et al. (2012), while the cyan dashed line and dashed dottedline give their best fit maximum likelihood results for their double power law BHMF and Schechter function BHMF model, respectively.
Our best fit results are shown by the black solid line. The ERDF is normalized by integrating the BHMF until M• = 108 M⊙ . We also
show our luminosity weighted binned estimates for the VVDS (open orange circles), zCOSMOS (open blue squares) and SDSS (open red
triangles).

turnover in the BHMF, which we do not confirm with ourdata. The reason for this discrepancy is not obvious. It might
be caused by cosmic variance between the small fields, differences in the model assumptions or in the target selection.
The SXDS sample is soft X-ray selected, while we use basically optically selected AGN samples.
In the right panel of Figure 17, we normalized the ERDF
by integrating the BHMF only until log M• = 8, to alleviate
the difference in normalization due to the BHMF difference
at lower M• . At log λ < −1, our work agrees well with
Nobuta et al. (2012), while they predict a much larger spacedensity at higher λ. This is probably caused by the larger
uncertainty at the high λ end, employing only the SXDS
survey, and the fact that they have not accounted for thevirial mass uncertainty in their fit, as we have done here.
Overall, our results are consistent with previous work,
but provide much more robust constrains on the bivariate
distribution function of M• and λ at 1 < z < 2 over a wide
range in M• and λ.

5
5.1

DISCUSSION
Cosmic evolution in black hole mass andaccretion rate

The cosmic evolution of the AGN luminosity function out toz ∼ 2 is by now well established (e.g. Richards et al. 2006a;
Bongiorno et al. 2007; Croom et al. 2009). The space density of the most luminous QSOs strongly increases towardsz ∼ 2 and peaks around this redshift. On the other hand, thespace density of lower luminosity AGN shows a weaker evolution. Their space density also increases, but peaks at lower
redshift. This behavior is known as ”AGN cosmic downsizing” (e.g. Hasinger et al. 2005). In the left panel of Figure 18 we illustrate this downsizing behavior of the AGN
LF by comparing our z ∼ 1.6 AGN LF with the AGN LF
at z = 0, as determined from the Hamburg/ESO Survey by
Schulze et al. (2009), shown by the dashed blue line. Therec 2014 RAS, MNRAS 000, 1–31

is not only a change in the normalization of the LF, but alsoa distinctive change in its shape, from being close to a single
power law at z = 0 to the presence of a prominent break atz ∼ 1.6.
With the additional information on the BHMF and
ERDF, we are able to disentangle the AGN downsizing behavior in its contribution due to black hole mass downsizing
and accretion rate downsizing. In the middle and left panelof Figure 18 we compare our z ∼ 1.6 BHMF and ERDF
with the distribution functions at z = 0, determined by
Schulze & Wisotzki (2010). Note however, that the comparison of our results and the z = 0 work by Schulze & Wisotzki(2010) may be affected by systematics in the black hole massestimates, since not the same broad lines (Mg ii vs. Hβ)have been used. Nevertheless, the general trends of evolution should be robust against these possible systematics.
The behavior of the active BHMF is similar to the LF.
We see a strong downsizing trend in the BHMF. At thehigh mass end, M• ≈ 109 M⊙ , there is a strong increase of
the space density from z = 0 to z ∼ 1.6, while at the low
mass end, M• ≈ 107 M⊙ , there is almost no change in thespace density. This mass downsizing has direct implications
for the active fraction (or duty cycle), as we will discuss inmore detail below.
In the Eddington ratio distribution function we see achange in the shape of the distribution function, largely corresponding to a flattening of the power law slope of the
Schechter function. There is a shift of the average accretion
rate towards higher values (see also Figure 20), i.e. if a supermassive black hole is active at z ∼ 1.6 it will on average
accrete at a higher rate than at z = 0.
To investigate these evolution trends further, in Figure 19 we directly show the cosmic evolution of the type-1
AGN space density in bins of bolometric luminosity (leftpanel), black hole mass (middle panel) and Eddington ratio(right panel) over the range 0 < z < 2.1. Here the lines at
1.1 < z < 2.1 give the prediction from our best fit model tothe full combined data set of VVDS, zCOSMOS and SDSS
over that redshift range. The squares give the z = 0 results

18

Schulze et al.
z =1.6z =0.0-4log Φ •

log Φ

-6

-4

log Φλ

-4

-6

-6

-8

-8

-8

-10
-20

-22

-24

-26

7

-28

8

Mi (z =2) [mag]

9

10

-1

-2

0

log λ

log M • [M ⊙ ]Figure 18. Comparison of the AGN distribution functions at z ∼ 1.6 with z = 0. Left panel: AGN luminosity function in Mi (z = 2),
Middle panel: active black hole mass function, Right panel: Eddington ratio distribution function. The black solid lines and blue dashedlines show the three distribution functions at z = 1.6 and z = 0, respectively.

log LBollog M •

log λ

-4

logΦ

-5

-6

-7

44.0
44.5
45.0
45.5
46.0
46.5
47.0

-8

-9
0.5

1.0

z

1.5

2.0

2.5

−2.0
−1.7
−1.4
−1.1
−0.8
−0.5
−0.2

7.0
7.5
8.0
8.5
9.0
9.5
10.0
0.5

1.0

z

1.5

2.0

2.5

0.5

1.0

z

1.5

2.0

2.5

Figure 19. Evolution of the type 1 AGN space density between z = 0 and z = 2.1 in bins of AGN bolometric luminosity (left panel),black hole mass (middle panel) and Eddington ratio (right panel). The lines at 1.1 < z < 2.1 give the prediction from our best fit
model over this redshift range. The triangles show the best fit to the SDSS data only in the redshift range 0.5 < z < 1.1 in z-bins of
0.2 dex. The squares give the z = 0 results from Schulze et al. (2009) and Schulze & Wisotzki (2010) and the dotted lines are simple
linear extrapolations between the z = 0 results and the our best fit model results at z = 1.1. A clear trend of AGN downsizing is present
in Lbol , M• and λ.

from the Hamburg/ESO survey (based on Hβ masses), while
the triangles show the results from our best fit to the SDSS
data only over bins in the redshift range 0.5 < z < 1.1 (using Mg ii for the SMBH mass estimation). The dotted lines
are simple linear extrapolations between z = 0 and z = 1.1
to help guide the eye. We normalized the ERDF by integrating the BHMF only to M• = 108 M⊙ , to allow a better
comparison with the SDSS only results at lower z.
Again, we confirm the well known behavior of AGN
downsizing in the AGN LF. The space density is strongly
decreasing between z = 2 and z = 0 at the bright end, withmuch less evolution at the faint end of the AGN LF. However, our restricted redshift range does not allow to probe
the turnover present at the bright end at z > 2. A very similar behavior is seen in black hole mass. At the high mass
end there is strong evolution in the space density, alreadywithin 1.1 < z < 2.1, but even more compared to z = 0,
with very little change in the space density at lower masses(M• < 109 M⊙ ). We also find evidence for evolution in Eddington ratio, though to a somewhat smaller extend. The

redshift dependence in the break of the Schechter function
ERDF in our best fit model implies a stronger decrease of thespace at the high λ end, visible in Figure 19. These trends
are also verified by the evolution of the mean M• and meanλ, as shown in Figure 20. The mean M• shows a strong decrease towards z = 0, while also the mean λ decreases fromz = 2 to z = 0.

We conclude that the downsizing in the AGN LF isdriven by both, the downsizing in the BHMF and evolution of the ERDF, while the BHMF evolution seems to be
the dominating term. At higher redshift more massive black
holes were active and these were on average accreting at
a higher rate. The AGN LF evolution is therefore a combination of simultaneous evolution in black hole mass andaccretion rate. Disentangling both processes can resolve degeneracies on the black hole growth history and accretion
rate evolution.
c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2-2-2

type 2type 1
total

7.7
log Φ •

7.6
7.5

-4

log Φλ-4log M •

19

-6

-6

7.4
-8

-8

7.3
7

8

9

log M • [M ⊙ ]

log λ

-1.4

10

-2

-1

0

log λ

Figure 21. Obsuration correction for the active BHMF (leftpanel) and the ERDF (right panel). The blue solid line shows
our best fit model for the active BHMF and ERDF of type 1
AGN. The black dashed line give the total BHMF/ERDF, inferred by using the luminosity dependent obscuration correction
from Merloni et al. (2014), while the red dashed-dotted line gives
the prediction for type 2 AGN.

-1.5
-1.6
-1.7
0.0

0.5

1.0

z

1.5

2.0

Figure 20. Evolution of the mean black hole mass (upper panel)and mean Eddington ratio (lower panel). The lower integration
limits for the computation of the mean are log M• = 7 andlog λ = −2. The blue solid line shows the respective mean value
from our best fit model to the SDSS+VVDS+zCOSMOS dataset, the blue triangles are for the SDSS sample only in the redshift range 0.5 < z < 1.1, and the blue square is the z = 0 resultsfrom the Hamburg/ESO Survey (Schulze & Wisotzki 2010). The
red open symbols and the red dashed line show the mean values for the obscuration corrected AGN population as discussed in
section 5.2.

5.2

Correction for obscured AGN

As noted above, our results for the active BHMF and ERDF
only refer to type 1 AGN, for which we can estimate virial
black hole masses. However, obscured AGN represent a significant fraction of the AGN population. This obscured fraction is luminosity dependent, with a decrease of the type 2
fraction with increasing luminosity (e.g. Ueda et al. 2003;
La Franca et al. 2005; Hao et al. 2005; Maiolino et al. 2007;
Hasinger 2008; Merloni et al. 2014). In the unified scheme,
this corresponds to a decrease of the opening angle of thetorus with increasing luminosity. But also AGN feedback
could be important, where stronger winds in more luminous AGN clean the line of sight quicker (Menci et al. 2008).
Thus, the total active BHMF and ERDF will differ in normalization and shape from the results for type-1 AGN which
we have discussed so far. However, for a better understanding of the physical processes at play and a comparison withtheoretical models, the bivariate distribution function for
the full AGN population would be desirable. We here apply
a simple correction for the obscured fraction to our type1 AGN results, to transform them to the expected BHMF
and ERDF of the total AGN population. While such a correction still suffers from significant uncertainties, it should
reveal the general trends and implications valid for the full
AGN population. Note that the total active BHMF still defines an AGN as having λ > 0.01, thus it does not include
low accretion rate or quiescent SMBHs.
Recently, Merloni et al. (2014) presented an obscured
c 2014 RAS, MNRAS 000, 1–31

fraction dependent on X-ray luminosity and consistent withno redshift evolution within 0 < z < 3.5 based on an Xray selected AGN sample from XMM-COSMOS. The massand Eddington ratio dependence of the obscured fraction isconsistent with a direct luminosity dependence. We use their
X-ray luminosity dependent obscuration correction
1
43.89 − log LX
Fobs (LX ) = 0.56 + arctan
,
(22)
π
0.46
to convert our bivariate distribution function Ψ(M• , λ) tothe distribution function for the full AGN population. To
reduce the large uncertainties in this correction at low luminosities, where the type 1 fraction is low, at log LX < 43
we fix the obscured fraction to the value at this luminosity, Fobs (1043 ) = 0.985 . We convert X-ray luminosities (210 keV) into bolometric luminosities via the bolometric correction from Marconi et al. (2004) and applied the correction factor above to Ψ(M• , λ).
The derived total active BHMF and ERDF are shownin Figure 21. This might still be an underestimate of the total AGN population, since we are missing heavily obscured,Compton thick AGN (CT AGN). The fraction of CT AGN
is currently highly uncertain, but expected to be between
10−30% (Gilli et al. 2007; Treister et al. 2009; Akylas et al.
2012; Alexander et al. 2013; Vignali et al. 2014). This would
further increase the space density of the total active BHMF
and ERDF by 0.05 − 0.15 dex. A possible M• and λ dependence of the CT AGN fraction is even more uncertain,
but recently Lanzuisi et al. (2014) found evidence that CT
AGN have on average smaller M• and higher λ with respect
to unobscured sources, suggesting that the inclusion of thispopulation can further change the shape of the active BHMF
and ERDF. However, given the large uncertainties involved
it is at the moment not warranted to apply such corrections
to our distribution functions for type 1 AGN.
At M• > 109 M⊙ the total active BHMF is dominated by type 1, with type 2 AGN becoming dominant
at lower masses. This leads to a significant steepening of
the low mass slope of the active BHMF. The inferred total
ERDF is also largely dominated by obscured AGN at thelow λ end. The accretion rate distribution function, as determined via the ERDF, links the physical quantity of themass of an SMBH to its instantaneous observable, the AGN

20

Schulze et al.

luminosity. It therefore contains information about the accretion process (Yu et al. 2005; Hopkins & Hernquist 2009;
Novak et al. 2011; Hickox et al. 2014). At z = 0 the accretion rate distributions obtained from type 1 and type 2 AGN
are consistent with each other (Schulze & Wisotzki 2010;
Yu et al. 2005). Our obscuration correction above implies
that at 1 < z < 2 both AGN populations show differencesin their accretion rate distribution. However, this indication needs to be verified by direct studies of the distributionfunctions for both AGN populations.
The flattened distribution of Eddington ratios foundhere for type 1 AGN and also inferred for the total AGN
population seems to be in contrast to recent studies of the
distribution of the specific accretion rate LX /M∗ . Aird et al.
(2012) determined the specific accretion rate distributionfunction of type 2 AGN at 0.2 < z < 1.0, defined as the
distribution function of the ratio of AGN luminosity overstellar mass. This quantity is effectively the convolution of
the ERDF with the M• − MBulge relation. They found a
power law distribution with slope −0.65 over the full redshift range, while the fit to the AGN LF implies a cut off
towards the Eddington limit, which is not directly probed
by their data (Aird et al. 2013). The different ranges in redshift, Eddington ratio (with λ < 0.1) and the different AGN
population probed makes a comparison to the results presented in this work difficult.
Using an X-ray selected COSMOS sample, including
both type 1 and type 2 AGN, Bongiorno et al. (2012) also
determined the distribution of the specific accretion rateover the redshift range 0.3 < z < 2.5. They also founda power law distribution with a steeper slope, γ ≈ 1, butalso find evidence for a sharp decrease of the space density
towards the Eddington limit, consistent with our results.
While a detailed comparison is again challenging, we notethat for the redshift range in common, their ERDF is only
well probed for λ > 0.1 (see their Figure 16). Over this λrange our ERDF, in particular the total ERDF shown in
Figure 21, is actually consistent with their data, while being
inconsistent with their power law
On the other hand, Lusso et al. (2012) did not findevidence for the steep power law Eddington ratio distribution functions reported by Aird et al. (2012) andBongiorno et al. (2012), based on a study of the XMMCOSMOS AGN sample. Their results rather suggest a lognormal distribution of Eddington ratios, both for type 1 andfor type 2 AGN. Future, more detailed studies of the BHMF
and the ERDF of both the type 1 and type 2 AGN population are required for a better understanding of the interrelation in the demographics of these two AGN populations.
5.3

The active fraction of black holes

The observed black hole mass downsizing has implications
for the evolution in the active fraction or duty cycle of
SMBHs. We here define the active fraction as the ratio ofactive black holes, according to our definition, to the total
SMBH population. Therefore ”active” only includes AGN
above an Eddington ratio of 1%. Our direct determinationof the BHMF is limited to type 1 AGN, missing obscured
AGN, while we will also adopt the obscuration correction
presented in the previous section to investigate the total,
obscuration corrected active fraction.

At z = 0, Schulze & Wisotzki (2010) found a decreaseof the type 1 AGN active fraction with increasing black holemass, implying that in the local universe the most massive
black holes are preferentially in a quiescent state, while lower
mass black holes are still actively growing (see blue dashedline in Figure 23). This suggests an enhanced black hole
growth episode of massive black holes at earlier times. With
our study we directly probe this epoch of enhanced growthof SMBHs.
The total BHMF, including quiescent black holes, isusually computed by convolving the stellar mass function(or galaxy luminosity function), converted to a spheroid
mass function, by the M• − MBulge relation. This total
BHMF is reasonably well established for the local universe (e.g. Marconi et al. 2004; Shankar et al. 2009). Theirextension to higher redshift has still large uncertainties(Li et al. 2011). First, the stellar mass function itself and thebulge-to-total ratios are less well established (Bundy et al.
2006; Fontana et al. 2006; Ilbert et al. 2010, 2013). Furthermore, the M• − MBulge relation is not well known atz ∼ 1.5. It has been suggested that there is redshift evolution in the M• − MBulge relation, with an increase of
the M• /MBulge ratio with redshift (e.g. Peng et al. 2006;
Decarli et al. 2010; Bennert et al. 2010). However, the observed apparent evolution is fully consistent with a nonevolution scenario, once sample selection effects are taken
into account (Schulze & Wisotzki 2011; Schulze & Wisotzki2014). We therefore assume no evolution in the M• − MBulge
relation out to z = 2.
We here use the stellar mass function from Ilbert et al.
(2013), without applying any conversion to a spheroid
mass function. Indeed, at z > 1 it is currently not clear
if SMBH mass correlates better with total stellar mass
or with spheroid mass. Taking into account the observational challenges of bulge-disk decomposition at high z,
current observations indicate that at these redshifts total stellar mass correlates at least as well with M• as
spheroid mass (Jahnke et al. 2009; Merloni et al. 2010;
Schramm & Silverman 2013). Furthermore, there is someevidence that even at lower redshift total stellar mass mightprovide an equally good correlation (Läsker et al. 2014;
Marleau et al. 2013).
The total black hole mass function is given by:(µ − α − βs)2exp −Φs (s, z) ds ,
2σ 2
(23)
where µ = log M• , s = log M∗ − 11, Φs (s, z) is the stellar
mass function at z and α, β and σ are the normalization,slope and intrinsic scatter of the M• − M∗ relation. We use
the M• −Mbulge relation from McConnell & Ma (2013), withα, β, σ = (8.46, 1.05, 0.34) to transform the galaxy stellar
mass function to the black hole mass function.
The resulting total BHMF at z = 1.6 is shown by theblack dashed line in the upper panels of Figure 22, withthe uncertainty associated with the underlying stellar massfunction shown by the grey area. The black solid line in theupper panels of Figure 22 gives the active BHMF based on
the modified Schechter BHMF. The left panel shows the results for type 1 AGN and the right panel gives our estimate
for the total AGN population, applying the obscuration correction from Merloni et al. (2014). The lower panels presentΦ• (µ, z) = √

1
2πσ

Z

c 2014 RAS, MNRAS 000, 1–31

log (active fraction)log Φ • [Mpc−3 logM •−1 ]The cosmic growth of the active black hole population at 1 < z < 2-2type 121

obs. corrected

-3
-4
-5
-6
-7

-8
0.0

-0.5
-1.0
-1.5
-2.0
7

8

10 7

9

8

log M • [M ⊙ ]910log M • [M ⊙ ]Figure 22. Upper panel: Comparison of the active black hole mass function at z ∼ 1.5 (black solid line for the modified Schechter BHMF
model and dotted line for the double power law BHMF model) with the total black hole mass function at this redshift (black dashed linewith 1σ confidence region given by the gray area). Lower panel: Active fraction of black holes as a function of black hole mass (dashedline for the modified Schechter BHMF and dotted line for the double power law BHMF). The gray area gives the uncertainty, basedon the stellar mass function uncertainty and the two active BHMF models. The blue dashed dotted lines gives the BHMF and active
fraction for star forming galaxies. The active fraction for type 1 AGN at z ∼ 1.5 is consistent with being black hole mass independent.
The left panels present the results for type 1 AGN, while the right panels use an obscuration correction to derive the active fraction of
the total AGN population.

log (active fraction)

0

-1

-2

z =0.1z =0.6z =0.8z =1.0z =1.2z =1.4z =1.6z =1.8z =2.0-3

type 1-4789log M • [M ⊙ ]7obs. corrected89log M • [M ⊙ ]Figure 23. Redshift evolution of the active fraction for the type 1 AGN population (left panel) and the total AGN population (rightpanel). The lines within 1.1 < z < 2.1 show the results for our best fit model to the combined VVDS, zCOSMOS and SDSS dataset.
The lines within 0.5 < z < 1.1 are based on a maximum likelihood fit to the SDSS data only in three redshift bins. The blue dashed line
gives the active fraction at z ∼ 0, based on the BHMF from Schulze & Wisotzki (2010). A clear AGN downsizing trend in black holemass is present in the active fraction, where we see the shutoff of AGN activity for the high M• population between z = 2 and z = 0.

the active fraction for these two populations by the black
dashed line. The active fraction for the type 1 AGN population is almost constant, at ∼ 3 %, over a wide mass range. At
higher masses the active fraction becomes more uncertain,
due to the uncertainties in the active BHMF, the galaxy stellar mass function and the intrinsic scatter in the M• − M∗
c 2014 RAS, MNRAS 000, 1–31

relation at this redshift. There is tentative evidence for an
increase of the active fraction around log M• ∼ 9.0, which
we will investigate in more detail below. The implied strong
decrease of the active fraction towards log M• ∼ 10.0 is not
significant, due to the statistical uncertainties at these large
masses. Already using the best fit double power law BHMF

22Schulze et al.

would remove this turnover, as shown by the dotted blackline.
The lower right panel shows the active fraction for the
total AGN population, showing a different mass dependence
than the type 1 AGN active fraction. The active fraction
decreases until log M• ∼ 8.5 and flattens out at higher
masses. Our obscuration corrected active fraction implies
that ∼ 30% of all galaxies at M• ∼ 107 are active, i.e accrete above 1% of the Eddington limit, while this fraction
is still ∼ 10 − 20% at M• > 109 , the latter being consistent
with observations of X-ray AGN in deep survey fields (e.g.
Brusa et al. 2009; Xue et al. 2010; Bongiorno et al. 2012).
The close connection between star formation and AGN
activity (e.g. Silverman et al. 2009; Rosario et al. 2012)
might imply that the star forming galaxy population mightprovide the main reservoir to supply the AGN population,
before AGN feedback is quenching star formation. To test
the implications of this hypothesis, we compare the shapeof the mass function of star forming galaxies with the AGN
population. We use the stellar mass function of star forming galaxies, based on the work by Ilbert et al. (2013), tocompute the respective BHMF, shown by the dashed dottedblue line in Figure 22. We caution that there are significant
uncertainties in the assumptions made for the estimation of
the black hole mass function for this subpopulation, while
the main trends should be relative robust. Since at theseredshift the star forming galaxy population is dominating
the total stellar mass function, the difference between thestar forming and total galaxy population is small, and hence
also their differences in comparison to the AGN population.
Thus, we cannot distinguish between the star forming andtotal galaxy population as the main source to supply the
AGN population.
In Figure 23 we study the redshift evolution of the active
fraction, both for the type 1 AGN population and for the
obscuration corrected AGN population. The individual lines
show the active fraction at different redshifts. To extend theredshift coverage of our study, we augment our results byour best fit BHMF to the SDSS DR7 QSO sample in three
redshift bins within 0.5 < z < 1.1, limited in the M• range
to M• > 108 M⊙ , as presented above (Note that these suffer from larger uncertainties already below log M• ≈ 8.5,
partly being responsible for the discontinuity seen betweenz = 1.2 and z = 1.0). Furthermore, we show by the blue
dashed line the active fraction at z ∼ 0, derived from thelocal active BHMF by Schulze & Wisotzki (2010). We computed the z = 0 active fraction using the McConnell & Ma(2013) M• − Mbulge relation and an estimate of the spheroid
mass function, based on the local early type and late typegalaxy stellar mass function from Bell et al. (2003), assuming B/T = 0.3 for the late type galaxy population (see also
Schulze & Wisotzki 2011).3
Figure 23 again reveals a strong mass dependence in theblack hole growth history. While there is only moderate evolution of the active fraction at the low mass end, we witness
the shutoff of black hole growth at the high mass end between z = 2 and z = 0. Focusing only on our 1.1 < z < 2.1

3

If we would instead use the total stellar mass function fromBell et al. (2003), the z = 0 active fraction would shift down by∼ 0.15 dex.

sample, we find that at the low redshift edge the active fraction is basically constant with black hole mass, i.e. black
holes are actively accreting independent of black hole mass
with an active fraction of ∼ 2 − 3%. With increasing redshift, we see an increase of an upturn of the active fraction
for the most massive black holes, M• > 109 M⊙ . These become more active than the low mass population. At the peak
epoch of AGN activity, around z ∼ 2, about 20−30 % of themost massive black hole are in an active state (i.e. accrete at
λ > 0.01). Adding the information from the SDSS at z < 1.1
and the z = 0 results, this trend continues, with the most
massive black holes shutting off towards z = 0 while there
is little evolution at low M• , leading to an active fraction
decreasing with increasing M• at z < 1. The trend is present
both in the type 1 AGN population and in the total AGN
population. We note that the strength of the upturn at highM• towards z = 2 would be reduced, but is still present, if wewould use the M• − Mbulge relation from Kormendy & Ho
(2013). Contrary, it would be much more pronounced for the
relation by Häring & Rix (2004). Furthermore, assuming alarger intrinsics scatter in the M• − Mbulge relation would
also flatten this upturn. The general redshift evolution trend
however is robust against these uncertainties.
5.4

Comparison to model predictions

Current theoretical models of galaxy evolution and blackhole growth are reasonably well able to reproduce theobserved AGN LF and thus, the cosmic downsizing ofAGN. These models range from large cosmological hydrodynamic simulations (Degraf et al. 2010; Khandai et al. 2014;
Hirschmann et al. 2014; Sijacki et al. 2014; Bachmann et al.
2014, Rosas-Guevara et al., in preparation), semianalytical models (Cattaneo et al. 2005; Marulli et al. 2008;
Bonoli et al. 2009; Fanidakis et al. 2012; Hirschmann et al.
2012; Menci et al. 2014; Neistein & Netzer 2014) to more
phenomenological models (Wyithe & Loeb 2003; Shen
2009; Conroy & White 2013; Veale et al. 2014; Caplar et al.
2014). However, since these models differ in many details,
in particular the trigger mechanisms for AGN activity, no
consistent picture on the physics of black hole growth can beinferred from these comparisons alone. In fact, the AGN LF
alone is degenerate to physically distinct black hole growthmodels (e.g. Veale et al. 2014). The addition of the observed
active BHMF and ERDF provides important independent
constraints for the theoretical models of galaxy and AGN
evolution. Due to varying assumptions for BH growth, thesemodels predict or use differently evolving Eddington ratio
distributions.
While a comprehensive comparison with theoretical predictions in the literature is beyond the scope of this work,
we here want to illustrate the additional constraints gained
from SMBH mass and Eddington ratio distribution for theoretical models. We compare our observations to the study
by Hirschmann et al. (2014). They analysed a subset ofthe cosmological, hydrodynamical simulation Magneticum
Pathfinder (Dolag at al., in preparation), based on an improved version of the SPH code GADGET3 (Springel 2005),
with a comoving box size of (500 Mpc)3 . Their model successfully reproduces the observed AGN LF at z < 3 andexplains their downsizing behavior by the evolution of thegas density in the vicinity of the SMBH. This gas reserc 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2
z =0.0log Φ •

-4
-5
-6
-7

z =0.5log Φ •

-4
-5
-6
-7

z =1.0log Φ •

-4
-5
-6
-7

z =1.5log Φ •

-4
-5
-6
-7

z =2.0log Φ •

-4
-5
-6
-7
8.0

8.5

9.0

9.5

10.0 -2.0

log M • [M ⊙ ]

-1.5

-1.0

-0.5

0.0

log λ

Figure 24. Comparison of our observed active BHMF and
ERDF, corrected for obscuration (black dashed line) with the results of the hydrodynamical simulation from Hirschmann et al.
(2014), shown by the red solid lines, for the redshift range
0 < z < 2. We also show our observed distribution functions
without obscuration correction for reference (black dashed dotted lines).

voir for accretion gets depleted over cosmic time, due tostar formation and AGN feedback. This decrease is stronger
for more massive black holes, i.e. in more massive galaxies,
which therefore shut off their high accretion phase at earlier
cosmic times than less massive black holes.
In Figure 24 we show the predictions of the simulation
studied in Hirschmann et al. (2014) for the active BHMF
and the ERDF by red solid lines. We here define an active
SMBH consistent with our work as accreting above 1% Eddington. Below log M• ≈ 8 their simulations under predict
the total SMBH space density (see their Figure 5) due totheir relatively low resolution in the large volume box, thus
we also expect an under prediction for the active BHMF forthat mass range. Therefore we restrict our comparison tolog M• > 8. The shown ERDF is also normalized to thislower mass limit. We compare the simulation with our distribution function, corrected for obscuration (black dashedline), while we also show our type 1 AGN result for reference (black dashed dotted line). At z > 1 and log M• < 9.5
the simulations of the active BHMF are in reasonably good
agreement with the observations, in particular at z = 2.
However, they do not match our observations at lower redshift or higher masses. At z = 1 and z = 1.5 the simulations predict a higher space density at the massive end
log M• > 9.5. This is probably related to the overestimation
also seen in the stellar mass function at these redshifts (seetheir Figure 4), likely a consequence of their implementac 2014 RAS, MNRAS 000, 1–31

23

tion of AGN feedback in the radio-mode, which is relatively
inefficient in locking the baryons in the hot gas halo andthus, in suppressing star formation and gas accretion onto
the central BH.
For z < 1 the simulations predict a strong decrease of
the AGN space density, especially at the high mass end, in
qualitative agreement with our results. However, quantitatively this shutoff of AGN activity is much stronger than weobserve, leading to a significant underestimation of the AGN
space density. The reason is likely again the relatively inefficient radio-mode AGN feedback implementation. This let
massive black holes grow too much at higher redshift andthus too much gas is consumed to allow active black hole
growth at low redshift.
The right panel in Figure 24 compares our observed
ERDF with the simulation result. In general, the shape of
the ERDF is in reasonable agreement with our observation,
especially at z = 1. They predict a wide distribution of Eddington ratios with a decrease towards the Eddington limit.
In detail there are differences at high and low λ. The simulations show a less steep decrease of the space density towards the Eddington limit. At the same time they tend tounderestimate the space density at lower Eddington ratio,
λ < 0.1 %, compared to our observations. At z < 1 also thenormalization is much lower, but this is a direct consequenceof the underestimation of the active BHMF.
Overall, the agreement of the hydrodynamical simulations with the observations of the active BHMF and ERDF
over certain ranges in parameters space is already encouraging. However, clear differences remain which can guide thefurther improvements of the theoretical models.
While numerical simulations potentially contain themost complete physical scenario, their interpretation canbe complex. Semi-empirical or phenomenological models complement these by allowing to easily test the effect of individual model parameters. We here compare
our results to model predictions based on the semiempirical model from Shankar et al. (2013). They model
black hole growth via a standard continuity equation approach (e.g. Small & Blandford 1992; Yu & Tremaine 2002;
Marconi et al. 2004; Merloni 2004; Shankar et al. 2009), using observational constraints of the AGN LF and the local
total BHMF together with assumptions on the ERDF. TheirERDF model goes beyond the single Eddington ratio often
adopted in these kind of studies and assumes a Gaussian
distribution, with potential dependence on redshift, M• or
the radiative efficiency ǫ. This Gaussian distribution differsfrom the Schechter function ERDF inferred in this work,
so differences in the details between these models and our
results would be anticipated.
In Figure 25 we compare different models from
Shankar et al. (2013), given by the red and cyan lines, with
our active BHMF with obscuration correction (black dashedline) at two characteristic redshifts. The type 1 BHMF isagain shown only for reference. We compare three different
models for the ERDF: (i) a Gaussian in log λ with dispersionof 0.3 dex and a mean at log λc = −0.6 [G], (ii) a Gaussianwith the same dispersion, but a mean decreasing with redshift as (1 + z)2.2 [G(z)], (iii) a z dependent Gaussian withan additional dependence on M• , with decreasing mean λwith increasing mass [G(z, M• )].
All models predict a lower space density at low masses

24

Schulze et al.

z =0.5-3log Φ •

-4
-5
-6
-7

z =1.6-3log Φ •

-4
-5
-6
-7
7

8

9

log M • [M ⊙ ]Figure 25. Comparison of our active BHMF with the semiempirical models of black hole growth from Shankar et al. (2013)at two characteristic redshifts, z = 0.5 and z = 1.6. The blacklines show our results for the obscuration corrected AGN population (dashed) and for type 1 AGN only (dashed dotted). The redand cyan lines present different models from Shankar et al. (2013)
always assuming a log-normal ERDF. They tested models with aconstant ERDF (red solid line), but also included a z dependence(red dashed line) and an additional M• dependence (cyan dasheddotted line).

compared to our obscuration corrected BHMF. This is likely
a consequence of the lower space density predicted by themodels for the total BHMF, compared to the estimate derived above from the stellar mass function. It should be
noted that the latter estimate at these redshifts is also
highly uncertain at the low mass end, so it is not clear
if this indicates a shortcoming in the models. Overall, we
find the best agreement for the redshift dependent Gaussian G(z) model at both redshifts. The G model provides afair match at z = 1.6, but underestimates the space density
at low redshift, which is consistent with the conclusions in
Shankar et al. (2013), based on a comparison to z ∼ 0. Allother models overestimate the space density at the high mass
end at z = 0.5 and partially also at z = 1.6. This apparent
disagreement at high masses is mainly a consequence of thepower law behavior of the BHMF model in Shankar et al.
(2013), which reflects the shape of the bolometric AGN LF,
imposed in the their fitting. At z = 1.6 the G(z) model
provides clearly a better agreement with our observations
than the M• dependent model, G(z, M• ), while this difference is less strong at z = 0.5. The M• dependent model was
motivated by some observational evidence for rising active
fractions of AGN with M• at a given z, a trend we do notconfirm with our results. Thus the G(z) model naturally
performs better in comparison to our work.

In total, we find that the relative simple redshift dependent Gaussian model presented in Shankar et al. (2013)
already provides a reasonable match to our observations. We
note that this agreement with their models which assume alog-normal ERDF does not directly support such a distribution function. The comparison is mostly valid around thebreak and above of the BHMF. This regime is mainly regulated by the highest λ sources, where a log-normal modeland a Schecher function model are basically consistent. The
disagreement we see at lower masses could be a consequenceof the mismatch in the ERDF.
This comparison rather highlights the potential of combining our observational constraints with such models forconstraining the black hole growth history. The standard
continuity equation approach uses only the AGN LF as observational constraint to study the black hole growth history.
While this provides strong constraints on integrated quantities like the black hole density (i.e. the Soltan argument),
the AGN LF is degenerate in constraining the details of
the accretion history (see e.g. Veale et al. 2014). Therefore
Shankar et al. (2013) already included additional observational constraints on active fractions (i.e. the active BHMF)
to derive more detailed models of black hole growth. Directly including the observational determination of the active BHMF and the ERDF presented in this work into thesekind of empirical models has the potential to refine thesestudies, in particular in regard to the detailed accretion andtriggering history.

6

SUMMARY AND CONCLUSIONSWe here present a census of the broad line active SMBH
population at redshifts 1.1 < z < 2.1. We employ 3 different
type 1 AGN surveys, together covering a wide range in AGN
luminosity. These are the VVDS epoch-2 AGN sample, thezCOSMOS 20k AGN sample and the well defined subset of
the SDSS DR7 QSO sample. All 3 samples have well defined selection functions, allowing the determination of the
AGN distribution functions. Besides the AGN luminosity
function, we here focus on the more fundamental distribution function, namely the bivariate distribution function of
black hole mass and accretion rate. The marginalization ofthis bivariate distribution function to one dimension gives
the active black hole mass function, the Eddington ratio distribution function or alternatively the AGN LF.
We obtain virial black hole mass estimates and Eddington ratios for the majority of objects either from the literature or from our own spectral fitting. Due to the bright flux
limit of the SDSS, this survey is limited to more luminous
AGN and to high black hole masses (M• & 108 M⊙ ) and
Eddington ratios (λ & 0.03). The deeper surveys from the
VVDS and zCOSMOS extend this range down by ∼ 2 dex in
AGN luminosity, ∼ 1 dex in black hole mass and ∼ 0.5 dex
in Eddington ratio. They fill in the space in the massluminosity plane not occupied by the SDSS due to its flux
limit. On the other hand, due to the large survey area SDSS
provides an excellent coverage of the bright end of the LF
and thus of the high mass, high Eddington ratio regime,
which is poorly constrained by VVDS and zCOSMOS due totheir limited volume. Therefore, it is essential to use a multilayer approach and combine the information from bright,
c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2large area surveys with deep, small area surveys, like wehave performed it in this work.
We employ a maximum likelihood fitting technique,based on our previous work (Schulze & Wisotzki 2010), but
generalized to be more flexible and applicable to our specific surveys. The approach consists of fitting a parametricmodel for the bivariate distribution function of black holemass and Eddington ratio and its redshift evolution to theobservations. It accounts for incompleteness due to the sample selection and also for the uncertainty in the virial blackhole mass estimates that broadens the observed mass andEddington ratio distributions. This way we determine theactive BHMF and ERDF. We stress that our definition ofan active black hole refers to type 1 (broad line) AGN accreting above an Eddington ratio of log λ = −2. Our best
fit model accounts for a black hole mass dependence in the
ERDF and allows for a linear redshift dependence in the normalization, the break mass in the BHMF and in the breakand low λ slope of the ERDF. Our combined bivariate distribution function of M• and λ, Ψ(M• , λ), is consistent withthe best fit result for the individual surveys and with previous work, but significantly improves the reliability and the
dynamical range of each individual study. The thereby predicted AGN LF is also in excellent agreement with theirdirect determination.
We investigate the cosmic evolution of the AGN population as a function of AGN luminosity, black hole massand Eddington ratio between z = 0 and z = 2. To extend
our results into the redshift range 0.5 < z < 1.1, we additionally use the SDSS QSO sample to consistently determine the bivariate distribution function within this redshift
range, down to M• = 108 M⊙ . The bivariate AGN distribution function Ψ(M• , λ) allows to disentangle the well known
downsizing behavior we see in the AGN LF into its contribution due to black hole mass downsizing and Eddington ratio evolution. The downsizing is well represented in
the active BHMF, with the most massive SMBHs experiencing the strongest decrease of their space density, while
the space density of low mass SMBHs is approximately constant throughout 0 < z < 2. The Eddington ratio distribution function shows a flattening at the low-λ end, with analmost constant distribution at log λ < −1, which is in contrast to the shape at z = 0. We find evidence for an increaseof the average Eddington ratio with redshift.
Further evidence for AGN black hole mass downsizing
comes from the cosmic evolution of the M• dependence of
the active fraction or duty cycle. Overall, we find an active fraction of type 1 AGN almost constant with black holemass, implying a phase of active accretion onto black holes
throughout the galaxy population. Furthermore, we see evidence for an upturn of AGN activity at the highest masses,
M• > 109 M⊙ between z = 1 and z = 2. At z ∼ 2, during
the peak period of AGN activity, the most massive black
holes reach a peak in their activity. At lower redshifts more
and more massive black holes shutoff their mass accretionand become inactive, while lower mass black holes continue
to accrete actively. At the same time the average accretion
rate shifts to lower values, as accretion at low rates becomes
more frequent.
While our results here strictly only apply to the type 1
AGN population, we make an attempt to account for obscured AGN. Including a correction for obscured sources
c 2014 RAS, MNRAS 000, 1–31

25

changes the AGN space density, while our main conclusions
remain unchanged.
Finally, we compare our results to model predictions.
First, to the large size, hydrodynamical simulation by
Hirschmann et al. (2014). We find a reasonably good agreement at z > 1 and log M• < 9.5, but also identify distinct differences between our observations and the simulation results. Robust observational results of the ERDF are
expected to provide important constraints on the implementation of the accretion process and their evolution for numerical simulations and semi-analytic models. Second, we
compared with the semi-empirical models by Shankar et al.
(2013), finding a reasonable agreement with their relative
simple model using a redshift dependent Gaussian ERDF to
trace black hole growth through cosmic time.
Our work strongly supports the picture of black holemass downsizing or ”anti-hierarchical” black hole growth,suggested by previous studies (Merloni 2004; Heckman et al.
2004; Shankar et al. 2009; Vestergaard & Osmer 2009;
Shen & Kelly 2012; Nobuta et al. 2012). The downsizing in
the AGN LF is driven by mass downsizing and by a changeof the accretion rate distribution function with time, suggesting a complex process driving the cosmic evolution of
black hole growth.

ACKNOWLEDGMENTS
We thank Knud Jahnke for fruitful discussions on this workand Micol Bolzonella for providing the zCOSMOS target
sampling rates. AS acknowledges support by the China Postdoctoral Scientific Foundation, grant 2012M510256. This
work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan.
ABs work is supported by the INAF-Fellowship Program
and through PRIN-INAF 2011 ”BH growth and AGN feedback through the cosmic time”. IG acknowledges support
from FONDECYT through grant 11110501.
Based on observations made with ESO Telescopes atthe La Silla or Paranal Observatories under programme ID
175.A-0839. We gratefully acknowledge the contribution of
the entire VVDS and COSMOS collaborations.
Funding for the SDSS and SDSS-II was provided by the
Alfred P. Sloan Foundation, the Participating Institutions,
the National Science Foundation, the U.S. Department ofEnergy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England.
The SDSS Web site is http://www.sdss.org/.

REFERENCES
Aird, J., Nandra, K., Laird, E. S., et al. 2010, MNRAS,
401, 2531
Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746,
90
Aird, J., Coil, A. L., Moustakas, J., et al. 2013, ApJ, 775,
41
Akylas, A., Georgakakis, A., Georgantopoulos, I., Brightman, M., & Nandra, K. 2012, A&A, 546, AA98

26

Schulze et al.

Alexander, D. M., Stern, D., Del Moro, A., et al. 2013,
ApJ, 773, 125Aller, M. C. & Richstone, D. 2002, AJ, 124, 3035
Anglés-Alcázar, D., Özel, F., & Davé, R. 2013, ApJ, 770, 5
Assef, R. J., Kochanek, C. S., Ashby, M. L. N., et al. 2011,
ApJ, 728, 56
Avni, Y., & Bahcall, J. N. 1980, ApJ, 235, 694
Babić, A., Miller, L., Jarvis, M. J., et al. 2007, A&A, 474,
755
Bachmann, L. K., Dolag, K., Hirschmann, M., Almudena
Prieto, M., & Remus, R.-S. 2014, arXiv:1409.3221
Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D.
2003, ApJS, 149, 289
Bennert, V. N., Treu, T., Woo, J., et al. 2010, ApJ, 708,
1507
Bentz, M. C., Peterson, B. M., Netzer, H., Pogge, R. W.,& Vestergaard, M. 2009, ApJ, 697, 160
Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, ApJ,
767, 149
Blandford, R. D. & McKee, C. F. 1982, ApJ, 255, 419
Bongiorno, A., Zamorani, G., Gavignaud, I., et al. 2007,
A&A, 472, 443
Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, 3103
Bonoli, S., Marulli, F., Springel, V., et al. 2009, MNRAS,
396, 423
Booth, C. M. & Schaye, J. 2011, MNRAS, 413, 1158
Boyle, B. J., & Terlevich, R. J. 1998, MNRAS, 293, L49
Boyle, B. J., Shanks, T., Croom, S. M., et al. 2000, MNRAS, 317, 1014
Brusa, M., Fiore, F., Santini, P., et al. 2009, A&A, 507,
1277
Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716,
348
Bundy, K., Ellis, R. S., Conselice, C. J., et al. 2006, ApJ,
651, 120
Caplar, N., Lilly, S., & Trakhtenbrot, B. 2014,arXiv:1411.3719
Cattaneo, A., Blaizot, J., Devriendt, J., & Guiderdoni, B.
2005, MNRAS, 364, 407
Conroy, C., & White, M. 2013, ApJ, 762, 70
Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996,
AJ, 112, 839
Croom, S. M., Smith, R. J., Boyle, B. J., et al. 2004, MNRAS, 349, 1397
Croom, S. M., Richards, G. T., Shanks, T., et al. 2009,
MNRAS, 399, 1755
Croton, D. J. 2006, MNRAS, 369, 1808
Decarli, R., Falomo, R., Treves, A., et al. 2010, MNRAS,
402, 2453
Denney, K. D., Pogge, R. W., Assef, R. J., et al. 2013, ApJ,
775, 60
Degraf, C., Di Matteo, T., & Springel, V. 2010, MNRAS,
402, 1927
Di Matteo, T., Colberg, J., Springel, V., Hernquist, L., &Sijacki, D. 2008, ApJ, 676, 33
Draper, A. R., & Ballantyne, D. R. 2012, ApJ, 751, 72
Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012,
MNRAS, 420, 2662
Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468,
33
Fanidakis, N., Baugh, C. M., Benson, A. J., et al. 2012,

MNRAS, 419, 2797
Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9Fine, S., Croom, S. M., Hopkins, P. F., et al. 2008, MNRAS,
390, 1413
Fiore, F., Puccetti, S., Grazian, A., et al. 2012, A&A, 537,
A16
Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A,
459, 745
Gallo, E., Treu, T., Marshall, P. J., et al. 2010, ApJ, 714,
25
Garilli, B., Le Fèvre, O., Guzzo, L., et al. 2008, A&A, 486,
683
Gavignaud, I., Bongiorno, A., Paltani, S., et al. 2006, A&A,
457, 79
Gavignaud, I., Wisotzki, L., Bongiorno, A., et al. 2008,
A&A, 492, 637
Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539,
L13
Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79
Goulding, A. D., Alexander, D. M., Lehmer, B. D., & Mullaney, J. R. 2010, MNRAS, 406, 597
Greene, J. E. & Ho, L. C. 2005, ApJ, 630, 122
Greene, J. E. & Ho, L. C. 2007, ApJ, 667, 131
Hao, L., Strauss, M. A., Fan, X., et al. 2005, AJ, 129, 1795
Häring, N. & Rix, H.-W. 2004, ApJ, 604, L89
Hasinger, G., Cappelluti, N., Brunner, H., et al. 2007,
ApJS, 172, 29
Hasinger, G. 2008, A&A, 490, 905
Hasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441,
417
Heckman, T. M., Kauffmann, G., Brinchmann, J., et al.
2004, ApJ, 613, 109
Hickox, R. C., Jones, C., Forman, W. R., et al. 2009, ApJ,
696, 891
Hickox, R. C., Mullaney, J. R., Alexander, D. M., et al.
2014, ApJ, 782, 9
Hirschmann, M., Somerville, R. S., Naab, T., & Burkert,
A. 2012, MNRAS, 426, 237
Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS,
442, 2304
Ho, L. C. 2008, ARA&A, 46, 475
Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ,
654, 731
Hopkins, P. F. & Hernquist, L. 2009, ApJ, 698, 1550
Ilbert, O., Salvato, M., Le Floc’h, E., et al. 2010, ApJ, 709,
644
Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013,
A&A, 556, A55
Jahnke, K., Bongiorno, A., Brusa, M., et al. 2009a, ApJ,
706, L215
Jin, C., Ward, M., & Done, C. 2012, MNRAS, 425, 907
Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, ApJ, 533,
631
Kauffmann, G. & Heckman, T. M. 2009, MNRAS, 397, 135
Kelly, B. C., Vestergaard, M., & Fan, X. 2009, ApJ, 692,
1388
Kelly, B. C., Vestergaard, M., Fan, X., et al. 2010, ApJ,
719, 1315
Kelly, B. C., & Shen, Y. 2013, ApJ, 764, 45
Khandai, N., Di Matteo, T., Croft, R., et al. 2014,arXiv:1402.0888
Kollmeier, J. A., Onken, C. A., Kochanek, C. S., et al. 2006,
c 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2
ApJ, 648, 128
Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581
Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007,
ApJS, 172, 196
La Franca, F., Fiore, F., Comastri, A., et al. 2005, ApJ,
635, 864
Lanzuisi, G., Ranalli, P., Georgantopoulos, I., et al. 2014,arXiv:1409.1867
Läsker, R., Ferrarese, L., van de Ven, G., & Shankar, F.
2014, ApJ, 780, 70
Lauer, T. R., Tremaine, S., Richstone, D., & Faber, S. M.
2007, ApJ, 670, 249
Le Fèvre, O., Mellier, Y., McCracken, H. J., et al. 2004,
A&A, 417, 839
Le Fèvre, O., Vettolani, G., Garilli, B., et al. 2005, A&A,
439, 845
Le Fèvre, O., Cassata, P., Cucciati, O., et al. 2013, A&A,
559, A14
Li, Y.-R., Ho, L. C., & Wang, J.-M. 2011, ApJ, 742, 33
Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS,
172, 70
Lilly, S. J., Le Brun, V., Maier, C., et al. 2009, ApJS, 184,
218
Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS, 425, 623
Lutz, D., Sturm, E., Tacconi, L. J., et al. 2008, ApJ, 684,
853
Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ,
115, 2285
Maiolino, R., Shemmer, O., Imanishi, M., et al. 2007, A&A,
468, 979
Marconi, A. & Hunt, L. K. 2003, ApJ, 589, L21
Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS,
351, 169
Marleau, F. R., Clancy, D., & Bianconi, M. 2013, MNRAS,
435, 3085
Marshall, H. L., Tananbaum, H., Avni, Y., & Zamorani, G.
1983, ApJ, 269, 35
Marulli, F., Bonoli, S., Branchini, E., Moscardini, L., &Springel, V. 2008, MNRAS, 385, 1846
McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184
McCracken, H. J., Radovich, M., Bertin, E., et al. 2003,
A&A, 410, 17
McGill, K. L., Woo, J.-H., Treu, T., & Malkan, M. A. 2008,
ApJ, 673, 703
McLure, R. J. & Dunlop, J. S. 2004, MNRAS, 352, 1390
Menci, N., Fiore, F., Puccetti, S., & Cavaliere, A. 2008,
ApJ, 686, 219
Menci, N., Gatti, M., Fiore, F., & Lamastra, A. 2014, A&A,
569, A37
Merloni, A. 2004, MNRAS, 353, 1035
Merloni, A. & Heinz, S. 2008, MNRAS, 388, 1011
Merloni, A., Bongiorno, A., Bolzonella, M., et al. 2010,
ApJ, 708, 137
Merloni, A., Bongiorno, A., Brusa, M., et al. 2014, MNRAS, 437, 3550
Miller, B., Gallo, E., Treu, T., & Woo, J.-H. 2012, ApJ,
747, 57
Neistein, E., & Netzer, H. 2014, MNRAS, 437, 3373
Netzer, H., Lutz, D., Schweitzer, M., et al. 2007, ApJ, 666,
806
c 2014 RAS, MNRAS 000, 1–31

27

Netzer, H. & Trakhtenbrot, B. 2007, ApJ, 654, 754
Nicastro, F. 2000, ApJ, 530, L65
Nobuta, K., Akiyama, M., Ueda, Y., et al. 2012, ApJ, 761,
143
Novak, G. S., Ostriker, J. P., & Ciotti, L. 2011, ApJ, 737,26
Onken, C. A., Ferrarese, L., Merritt, D., et al. 2004, ApJ,
615, 645
Peng, C. Y., Impey, C. D., Rix, H., et al. 2006, ApJ, 649,
616
Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721,
193
Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004,
ApJ, 613, 682
Park, D., Woo, J.-H., Treu, T., et al. 2012, ApJ, 747, 30
Park, D., Woo, J.-H., Denney, K. D., & Shin, J. 2013, ApJ,
770, 87
Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008,
ApJ, 675, 234
Richards, G. T., Fan, X., Newberg, H. J., et al. 2002, AJ,
123, 2945
Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al.
2006a, ApJS, 166, 470
Richards, G. T., Strauss, M. A., Fan, X., et al. 2006b, AJ,
131, 2766
Rosario, D. J., Santini, P., Lutz, D., et al. 2012, A&A, 545,
A45
Rosario, D. J., Trakhtenbrot, B., Lutz, D., et al. 2013,
A&A, 560, A72
Ross, N. P., McGreer, I. D., White, M., et al. 2013, ApJ,
773, 14
Runnoe, J. C., Brotherton, M. S., & Shang, Z. 2012, MNRAS, 422, 478
Salucci, P., Szuszkiewicz, E., Monaco, P., & Danese, L.
1999, MNRAS, 307, 637
Schechter, P. 1976, ApJ, 203, 297
Schmidt, M. 1968, ApJ, 151, 393
Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010,
AJ, 139, 2360
Schramm, M., & Silverman, J. D. 2013, ApJ, 767, 13
Schulze, A., Wisotzki, L., & Husemann, B. 2009, A&A, 507,
781
Schulze, A. & Wisotzki, L. 2010, A&A, 516, A87+
Schulze, A. & Wisotzki, L. 2011, A&A, 535, A87
Schulze, A. 2011, PhD thesis, urn:nbn:de:kobv:517-opus54464
Schulze, A., & Wisotzki, L. 2014, MNRAS, 438, 3422
Sijacki, D., Vogelsberger, M., Genel, S., et al. 2014,arXiv:1408.6842
Silverman, J. D., Green, P. J., Barkhouse, W. A., et al.
2008, ApJ, 679, 118
Silverman, J. D., Lamareille, F., Maier, C., et al. 2009, ApJ,
696, 396
Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2009,
ApJ, 690, 20
Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2013,
MNRAS, 428, 421
Shen, Y., Greene, J. E., Strauss, M. A., Richards, G. T., &Schneider, D. P. 2008, ApJ, 680, 169
Shen, Y. 2009, ApJ, 704, 89
Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS,
194, 45

28

Schulze et al.

Shen, Y., & Kelly, B. C. 2012, ApJ, 746, 169
Shen, Y. 2013, Bulletin of the Astronomical Society of India, 41, 61
Small, T. A., & Blandford, R. D. 1992, MNRAS, 259, 725
Soltan, A. 1982, MNRAS, 200, 115
Soria, R., Fabbiano, G., Graham, A. W., et al. 2006, ApJ,
640, 126
Springel, V. 2005, MNRAS, 364, 1105
Steinhardt, C. L., & Elvis, M. 2010, MNRAS, 406, L1
Trakhtenbrot, B., & Netzer, H. 2012, MNRAS, 427, 3081
Treister, E., Urry, C. M., & Virani, S. 2009, ApJ, 696, 110
Tremaine, S., Gebhardt, K., Bender, R., et al. 2002, ApJ,
574, 740
Trump, J. R., Impey, C. D., Kelly, B. C., et al. 2009, ApJ,
700, 49
Trump, J. R., Impey, C. D., Kelly, B. C., et al. 2011, ApJ,
733, 60
Ueda, Y., Akiyama, M., Ohta, K., & Miyaji, T. 2003, ApJ,
598, 886
Ueda, Y., Watson, M. G., Stewart, I. M., et al. 2008, ApJS,
179, 124
Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104
Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001,
AJ, 122, 549
Vasudevan, R. V., & Fabian, A. C. 2007, MNRAS, 381,
1235
Veale, M., White, M., & Conroy, C. 2014, MNRAS, 445,
1144
Vestergaard, M., & Wilkes, B. J. 2001, ApJS, 134, 1
Vestergaard, M. 2002, ApJ, 571, 733
Vestergaard, M. & Peterson, B. M. 2006, ApJ, 641, 689
Vestergaard, M., Fan, X., Tremonti, C. A., Osmer, P. S.,& Richards, G. T. 2008, ApJ, 674, L1
Vestergaard, M. & Osmer, P. S. 2009, ApJ, 699, 800
Vignali, C., Mignoli, M., Gilli, R., et al. 2014, A&A, 571,
AA34
Wang, J.-M., Chen, Y.-M., & Zhang, F. 2006, ApJ, 647,
L17
Wang, J.-G., Dong, X.-B., Wang, T.-G., et al. 2009, ApJ,
707, 1334
Wisotzki, L., Christlieb, N., Bade, N., et al. 2000, A&A,
358, 77
Wolf, C., Wisotzki, L., Borch, A., et al. 2003, A&A, 408,
499
Woo, J., Treu, T., Barth, A. J., et al. 2010, ApJ, 716, 269
Woo, J.-H., Schulze, A., Park, D., et al. 2013, ApJ, 772, 49
Wyithe, J. S. B., & Loeb, A. 2003, ApJ, 595, 614
Xue, Y. Q., Brandt, W. N., Luo, B., et al. 2010, ApJ, 720,
368
Yu, Q. & Lu, Y. 2004, ApJ, 602, 603
Yu, Q. & Lu, Y. 2008, ApJ, 689, 732
Yu, Q., Lu, Y., & Kauffmann, G. 2005, ApJ, 634, 901
Yu, Q. & Tremaine, S. 2002, MNRAS, 335, 965
Yuan, F. 2007, The Central Engine of Active Galactic Nuclei, 373, 95

APPENDIX A: SELECTION FUNCTIONS
For the study of the statistical properties of the AGN samples and the determination of the distribution functions(AGN LF, BHMF, ERDF) we need to statistically accountfor objects that are not observed due to our survey layout.
This information is included in the selection function. Besides the flux limit of the survey, there are several levels
where we can loose objects. We here largely adopt the terminology used within the VVDS and zCOSMOS survey.
Target sampling rate (TSR): Based on the photometric input catalog only a fraction of the objects above
the survey limit have been observed spectroscopically. We
account for this fraction by correcting the effective areaof the respective survey j by this target sampling ratejfTSR
= Nspec /Nphot .
Spectroscopic success rate (SSR): For the objects
that have been observed spectroscopically the spectrum has
to allow the secure identification as a broad line AGN andthe measurement of a redshift. The SSR gives the probability of successful AGN identification which depends on theapparent magnitude, redshift and SED of the object. We include the flux limit of the respective survey into the SSR.
The SSR of the respective survey is a function of z andjmagnitude, fSSR
(m, z) = Nz /Nspec .
Mass measurement rate (MMR): For the determination of the BHMF we additionally require a reliable blackhole mass estimate for the AGN in the sample. The MMR
is the probability that for an object such a mass estimate ispossible. It can be defined by the fraction of objects in thesample that have a M• estimate. In general it will depend on
redshift and the S/N (or magnitude) of an object, i.e on thelocation of the broad emission line in the spectral coverage
and on the quality of the spectrum. As discussed below, we
here only account for a redshift dependence of the MMR,jfMMR(z) = NM• /Nz .
Redshift degeneracy weight (Wzd ): For the VVDS
sample some objects have a degenerate redshift and thus
only a certain probability fzd to be intrinsically within theredshift range we study. For all sources with a secure redshift
fzd = 1. We account for this redshift degeneracy weight via
an individual weight factor applied to each object in the
VVDS.
The combined selection function times area for each surjj
vey then follows as Ωje (m, z) = Ωj × fTSR
× fSSR(m, z) ×jfMMR (z).
A1VVDS

The target sampling rate (TSR) is given by the mean
sampling rate of a specific VVDS field, taken from
Garilli et al. (2008), namely 0.23, 0.24, 0.24, 0.22 and0.22 for CDFS, VVDS-0226-04, VVDS-1003+01, VVDS1400+05 and VVDS-2217+00, respectively. They also provide the effective area of the fields with 0.13, 0.48, 0.6,
0.9 and 3.0 deg2 , respectively. The corrected total area
is given
P by i the iproduct of area and TSR over all fieldsΩ= N
i=1 fTSR Ω .
The spectroscopic success rate (SSR) has been estimated in Gavignaud et al. (2006) by simulations of VIMOS
pointings and the identification of broad line AGN in thesimulated spectra.
For the mass measurement rate (MMR) we computed
the ratio of AGN with successful mass measurements overall detected AGN over several bins in redshift or magnitude.
We could not find any evidence in the sample for a z or magc 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2nitude dependence, thus we assume that any unaccounted
mild dependence is negligible for our study. We therefore use
a constant weight factor for the MMR, fMMR = 0.94.
To statistically account for the redshift degenerate objects we include them into the sample but weight them by
the probability that the broad line is in fact Mg ii. For someof the objects with an initially degenerate redshift a secure
redshift could be assigned from additional spectroscopy. The
broad line was resolved into Mg ii in 17 cases, C iii in 7 cases
and C iv for 4 objects. We use these empirical numbers tocompute the weight factorPfor the individual redshift degenerate objects via fzd =NMgII /Nline , where the sum isover all possible redshift solutions for the object, as given in
Gavignaud et al. (2008).

A2

zCOSMOSThe TSR is determined by dividing the number of objects
spectroscopically observed and the number of objects in therespective photometric input catalog. For the zCOSMOS
field the TSR is 0.966 for the compulsory targets and 0.550
for the random targets. The total TSR for our sample is thengiven by fTSR = Ntarget /Ntotal = (Nc + Nr )/(Nc /T SRc +Nr /T SRr ), where Nc and Nr are the number of compulsory
and random targets in our sample. This yields fTSR = 0.67
for our zCOSMOS sample. The effective area of the survey
is 1.648 deg2 .
zCOSMOS used the same instrument (VIMOS) and thesame integration time as the VVDS-wide survey. The maindifference is the different, higher resolution grism used (MRvs. LRRED), which will not significantly effect the detection probability of broad line AGN. Therefore we can adopt
the SSR for the VVDS-wide, as given in Gavignaud et al.
(2006), also for the zCOSMOS sample. We tested the consistency of the SSR, based on the X-ray detected AGN sample
from XMM-COSMOS. The location in the z-IAB diagramof X-ray detected type 1 AGN that have been targeted byzCOSMOS and are detected vs. not detected in the spectroscopy is fully consistent with our SSR, confirming ourused SSR, while the number statistics are too low to usethem for a direct empirical determination of the zCOSMOSSSR.
We estimate fMMR for our zCOSMOS sample as theratio of targets with mass measurements over all targets inseveral z bins. This mass success rate shows a redshift dependence in our sample, with a lower success rate at the lower
and upper redshift edges, while we achieve high completeness in between. We find a mean of 0.94, varying between
1.0 and 0.79 over our redshift range.

(Richards et al. 2006a). We compute the MMR over severalz bins. The mean is 0.97, with a constant value over most of
the redshift range and a slight decrease towards the upperedge.

APPENDIX B: RESULTS FROM THE
LUMINOSITY WEIGHTED 1/VMAX METHOD
We here discuss our results, adopting the directly the 1/Vmax
method. We re-emphasize that this approach introduces unaccounted sample incompleteness at the low mass end of the
BHMF and a potential overestimation of the space density
at the high mass end, due to not accounting for the effect ofvirial mass uncertainties.
In general the 1/Vmax method would be able to determine the intrinsic BHMF and ERDF when the proper
volume weights are applied (Schulze & Wisotzki 2010). The
accessible volume of an AGN in the sample is given by:
Z zmaxdVVmax =Ω(L, z)
dz ,
(B1)
dz
zmin
where Ω(L, z) is the luminosity selection function of the respective survey, defined in Appendix A. Using Equation B1
gives the proper volume weights for the determination of theluminosity function, but in general not for the determination of the BHMF or ERDF. For these, we would have to useΩ(M• , z), i.e. compute mass weighted volumes, or Ω(λ, z),compute λ weighted volumes. With Ω(L, z) = Ω(M• , λ, z),these volume weights are given by weighting Ω(L, z) with
the ERDF (for Ω(M• , z)) or the BHMF (for Ω(λ, z)). Because these distribution functions are a priori unknown beforehand, this approach is not practical, but at best serves
as a consistency check (see Schulze & Wisotzki 2010). Alternatively, given a sufficiently large sample size, the bivariate
distribution function can be determined from the 1/Vmax
method by binning in both, M• and λ (e.g. Yu et al. 2005;
Kauffmann & Heckman 2009).
On the other hand, the determination of the binned
BHMF and ERDF using the luminosity weighted 1/Vmax
method can serve as a useful tool. It provides a nonparametric and model independent estimate of the spacedensity, even if only as a lower limit, and therefore can guidethe refined but parametric analysis presented in Section 4.
The binned distribution function of the property X isgiven by (Schmidt 1968):
1 X 1
Φ(X) =.
(B2)∆X i Vmax,i
B1

A3

SDSS

For the SDSS the target sampling rate is already included in
the effective area, thus we can assume fTSR = 1. In contrast
to the VVDS and zCOSMOS samples the SDSS QSO sample
involves a color pre-selection. This pre-selection significantly
reduces the selection efficiency at redshifts 2 < z < 3.5
where QSO colors move into the stellar locus. For the redshift range studied in this paper the QSO selection efficiency
is not strongly affected by the color selection, i.e. the selection function is unity within the flux limit, 15.0 < i < 19.1
c 2014 RAS, MNRAS 000, 1–31

29

The AGN luminosity function

We first compute the AGN LF for the 3 samples, which
can be compared to previous results on the SDSS QSO
sample (Richards et al. 2006a; Shen & Kelly 2012) and the
VVDS AGN sample (Bongiorno et al. 2007). We derive the
AGN LF for absolute magnitudes in the SDSS i band defined at z = 2, i.e. X = Mi (z = 2) (Richards et al. 2006a;
Ross et al. 2013). As discussed in 3.2, we employ the SDSS
i-band K-correction from Richards et al. (2006a), where wetransform it to the CFHT/CFH12K I-band for VVDS andthe HST/ACS F814W filter for zCOSMOS.

30

Schulze et al.

1.1 <z <1.6

1.6 <z <2.1

log Φ

-6

-8

-10
-20

-22

-24

-26

-28

-20

Mi (z =2) [mag]

-22

-24

-26

-28

Mi (z =2) [mag]Figure B1. AGN LF in Mi (z = 2) from the SDSS sample (red triangles), the VVDS sample (orange circles) and the zCOSMOS sample(blue squares). The AGN LF is shown in two redshift bins, 1.1 < z < 1.6 (left panel) and 1.6 < z < 2.1 (right panel). The black solidline and black dashed line show our best fit PLE model and LDDE model to the combined sample. The open green squares show theAGN LF from Bongiorno et al. (2007), based on VVDS epoch 1.

Table B1. Best fit model parameters for the combined AGN LF from SDSS, VVDS and zCOSMOS over 1.1 < z < 2.1. We also compare
our maximum likelihood fit models to the binned AGN LFs of the three samples in three redshift bins via a χ2 test.
Evolution Modellog(Φ∗ )Mi∗ (z = 0)αβk1 / p 1

k2 / p 2

γ

zc,0

Mc,0

χ2 /ν

PLE

−6.14

−22.64

−1.63

−3.66

1.469

−0.317−−−

129/70

LDDE

−8.61

−26.57

−2.20

−3.64

5.77

−1.39

0.19

2.0

−27.19192/67The binned AGN LFs for the 3 samples are shown in
Figure B1 in two redshift bins. We find excellent agreement with the original SDSS DR7 AGN LF by Shen & Kelly(2012). We also find good agreement with the previous
VVDS epoch 1 AGN LF by Bongiorno et al. (2007), when
accounting for the difference in the applied K-correction.
The zCOSMOS AGN LF is in good agreement with the
VVDS result, spanning a similar luminosity range.
Our data improves the statistical accuracy both on the
faint end and on the bright end of the AGN LF, comparedto the work of Bongiorno et al. (2007). Thus, we here alsopresent a determination of the combined AGN LF for thethree samples. We fit a parametric AGN LF model to thefull unbinned data over 1.1 < z < 2.1 by performing amaximum likelihood fit (Marshall et al. 1983). As discussed
in section 4.1, we combine the three surveys, following thestrategy outlined by Avni & Bahcall (1980).
Following common practice we fit the AGN LF as a
double power law:Φ∗
ΦL (M, z) = 0.4(α+1)(M −M ∗ (z)) M 0.4(β+1)(M −M ∗ (z)) .
10
+ 10
(B3)
The normalization of the AGN LF is determined by the ratio of the number of observed AGN to the number of AGN
expected from the best fit model, Φ∗M = Nobs /Nmod . Overthe restricted redshift range we are probing (1.1 < z < 2.1)
the redshift evolution in the AGN LF can be well approxi-

mated with a simple pure luminosity evolution (PLE) model(Croom et al. 2004; Ross et al. 2013). We follow Boyle et al.
(2000) by modeling the redshift evolution in the break of the
LF with a second order polynomial:Mi∗ (z) = Mi∗ (z = 0) − 2.5(k1 z + k2 z 2 ) .

(B4)The best fit PLE model is shown as solid black line in Figure B1 and the best fit parameters are given in Table B1.
Additionally, we also adopt a luminosity dependent density
evolution model (LDDE), following the parameterization in
Bongiorno et al. (2007) (see also Hasinger et al. 2005). The
best fit LDDE model is shown as dashed black line in Figure B1 and we give the best fit parameters in Table B1. Bothmodels provide a reasonably good fit to our data. Since wehere only probe a limited redshift range, further discussions
on the redshift evolution of the AGN LF are beyond thescope of this work.

B2

The binned BHMF and ERDF

Next, we determine the binned, luminosity weighted BHMF
(X = log M• ) and ERDF (X = log λ) for the 3 samples. The results are shown in Figure B2 and Figure B3.
The binned BHMF and ERDF for the SDSS sample shows
an apparent turnover towards low masses and low Eddington ratio. This is an artifact of the incompleteness introduced to the binned estimate due to the use of luminosityc 2014 RAS, MNRAS 000, 1–31

The cosmic growth of the active black hole population at 1 < z < 2
1.1 <z <1.6

1.6 <z <2.1

log Φ •

-4

-6

-8

7

8

9

10

7

log M • [M ⊙ ]

8

9

10

log M • [M ⊙ ]Figure B2. Binned, luminosity weighted active black hole massfunction (BHMF) in two redshift bins, 1.1 < z < 1.6 and1.6 < z < 2.1, for the SDSS sample (red triangles), VVDS sample (orange circles) and zCOSMOS sample (blue squares). While
all three distribution function estimates are affected by inherent
incompleteness, this is most severe for the SDSS sample, due tothe brighter flux limit.

1.1 <z <1.6

1.6 <z <2.1

log Φλ

-4

-6

-8

-2

-1

0

log λ

-2

-1

0

log λ

Figure B3. Binned, luminosity weighted Eddington ratio distribution function (ERDF) in two redshift bins. Symbols are thesame as for the BHMF in Figure B2.

weighted volumes instead of mass weighted or Eddington
ratio weighted volumes. This becomes evident when adding
the VVDS and zCOSMOS samples. Due to the deeper flux
limit they are less affected by this kind of incompletenessand thus have a higher space density at low masses and Eddington ratios. Thus, incompleteness corrections are smaller
than for the SDSS, while still being required. In particular
the BHMFs from the VVDS and zCOSMOS also show anapparent turnover, which occurs at lower masses than for
the SDSS sample. Based on the luminosity weighted binned
estimates it is less clear if this turnover is intrinsic or due
to incompleteness in the luminosity weighted BHMF. The
binned zCOSMOS BHMF has a more prominent turnoverat log M• < 8 than the VVDS BHMF, caused by the deep
fields included in the VVDS with a flux limit of IAB = 24.
Restricting the VVDS sample to the wide fields only gives aluminosity weighted BHMF fully consistent with the zCOSMOS result.
A further limitation of the binned approach is the factthat it does not correct for the uncertainty in the virial blackhole mass estimates (Kelly et al. 2009; Shen & Kelly 2012).
This will generally broaden the observed M• and λ distributions and in particular lead to an overestimate of the spacedensity at the steeply decreasing parts of the distributionfunctions, i.e. at the high mass end in the BHMF and athigh λ in the ERDF. In our maximum likelihood approach
this effect has been taken into account.
c 2014 RAS, MNRAS 000, 1–31

31