{'fraction_non_alphanumeric': 0.03628623954226841, 'fraction_numerical': 0.01786293612796036, 'mean_word_length': 4.645152339197514, 'pattern_counts': {'":': 1, '<': 21, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A synthesis model for AGN evolution: supermassive black holes growth and feedback modes

We present a comprehensive synthesis model for the AGN evolution and the growth of supermassive black holes in the Universe. We solve the continuity equation for SMBH mass function using the locally determined one as a boundary condition, and the hard X-ray luminosity function as tracer of the AGN growth rate distribution, supplemented with a luminosity-dependent bolometric correction and an absorbing column distribution. Differently from most previous semi-analytic and numerical models, we do not assume any specific distribution of Eddington ratios, rather we determine it empirically by coupling the mass and luminosity functions. SMBH show a very broad accretion rate distribution, and we discuss the consequences of this fact for our understanding of observed AGN fractions in galaxies. We confirm previous results and demonstrate that, at least for z<1.5, SMBH mass function evolves anti-hierarchically, i.e. the most massive holes grew earlier and faster than less massive ones. For the first time, we find hints of a reversal of such a downsizing behaviour at redshifts above the peak of the black hole accretion rate density (z~2). We also derive tight constraints on the (mass weighted) average radiative efficiency of AGN: we find that 0.065<\xi_0 \epsilon_{rad}<0.07$, where \xi_0 is the local SMBH mass density in units of 4.3x10^5 M_sun Mpc^{-3}. We trace the cosmological evolution of the kinetic luminosity function of AGN, and find that the overall efficiency of SMBH in converting accreted rest mass energy into kinetic power, ranges between 3 and 5 \times 10^{-3}. Such a ``kinetic efficiency'' varies however strongly with SMBH mass and redshift, being maximal for very massive holes at late times, as required for the AGN feedback by many galaxy formation models in Cosmological contexts. (Abriged)


INTRODUCTION
Soon after their discovery, it was realized that Quasars (QSOs) were a strongly evolving Cosmological population, that could effectively be used as a tracer to study the structural properties of the evolving Universe (Longair 1966). The physical understanding of QSOs and of their lower-luminosity counterparts (both generally called Active Galactic Nuclei, AGN) as accreting supermassive black holes (SMBH) immediately led to speculations about the presence of their dormant relics in the nuclei of nearby galaxies, with Soltan (1982) first proposing a method to estimate the mass budget of SMBH based on the demographics and evolutionary paths of observed QSO/AGN.
The search for the local QSO relics via the study of their dynamical influence on the surrounding stars and gas carried out since the launch of the Hubble Space Telescope (see Richstone et al. 1998;Ferrarese et al. 2008, and reference therein) led ultimately to a shift of paradigm in our understanding of the role of black holes in the formation of structures in the Universe. The discovery of tight scaling relations between SMBH masses and properties of the host galaxies' bulges (Magorrian et al. 1998 Hopkins et al. 2007a) clearly pointed to an early co-eval stage of SMBH and galaxy growth. However, local scaling relations have proved themselves unable to uniquely determine the physical nature of the SMBHgalaxy coupling, and a large number of feedback models have been proposed which can reasonably well reproduce these relations (Silk & Rees 1998;Fabian 1999;Kauffmann & Haehnelt 2000;Umemura 2001;Cavaliere & Vittorini 2002;Wyithe & Loeb 2003;Granato et al. 2004;King 2005;Sazonov et al. 2005;Escala 2006). Most likely, it will only be by studying in detail the redshift evolution of these scaling relations that we will gain critical insight on the nature of the interaction between growing black holes and their host galaxies and on the nature of the accretion process and of jet formation in various regimes. In fact, many have recently attempted to push the study of scaling relations to z > 0, despite the large observational biases (Merloni, Rudnick & Di Matteo 2004;McLure et al. 2006;Shields et al. 2006;Peng et al. 2006;Wu 2007;Treu et al. 2007;Salviander et al. 2007;Shen et al. 2008;Schramm, Wisotzki & Jahnke 2008), but the question of if, and, in case, by how much did these relations evolve with redshift has so far eluded any unambiguous answer.
For all these reasons, it seems to us necessary to make the best effort possible to reconstruct the history of SMBH accretion in order to follow closely the evolution of the black hole mass function. That will serve as a benchmark to test various models for SMBH cosmological growth as well as those for the black hole-galaxy coevolution.
As opposed to the case of galaxies, where the direct relationship between the evolving mass functions of the various morphological types and the distribution of star forming galaxies is not straightforward due to the never-ending morphological and photometric transformation of the different populations, the situation is much simpler in the case of SMBH. For the latter we can assume their evolution is governed by a simple continuity equation (Cavaliere, Morrison & Wood 1971;Small & Blandford 1992;Yu & Tremaine 2002), where the mass function of SMBH at any given time can be used to predict that at any other time, provided the distribution of accretion rates as a function of black hole mass is known (see section 3 for more details). An analogous continuity equation could in principle be solved for the evolution of black hole spin (and thus for the overall accretion efficiency), whose time derivative, however, depends both on the accretion rate and on the merger rate evolution (Berti and Volonteri 2008;King, Pringle and Hofmann 2008).
Consequently, the simple physical and mathematical structure of a black hole implies that their intrinsic properties are essentially determined by three physical quantities only: mass, spin and accretion rate. Leaving aside a few open issues related to the physics of strongly super-Eddington flows, accretion theory is a mature enough field to allow a relatively robust mapping of observables (spectral energy distributions, variability patterns, etc.) into these physical quantities associated to SMBH growth. This reduces the problem of unveiling the census of growing black holes to that of a wavelength-dependent assessment of selection biases, bolometric corrections, and accurate determinations of the distribution of extrinsic observational properties, the level of intervening obscuration being the main one.
The fact that obscuration is a crucial factor in determining the observed properties of AGN was recognized early after the discovery of the X-ray back-ground (XRB) radiation (Setti & Woltjer 1989). Subsequent generation of synthesis models of the XRB (Madau, Ghisellini & Fabian 1993;Comastri et al. 1995;Wilman & Fabian 1999;Gilli, Salvati & Hasinger 2001;Gilli, Comastri & Hasinger 2007), following closely on the footsteps of increasingly larger and deeper surveys, have progressively reduced the uncertainties on the absorbing column density distribution, which, coupled to the observed X-ray luminosity functions, provide now an almost complete census of the Compton thin AGN (i.e. those obscured by columns NH < σ −1 T ≃ 1.5 × 10 24 cm −2 , where σT is the Thomson cross section). Indeed, this class of objects dominates the counts in the X-ray energy band where almost the entire XRB radiation has been resolved into individual sources (Hasinger et al. 1998;Rosati et al. 2002;Alexander et al. 2003;Moretti et al. 2003;Worsley et al. 2005). Accurate determinations of the XRB intensity and spectral shape (for the most recent results see Churazov et al. 2007;Ajello et al. 2008), coupled with the resolution of this radiation into individual sources , allow very sensitive tests of how the AGN luminosity and obscuration evolve with redshift (the so-called 'X-ray background synthesis models' see e.g. Gilli et al. 2007, and references therein). In much the same way, accurate determinations of the local SMBH mass density and of the AGN (bolometric) luminosity functions, coupled with accretion models that specify how the observed AGN radiation translates into a black hole growth rate, allow sensitive tests of how the SMBH population (its mass function) evolves with redshift. This is the aim of this work, which, by analogy, we name 'AGN synthesis modelling'.
Reconstructing in detail the history of accretion onto supermassive black holes is also of paramount importance for understanding the effects that the energy released by growing black holes has on its environment. As we mentioned before, the various relations between black hole mass and the physical properties of host galaxy's bulge suggest that the growth of black holes and large scale structure is intimately linked. Moreover, X-ray observations of galaxy clusters directly show that black holes deposit large amounts of energy into their environment in response to radiative losses of the cluster gas (Bîrzan et al. 2004;Fabian et al. 2006;Allen et al. 2006). From these (and other) studies of the cavities, bubbles and weak shocks generated by the radio emitting jets in the intra-cluster medium (ICM) it appears that AGN are energetically able to balance radiative losses from the ICM in the majority of cases (Rafferty et al. 2006). Also, numerical simulations of AGNinduced feedback have recently shown that mechanical feedback from black holes may be responsible for halting star formation in massive ellipticals, perhaps explaining the bimodality in the color distribution of local galaxies , as well as the size of the most massive ellipticals (De Lucia & Blaizot 2007).
All these arguments, however, hinge on the unknown efficiency with which growing black holes convert accreted rest mass into kinetic (ǫ kin ) and/or radiative (ǫ rad ) power. Constraints on these efficiency factors are therefore vital for all models of black hole feedback, and we aim here at pinning them down as accurately as currently possible.
The structure of this paper is the following: we will first start with a description of our current knowledge about the various modes in which a black hole may grow ( § 2). This will be used to introduce the function that relates the observed (bolometric) luminosity of an AGN to the true accretion rate. In section 3 we will outline the methods used to calculate the evolution of the SMBH mass (and accretion rate) function, while in section 4 we will spell out and discuss the fundamental assumptions made in this work, in particular those about the local mass function of SMBH, the adopted bolometric correction and the incidence of obscured AGN. Our results on the growth of the SMBH population will be presented in § 5. There, we will first discuss "integral" constraints ( § 5.1) and then "differential" ones ( § 5.2-5.4), which will reveal interesting features of the so-called anti-hierarchical (or 'downsizing') evolution of SMBH. Section 6 will instead be devoted to the study of the kinetic luminosity function of AGN and of its evolution. At the end of it ( § 6.2) we will be able to assess the relative importance of kinetic and radiative feedback from growing black holes as a function of both mass and redshift. Finally, we will discuss our results and draw our conclusions in section 7.

THE MODES OF BLACK HOLE GROWTH
Galactic (stellar mass) black holes in X-ray binaries, either of transient or persistent nature, are commonly observed undergoing so-called transitions, i.e. dramatic changes in their spectral and variability properties associated to changes in the way matter accretes onto them. The study of these various spectral states, either in individual objects or in ensemble of sources, has solidified in the last few years into a general framework of the various physical modes in which accretion onto a black hole can take place (Remillard & McClintock 2006;Malzac 2007;Done, Gierliński & Kubota 2007). We discuss these modes, and their AGN analogues, in this section, paying particular attention to their partition in terms of the value of the Eddington ratio λ ≡ L bol /L Edd , where L bol is the bolometric luminosity and L Edd = 4πGMBHmpc/σT ≃ 1.3 × 10 38 (MBH/M⊙) erg s −1 is the Eddington one.
There are at least three well defined spectral states. In the low/hard state the emission is dominated by a hard X-ray power-law with an exponential cutoff at about few 100 keV. Radio emission is always detected in this state (Fender 2006) usually with flat or inverted spectrum, which is commonly interpreted as due to a compact, persistent, self-absorbed jet of low intrinsic power. Interestingly, the radio luminosity of these jet cores seems to correlate quite tightly with the X-ray luminosity (Gallo, Fender, & Pooley 2003). The spectrum of the high/soft state, instead, is dominated by a thermal component likely originated in a standard Shakura & Sunyaev (1973) accretion disc, while in the intermediate (or very high) state, associated with transitions between hard and soft states, and often occurring at a source's highest flux level, both a thermal and a steep power-law component substantially contribute to the spectrum 1 . Maccarone (2003) has shown that the state change in these systems generally occurs at (bolometric) luminosities of about λ ≃ a few ×10 −2 . The hard-tosoft transition is also accompanied by a "quenching" of the steady radio emission observed in the low/hard state, following rapid flaring events in the radio band, often associated with fast ejection of matter moving at relativistic speeds and emitting optically thin synchrotron radiation (Fender, Belloni & Gallo 2004).
Given the many analogies between the high-energy spectra of galactic black holes and AGN, a similar phenomenology has long been searched for in active supermassive black holes. Merloni, Heinz and Di Matteo (2003;MHD03), following earlier suggestions by Falcke and Biermann (1996), presented a thorough discussion on the scale-invariant properties of black hole coupled accretion/jet systems, based essentially on a generalization of the radio-X-ray correlation of low/hard state X-ray binaries to SMBH of varying masses. This led to the discovery of a so-called "fundamental plane" of active black holes (see also Falcke, Körding and Markoff 2004;Körding, Falcke and Corbel 2006;Merloni et al. 2006;Fender et al. 2007). This correspondence between low/hard state X-ray binaries and low-luminosity AGN is confirmed by a number of recent studies. Detailed multi-wavelength observations of nearby galaxies have revealed a clear tendency for lower luminosity AGN to be more radio loud as the Eddington-scaled accretion rate decreases (see Ho 2002, and references therein;Nagar, Falcke and Wilson 2005). The above mentioned scaling relations with black hole X-ray binaries have thus helped to identify AGN analogues of low/hard state sources, in which the radiative efficiency of the accretion flow must be low (Churazov et al. 2005;Pellegrini 2005a;Chiaberge, Capetti & Macchetto 2005;Hardcastle, Evans & Croston 2007).
In this work, we will use the kinetic jet power as a key physical variable that identifies and characterizes the nature of various accretion modes. In fact, Merloni and Heinz (2007), using a sample of nearby AGN for which BH masses, nuclear radiative and kinetic luminosities and outer boundary condition (i.e. accretion rate at the Bondi radius) were simultaneously known, have found a clear relationship between Eddington-scaled kinetic power and bolometric luminosity, given by: log(W kin /L Edd ) = (0.49 ± 0.07) log λ − (0.78 ± 0.36). The measured slope suggests that these objects are in a radiatively inefficient accretion mode. Moreover, the observed correlations are in very good agreement with theoretical predictions of adiabatic accretion models with strong outflows, the relative strength of which increases with decreasing accretion rate (Blandford & Begelman 1999;Blandford & Begelman 2004). This is a very important finding, as the existence of jet dominated accreting SMBH bears important consequences for our understanding of the AGN feedback in a cosmological context. The specific importance of mechanical (as opposed to radiative) feedback for the heating of baryons within the deepest dark matter potential wells has recently been acknowledged by semi-analytic modelers of cosmological structure formation (Croton et al. 2006;Bower et al. 2006). Within these schemes, because the bright quasar population peaks at too early times, a socalled "radio mode" of SMBH growth is invoked in order to regulate both cooling flows in galaxy clusters and the observed sizes and colors of the most massive galaxies in the local Universe Croton et al. 2006).
Here, we also consider such a mode of accretion, where most of the released energy is in kinetic form, but we uniquely associate it to black holes accreting below a certain critical rate in Eddington units, rather than to the state of the fuelling gas (Hardcastle, Evans & Croston 2007), or to the properties of the very large scale environment (as, for example, the host dark matter halo size). In order to avoid confusion, we call this "low kinetic (LK)" mode, to distinguish it not only from the bright "high radiative (HR)" mode (the 'quasar' mode of the recent cosmology jargon, associated to radio quiet QSOs, analogous to the high/soft state of X-ray binaries), but also from the bright radio loud quasars (of which 3C273 is the prototype), characterized by both radiatively efficient accretion flows and powerful jets, that we term "high kinetic (HK)" mode, as we discuss below.
The most luminous sources, as those falling into the standard definition of Quasars and broad lined AGN, are most likely accreting at a high rate, close to the Eddington limit (McLure & Dunlop 2004;Kollmeier et al. 2006). Also their spectral energy distribution, usually dominated by the so-called Big Blue Bump (BBB: quasi thermal UV emission from an optically thick standard accretion disc, see e.g. Malkan 1983;Laor 1990), indicates that above a certain critical value of the ratio λcr, an accretion mode transition should take place to what is usually described as a standard, geometrically thin and optically thick accretion disc (Shakura & Sunyaev 1973). Evidence for such a transition taking place in AGN has been presented by Maccarone, Gallo & Fender (2003), Jester (2005) and Greene, Ho and Ulvestad (2006).
Thus, according to our scheme, black holes accreting above the critical rate come into two different physical states, one with (HK), and one without powerful jets (HR). The reason for this dichotomy (and even its reality, see e.g. Cirasuolo et al. 2003) has been subject to countless arguments since the very early years of the first radio loud QSOs discoveries, but it is still matter of debate (for example, see Sikora, Stawarz & Lasota 2007;Kaiser & Best 2008;Blundell 2008 for some recent ideas), with the two most popular options for the observed dichotomy being the influence of the black hole spin and a change in accretion mode.
As a working hypothesis, in line with what originally proposed in MHD03 (see also Nipoti et al. 2005;Körding, Jester and Fender 2006;Blundell 2008), we consider here that powerful radio jets are episodic and transients events during the time a black hole spends accreting above the critical rate (i.e. close to the Eddington limit). This idea, that is easy to incorporate into any population synthesis scheme, bears resemblance with what actually observed for transients black hole X-ray binaries, which often displays very powerful and rapid radio flares at the peak of their outbursts (see Fender et al. 2004, and references therein).

A schematic view of radiative, kinetic and accretion efficiency
We adopt a simple, physically motivated functional form for the λ vs.ṁ function (whereṁ is the accretion rate onto the black hole in Eddington unitsṁ = ηṀ c 2 /L Edd , with η accretion efficiency) of a broken power-law, bridging the low accretion rate (radiatively inefficient) regime (λ ∝ṁ 2 ) and the high accretion rate one. For radiatively efficient sources, we should assume, by definition, that the bolometric luminosity is simply proportional to the accretion rate, λ ∝ṁ. The overall normalization is found by imposing continuity at λcr and depends on the adopted value of this critical accretion rate. Based on the results of Maccarone et al. (2003), Greene et al. (2006), Merloni and Heinz (2007), we adopt here a value of λcr = 3 × 10 −2 (see also the discussion in Hopkins, Narayan & Hernquist 2006). A schematic view of these various modes in the accretion rate vs. released power plane is shown in Fig. 1. The Figure represents qualitatively the main features of the three-modes scheme we have outlined here: below the critical rate, only one mode (LK) is allowed, where most of the power emerges in kinetic form and there is a quadratic relationship between accretion rate and (bolometric) radiative luminosity; above it, two different modes are possible: one where kinetic and radiative power are comparable and high (HK), and one where kinetic power is quenched (HR). The exact slopes of the kinetic power vs. accretion rate relation are set, implicitly, by the relationship between radio core luminosity, kinetic power and X-ray luminosity (the "fundamental plane" of active black holes), as well as by the hard X-ray bolometric correction function, that we discuss in § 3 and 4. This simple relation, together with a bolometric correction, allows the determination of the dimensionless accretion rate for any object for which both the nuclear luminosity, Li, (in any band) and the black hole mass are measured. Conversely, the accretion luminosity of a black hole of mass M , growing at a rateṀ can be obtained from the above relationship. In § 3 we will describe in detail the method we employ to unveil the evolution of the black hole mass function based on the evolution of the intrinsic hard X-ray (2-10 keV) luminosity function, making use of the relation Li = Li(Ṁ , M ) discussed above.
As a final remark in this section, we come back here to the well known distinction between accretion efficiency, and radiative efficiency, that will be important in the following. The former represents the maximal amount of potential energy that can be extracted, per unit rest mass energy, from matter accreting onto the black hole. This quantity, η(a), depends on the inner boundary condition of the accretion flow only, and, in the classical no-torque case (Novikov & Thorne 1973) is a function of BH spin, a, only, ranging from η(a = 0) ≃ 0.057 for Schwarzschild (non-spinning) black holes to η(a = 1) ≃ 0.42 for maximally rotating Kerr black holes.
On the other hand, the radiative efficiency, ǫ rad ≡ L bol /Ṁ c 2 depends both on the accretion efficiency (i.e. on the inner boundary condition of the accretion flow) and on the nature of the accretion flow itself. In the case of a simple broken power-law shape that we have assumed for the λ −ṁ relation, we have that: whereṁcr = λcr is the critical Eddington-scaled accretion rate above which the disc becomes radiatively efficient. As we will show in § 5, the classical argument due to Soltan (1982) can be used to put constraints on the (mass weighted) average value of ǫ rad , by comparing the evolution of AGN/QSO luminosity function with estimates of the local relic black hole mass functions. Constraints on η(a) , and thus on the average spin of the SMBH population, can then be derived from the use of eq. (1) (Hopkins, Narayan & Hernquist 2006). The relationship between these average values η and ǫ rad depends both on the value of λcr (higher values of the critical accretion rate increase the number of objects which are in a radiatively inefficient regime, thus reducing the ratio between radiative and accretion efficiencies, and viceversa) and on the shape and evolution of the luminosity functions themselves. However, in the remain of the paper we will assume that eq. (1) holds, and will not discuss explicitly the dependence of the η / ǫ rad ratio on λcr. Moreover, for the sake of clarity, we will refer in the following to the mass weighted averages of the efficiencies simply with the symbols η and ǫ rad .

METHODS
In the simplest case of purely accretion driven evolution (i.e. assuming mergers do not play an important role in shaping the black holes mass function, but see § 7 for more on this point), the simultaneous knowledge of the black hole mass function and of the Eddington ratio distribution function, Φ M,λ (z) (i.e. the number of objects with given mass M and Eddington ratio λ per comoving volume and logarithmic interval), would completely determine the evolutionary solution for the mass function, as the two distributions must be coupled via a continuity equation (Small & Blandford 1992;Yu & Tremaine 2002;Marconi et al. 2004;Merloni 2004, hereafter M04;Shankar, Weinberg & Miralda-Escudé 2008) ∂ΦM(M, t) ∂t to be solved given an appropriate boundary condition. Here ΦM is the black hole mass function, and Ṁ (M, t) represents the average accretion rate (in physical units) for a SMBH of mass M at the time t, and can be uniquely determined once the black hole Eddington ratio distribution function, Φ M,λ (z), is known. The method used to solve eq. (2) follows closely that applied in M04, with some modifications and improvements, as we describe it in detail here. As in M04, we argue that the most robust constraint on the SMBH mass function at any redshift can be obtained by integrating the continuity equation backwards in time, using the local mass function as a boundary condition. Moreover, as opposed to all other previous calculation of the SMBH mass function evolution, we do not make any simplifying assumption about the distribution of the accretion rates for different black hole masses and redshifts, but we do calculate these distributions from the multi-wavelength luminosity functions of the AGN population. The guiding principle is that, within each of the three physically different modes of accretion described in the previous section (LK, HR, HK), there exist well defined relationships between an object luminosity in various bands and the physical parameters needed to solve eq. (2), i.e. the black hole mass and its accretion rate. In section 2 we have already presented the relationship between accretion rate, mass and bolometric luminosity in these various modes. The second relationship, needed in order to determine the evolution of the AGN jet (and thus of their kinetic energy feedback) is that between radio luminosity, bolometric luminosity and BH mass, that we discuss here.
In the 'low kinetic' (LK) mode, i.e. for all sources with Eddington ratios λ less than the critical value λcr, we adopt the now well established result, first described in MHD03, that a 'fundamental plane' relation between (core radio) jet luminosity, (accretion-powered) X-ray luminosity and black hole mass can indeed be defined for those objects (but see, for example, Panessa et al. 2007 for a contrasting view). We write it here as where M is the black hole mass in solar units, LR is the intrinsic (i.e. un-beamed) radio luminosity of the jet core at 5 GHz and LX stands for the intrinsic (i.e. unabsorbed) X-ray luminosity in the 2-10 keV band (both in units of erg s −1 ), which we uniquely relate to the bolometric luminosity by the relation of eq. (21) of Marconi et al. (2004). We adopt for the correlation coefficients the values recently determined in Merloni and Heinz (2007) based on a more accurate study of the effects of relativistic beaming, as well as a more accurate selection of the parent sample, as first discussed in Körding, Falcke and Corbel (2006). Specifically, we assume ξRX = 0.62, ξRM = 0.55 and βLK = 8.6. We also assume that this relationship has an intrinsic scatter of 0.6 dex, a value smaller than that found in MHD03, but in line with that discussed in  for samples of jet dominated objects only.
Let us now discuss the relationship between radio, X-ray luminosity and black hole mass for objects accreting at high rates. Irrespective of the physical explanation preferred for the existence of a radio loud/quiet (HK/HR) dichotomy, the problem we face for our method is that there are no well defined relationships to be found in the literature between radio luminosity, mass and bolometric luminosity which are specific for black holes accreting above the critical rate, λcr. One reason for this, as we have argued before in Merloni et al. (2005), is that, because such modes can exist only in very limited range of the parameter space (i.e. for λcr < ∼ λ < ∼ 1, where λ = 1 corresponds to the Eddington limit), it is almost impossible to unravel the mass from the accretion rate dependence, given the usual uncertainties on SMBH mass estimates and the likely presence of intrinsic scatter in any such relation. Thus, our treatment of the HK sources (classical radio loud QSOs) should necessarily be considered at best as an educated guess. For the sake of simplicity, we assume that their intrinsic radio core luminosity scales with the bolometric luminosity, so we fix log LR = log L bol + ξHM log M + βHK,HR. (4) We work out the mass scaling ξHM (the same for HK and HR modes) and the normalization of the HK (radio loud) mode assuming continuity with the fundamental plane relation at the critical accretion rate, while for the HR (radio quiet) mode, we simply take a luminosity 10 3 times smaller (βHR = βHK − 3), so that the radio luminosity in the HR (radio quiet) mode is effectively quenched. We obtain ξHM = 0.02 (mass dependency is negligible), and βHK = −4.9, in relatively good agreement with the standard radio loud SED of Elvis et al. (1994). Also for this relationship we assume an intrinsic scatter of 0.6 dex. Once fixed the relationship between LR, LX and M for all values of the accretion rate (i.e. for the three main modes LK, HR and HK), we can proceed in extracting the distribution of accretion rate as a function of SMBH mass and redshift. In order to do so, we take the observed luminosity functions for LR (intrinsic radio core emis-sion) and LX (intrinsic 2-10 keV luminosity) and the SMBH mass function. We then look for the joint distribution ΨC(LR, LX, M ), i.e. the function that quantifies the differential number of objects (per unit comoving volume) with given mass M , radio luminosity LR and X-ray luminosity LX. We do so by imposing the relations (3) and (4), with their own intrinsic scatter, and minimizing the differences between the projections of ΨC(LR, LX, M ) onto the three axes identified by LR, LX and M and the observed luminosity and mass functions. For high accretion rate objects (i.e. those with λ > λcr), we assume that the HK (radio loud) objects represent only 10% of the total, so that the number ratio of HR/HK AGN is about 9/1. We also impose the Eddington limit by assuming that the number of objects with λ > 1 are exponentially suppressed. By construction, then, our joint distribution function ΨC (LR, LX, M ) reproduces the observed luminosity (in the radio and X-ray bands) and mass functions, and is consistent with the scaling relations (3) and (4).
Once the joint distribution ΨC is found, it is straightforward to infer the accretion rate distribution function: each combination (LR, LX and M ) correspond to a unique value of the accretion rateṀ , and the number of objects per unit comoving volume (per logarithmic interval) with a given value of M andṀ is given by: (5) where LX = LX(M,Ṁ ) is the relation between intrinsic X-ray luminosity and mass and accretion rate discussed in section 2. Thus, the knowledge of the joint distribution function ΨC , complemented with the theoretical relationship between luminosity and accretion rate for the various modes (X(M,Ṁ ), see § 2) is sufficient to derive the distribution of accretion rates for black holes of different masses. The average Ṁ (M ) is what is then needed to solve eq. (2). In this way, the combined knowledge of the SMBH mass function and of the AGN luminosity function in at least one band at a given moment in time completely determine the boundary conditions needed to integrate the continuity equation. In practice, starting from the mass function and the accretion rate distribution at a given redshift zi, it is possible to derive the new black hole mass function at redshift zi + dz, ΦM(M, zi + dz), by just subtracting the mass accreted in the time interval dt = dz(dt/dz) calculated according to the accretion rate function of redshift zi. This new mass function can then be used together with the radio and X-ray luminosity functions at the same redshift, φR(LR, z + dz) and φX(LX, z + dz), to obtain the new joint distribution function ΨC(LR, LX, M, z+dz), and therefore the new accretion rate function, and so on, until a solution of eq. (2) is found in the interval zi + ∆z, where the value of ∆z depends solely on the availability of good enough data (luminosity functions) to meaningfully determine the joint distribution ΨC.

ASSUMPTIONS AND OBSERVABLES: MASS AND LUMINOSITY FUNCTIONS
According to the procedure we have outlined in § 3 above, in order to obtain meaningful constraints on the accretion rate distribution, we have to make the following assumptions: • The local BH mass function is known accurately down to a certain SMBH mass; • The bolometric corrections are known at all redshifts for AGN of all luminosities, so that meaningful bolometric luminosity functions can be derived; • The selection function of hard X-ray surveys provides us with a complete census of the AGN population, i.e. is not strongly biased against obscured sources. Although this is obviously not the case for the so-called Compton Thick AGN, i.e. those obscured by such a large column (NH > σ −1 T ) that also hard X-rays cannot be transmitted to the observer, we will proceed under the assumption that the Compton Thick AGN fraction can be meaningfully constrained on the basis of the X-ray background models, and verify, a posteriori that their contribution to the overall growth is small; • The bolometric LF of AGN, once extrapolated to low enough luminosities, recovers the entirety of the local SMBH population, as described by the adopted mass function. This is equivalent to the statement that the AGN fraction approaches unity towards very low nuclear luminosities, once incompleteness of any specific LF is accounted for.
We will discuss these issues in turn below.

Local black hole mass function
As it is well known, uncertainty in the local mass function determination propagates linearly in the computation of the average accretion efficiency in Soltan-type of arguments. This issue has been dealt with extensively in the literature (Soltan 1982;Marconi et al. 2004;Merloni, Rudnick & Di Matteo 2004). Recently, progress have been made towards a more secure computation of the local mass function of supermassive black holes, based on the various scaling relations. For example, Graham and Driver (2007) have shown how most previous inconsistencies among local black hole mass densities, ρBH,0, computed using different scaling relations were due to incorrect accountancy for the Hubble constant dependence. Here, we will not dwell into a lengthy discussion about the available constraints on the local distribution of SMBH masses, for which we refer the reader to the comprehensive treatment presented in Shankar et al. (2008b), who have carried out a detailed comparison among SMBH mass functions derived with different methods. Rather, we briefly describe in the following our choice of the initial condition for eq. (2). Although uncertainties do remain, with the local black hole mass density value estimated within the range ρBH,0 = (3.2 − 5.4) × 10 5 M⊙ Mpc −3 (for H0 = 70 km s −1 Mpc −1 ), the agreement between various measurements is encouraging. We will adopt a specific analytic expression for the SMBH mass function, coincident with the central value within the uncertainty range in Shankar et al. (2008b). Our adopted mass function is computed as the convolution of a Schechter function ΨM with a Gaussian scatter of 0.3 dex, to account for the intrinsic scatter of the scaling relations from which the mass function is derived. For the Schechter function, we adopt the following parametrization with parameters: φ * = 10 −3 , log M * = 8.4 and α = −1.19. With this choice, the local black hole mass density is ρBH,0 = 4.3 × 10 5 M⊙ Mpc −3 and the total number density of black holes above 5 × 10 5 solar masses is about nBH ≃ 1.3 × 10 −2 Mpc −3 .

Figure 2.
Redshift evolution of the accretion rate density onto supermassive black holes Ψ BHAR (z), multiplied by a constant (equal to 10 3 ) is plotted as a shaded area alongside the best fit to a large compilation of data for the star formation rate (SFR) density (blue dashed, from Hopkins and Beacom 2006). The shaded area represents the uncertainty in the black hole accretion rate density deriving from the observational uncertainty on the AGN luminosity functions.

Bolometric corrections and luminosity functions
The cosmic evolution of the SMBH accretion rate and its associated mass density can be calculated from the bolometric luminosity function of AGN: φ(L bol , z), where L bol = ǫ radṀ c 2 is the (intrinsic, i.e. before absorption) bolometric luminosity produced by a SMBH accreting at a rate ofṀ with a radiative efficiency ǫ rad . In practice, the accreting black hole population is always selected through observations in specific wavebands. Crucial is therefore the knowledge of two factors: the completeness of any specific AGN survey, and the bolometric correction needed in order to estimate L bol from the observed luminosity in any specific band.
As for the accuracy of our knowledge of the bolometric correction, we refer the reader to the studies of Marconi et al. (2004); Richards et al. (2006); Hopkins et al. (2007b) and Shankar et al. (2008b). All of them consistently demonstrate that a luminosity dependent bolometric correction is required in order to match type I (unabsorbed) AGN luminosity functions obtained by selecting objects in different bands 2 . No significant trend with redshift is instead evident, similarly to what is known for the spectral properties of type I AGN in the X-ray band alone (Vignali, Brandt & Schneider 2003;Tozzi et al. 2006;Steffen et al. 2006). Here we will adopt, as fiducial inputs, the bolometric corrections of Marconi et al. (2004) and the latest observed 2-10 keV luminosity function of Silverman et al. (2008), which provides the best available constraint on X-ray selected AGN at high redshift. In order to account for the measurement uncertainties in a simple and straightforward way, we will use two different analytical parametrization of the observed 2-10 keV luminosity function, as discussed in Silverman et al. (2008): a luminosity dependent density evolution model (LDDE; Ueda et al. 2003;) and a modified-pure luminosity evolution (MPLE; Hopkins et al. 2007b). The analytic expressions for these models can be found in eqs. (5-11) and in Table 4 of Silverman et al. (2008).

Accounting for obscured sources
Hard X-ray selected samples of AGN represent the most unbiased census of accretion activity in the Universe, due to the limited effect that (Compton thin) gas obscuration has on the emerging luminosity in this band. However, observed 2-10 keV luminosities of AGN are still affected by absorption, and we need to quantify this effect when looking for the intrinsic bolometric luminosity.
Assuming an average X-ray spectral model of the Gilli et al. (2007) X-ray background synthesis model, with fixed power-law slope of Γ = 1.9, we can compute the ratio of emerging to intrinsic 2-10 keV luminosity, g ≡ L X,obs /LX, as a function of the absorbing column density of cold gas, NH (in cm −2 ). Using XSPEC (Arnaud 1996) we have computed the effect of various obscuring columns on the emerging spectrum and then fitted the observed relation between g and NH with a fourth order polynomial in NH = log NH −21, and found the following analytical expression, accurate to better than 20% in the whole range 10 21 < NH < 10 24 : For Compton thick objects (log NH > 24), we have simply assumed that only a scattered component of 0.02 times the intrinsic flux is detectable (see Gilli et al. 2007).
Then, for any given observed 2-10 keV luminosity function φ ′ (LX), we can deconvolve the effect of obscuration if the distribution of column densities f (NH, LX) as a function of luminosity is known. It is now well established that the fraction of absorbed sources depends on the intrinsic X-ray luminosity of an AGN (Ueda et al. 2003;La Franca et al. 2005;Treister & Urry 2005;Barger et al. 2005;Hasinger 2008). Following again the latest incarnation of the XRB synthesis model of Gilli et al. (2007), we can write the ratio of absorbed to unabsorbed sources, R(LX) as with Lc = 10 43.5 , so that the ratio of obscured to unobscured sources tends to unity at high (QSO) luminosities and to 3.7 at low luminosities.
The distribution of absorbing columns is taken to be an increasing function of the column density for Compton thin objects, and, for the lack of better observational constraints, flat for Compton thick sources, once again fully consistent with the XRB synthesis model of Gilli et al. (2007). We thus define where A is fixed by the overall normalization Finally, we obtain the following expression for the intrinsic 2-10 keV AGN luminosity function as a function of the observed (un-corrected) one: (10) and this is what we use in our calculation of ΨC (section 3).

Do all SMBH have an "active" counterpart?
Traditionally, the presence of an accreting supermassive black hole in the nucleus of a galaxy have been associated to a number of unambiguous signs of activity: high excitation broad and/or narrow emission lines, hard X-ray luminosity exceeding 10 42 erg s −1 , the presence of relativistic jets, etc. Within this view, it is commonly held belief that the fraction of AGN in any randomly sampled population of galaxies is relatively small. This is in obvious contrast with the idea that SMBH are ubiquitous in the nuclei of galaxies (at least in the nearby Universe), and have prompted many authors to search for the "trigger" that activate an otherwise quiescent black hole. However, it must be obvious that every observational definition of AGN, like those mentioned above, suffers from strong selection biases, in particular at low intrinsic luminosities. It is in fact well known that high-resolution, highsensitivity studies of nearby galaxies not considered to be powerful AGN reveals that a very high fraction of them do indeed harbour SMBH showing clear signs of accretion. For example, Nagar et al. (2005) and Filho et al. (2006) have used high-resolution radio imaging of a magnitude limited sample of bright nearby galaxies (the Palomar sample) to show that at least a quarter (and possibly as much as half) of them are true AGN, albeit of very low luminosity. Similar results have been obtained by optical spectroscopic studies (Ho, Filippenko, & Sargent 1997) and by sensitive Chandra X-ray observations (Ho et al. 2001;Gallo et al. 2008;Santra, Sanders & Fabian 2007).
Furthermore, we notice that at z = 0, the low-luminosity slopes of both radio (Nagar, Falcke & Wilson 2005) and X-ray (Silverman et al. 2008) luminosity functions of AGN are steeper than the low mass end slope of the SMBH mass function. Therefore, if we limit ourselves to the study of SMBH with (relic) masses larger than a given value, only a minor extrapolation of the observed luminosity functions towards low luminosity is sufficient to reach total number densities comparable to that of the local SMBH.
For example, for our adopted luminosity functions, the number of SMBH with M > 5 × 10 5 M⊙ is equal to the number of objects with intrinsic (i.e. un-absorbed) 2-10 keV luminosity larger than a few ×10 40 erg s −1 and with intrinsic (i.e. debeamed, see § 6.1) 5 GHz radio core luminosity greater that a few ×10 35 erg s −1 . From hard X-ray surveys of the local Universe (Sazonov et al. 2007), as well as from XRB synthesis models (Gilli, Comastri & Hasinger 2007), we know that at these very low luminosities the number of X-ray AGN is dominated (in a ratio of approximately 4:1) by heavily obscured objects, so we expect that sensitive Chandra surveys reaching down to these luminosities should unveil nuclear accreting black holes in a large fraction of nearby galaxies, as indeed found by Gallo et al. (2008). On the other hand, radio surveys are insensitive to dust obscuration (although affected by relativistic beaming), and indeed, based on the works of Nagar et al. (2005) and Filho et al. (2006) it is reasonable to expect that essentially all nearby galaxies will show compact, high brightness temperature cores with L5GHz > 10 35 erg s −1 .
In our view, this justifies the assumption, implicit in all our subsequent calculations, that there is no distinction between SMBH and AGN, in the sense that every supermassive black hole is active, at some level, like the SMBH at the center of our galaxy clearly demonstrates. Thus, the SMBH population coincides with that of all actively accreting black holes, once the broad distribution of accretion rates is accounted for. This is a very important point, at variance with most current analytic, semi-analytic and numerical models of SMBH growth that often assume a fixed Eddington ratio, and the consequences thereof will become manifest more than once in the remain of the paper. In particular, we will discuss the observational implications for the determination of AGN fraction in different surveys in section 5.4.

UNVEILING THE GROWTH OF SUPERMASSIVE BLACK HOLES
In this section, we concentrate on the growth of the supermassive black hole population. We will proceed from the most general to the more detailed results, namely from the study of integrated quantities (SMBH mass and accretion rate density, radiative energy density, average radiative efficiency, etc.) and their redshift evolution, to differential ones, focusing on the full mass and accretion rate functions, and their evolution. The first part ( § 5.1) will shed light on some very general properties of the cosmological growth of SMBH, while the second ( § 5.2-5.4) will reveal interesting features of their so-called anti-hierarchical behaviour.
In this section, as well as in the next one, all results will be shown accounting for the intrinsic uncertainties of the adopted luminosity functions. We estimate that these uncertainties can be evaluated by comparing different analytic parametrization of the same data sets (the LDDE and MPLE parametrization for the hard X-ray luminosity function of Silverman et al. 2008, § 4.2 and the two alternative parametrizations for the flat-spectrum radio luminosity function of Dunlop and Peacock (1990) and De Zotti et al. (2005), see § 6). Then, we will solve the continuity equation for the SMBH evolution according to the method described in section 3 for the four possible combinations of X-ray and radio luminosity functions, and present the results as shaded areas encompassing the entire range spanned by these four solutions.

Integral constraints, mass density and radiative efficiency
The redshift evolution of the black holes accretion rate (BHAR) density can be easily calculated from the intrinsic X-ray luminosity function (10) as follows: where LX is the intrinsic X-ray luminosity in the rest-frame 2-10 keV band, and the bolometric correction function L bol (LX) is given by eq. (21) of Marconi et al. (2004).
In Figure 2 we show the redshift evolution of ΨBHAR(z), multiplied by a constant (equal to 10 3 ) in order to visually compare its shape with that of the star formation rate (SFR) density, which we take from the best fit to a large compilation of data from Hopkins and Beacom (2006). The evolution of the BHAR density of Fig. (2) has been computed for ǫ rad ≃ 0.07, but changing this value would simply result in a change of the overall normalization.
By construction, the computed SMBH accretion rate density includes both absorbed and unabsorbed sources as well as Compton thick AGN ( § 4.3), based on the prescription of Gilli et al. (2007) that uses, as almost unique observational constraint, the shape and normalization of the XRB radiation, and is thus somewhat degenerate in the luminosity-redshift plane (see, however, Worsley et al. 2005). As a sanity check, we can compute the total contribution of Compton thick sources to the total mass accumulated in SMBH today, which turns out to be limited to less than about 20%. This Bottom panels, the evolution of the average Eddington ratio of the SMBH population, defined in eq. (13). On the left we show the results of a calculation performed fixing the accretion efficiency to η = 0.1, which, for the particular model of accretion modes we assume correspond to an average radiative efficiency of ǫ rad ≃ 0.086, while on the right we show the case of η = 0.083, corresponding to ǫ rad ≃ 0.07. See text for details.
confirms that our results are not too sensitive to the uncertainties in the luminosity and redshift distribution of Compton thick AGN.
The integrated SMBH mass density as a function of redshift can then be calculated starting from a local value for ρBH,0 (see section 4.1): In this way, for any given φ(LX, z) and bolometric correction (thus ΨBHAR(z)), the exact shape of ρBH(z) then depends only on two numbers: the local black holes mass density ρBH,0 and the (average) radiative efficiency ǫ rad . This fact has an interesting, and somewhat surprising, consequence: even before making any assumption about the distribution of the Eddington ratio in the AGN population, we are able to constrain the (mass weighted) average radiative efficiency, in the following way. A firm lower limit to ǫ rad can be established using the Soltan argument. Specifically, it is enough to look for the value of the radiative efficiency that makes the integral z 0 Ψ BHAR (z ′ ) ρ BH,0 dt dz ′ dz ′ of eq. (12) larger than unity. For our choice of the local SMBH mass function and AGN bolometric luminosity function this lower limit is approximately equal to 0.065.
However, when starting from the local mass density to compute its evolution, the shape of ρBH(z) turns out to be very sensitive to the average radiative efficiency. The very fact that, in much the same way as the SFR density, also the black hole accretion rate density, ΨBHAR(z) has a well pronounced maximum (around z ∼ 2), implies that, if the adopted value of ǫ rad were too high, the black hole mass density would not decrease any longer towards higher redshift. This is clearly shown in the upper two panels of Figure 3, where we show the evolution of the total SMBH mass density com-puted for two different values of the accretion efficiency (and, correspondingly, of the radiative efficiency): η = 0.1 (corresponding to ǫ rad ≃ 0.086) and η = 0.083 (ǫ rad ≃ 0.07).
The same issue can be considered from a different point of view. Let us now define a global, mass weighted average Eddington ratio for the entire SMBH population as: The redshift evolution of Λ(z) is shown in the two bottom panels of Fig. 3 3 . There it is apparent that, while for the low radiative efficiency case (lower right panel) the average Eddington ratio is allowed to increase up to the highest redshift for which reliable X-ray luminosity data are available (z ∼ 4), increasing the efficiencies by less than 20% modifies the evolution of the Λ substantially (see bottom left panel), being it almost flat, and possibly decreasing at z > ∼ 3, an obvious consequence of the flattening of the corresponding mass density (top left panel). The physical explanation for this is straightforward: if we calculate the evolution of the SMBH mass density with a too high radiative efficiency, we fail to account for the mass locked up in SMBH in the local Universe, and we are forced to explain a "primordial" population of very weakly active (or even dormant) black holes, of substantial density in place by z ∼ 4 (in the example of Fig. 3, left panel, amounting to about 10 5 M⊙ Mpc −3 ). This puts already some strain on existing models for the black holes seed population (see e.g. Trenti & Stiavelli 2006;Volonteri, Lodato and Natarayan 2007). Models in which the primordial black holes originate from stellar mass progenitors (remnants of the first, PopIII, stars; Abel, Bryan & Norman 2000;Bromm, Coppi & Larson 2002), predict negligible BH mass density at high redshift (say, less than 10 4 M⊙ Mpc −3 at z > 5), and are incompatible with an average, mass weighted radiative efficiency larger than about 0.075 (in very good agreement with the recent estimates of Shankar et al. 2008b), or, correspondingly, with an accretion efficiency η > ∼ 0.09. Such a tight upper limit on the radiative efficiency can be somewhat relaxed if: (a) the local SMBH mass density has been overestimated (see the discussion in Merloni et al. 2004); (b) the Compton Thick population has a very different redshift distribution from that assumed in the Gilli et al. (2007) XRB synthesis models (with more heavily obscured high luminosity objects at high redshift); (c) the primordial BH seeds are indeed very massive, resulting from the direct collapse of supermassive stars (Koushiappas, Bullock & Dekel 2004;Begelman, Volonteri & Rees 2006;Lodato & Natarayan 2007).
On the other hand, the strict lower limit on the efficiency could also be relaxed if the local SMBH mass density had been underestimated, or if there was a substantial population of "wandering" SMBH, i.e. objects that have been kicked out of their galactic host, most likely for the recoil experienced due to asymmetric gravitational wave emission subsequent to a binary black hole merger event, after having accreted substantial amount of matter. Those black holes would then be computed in the AGN luminosity functions, but not among the local relic population that contributes to ρBH,0. A discussion of the relevance of this effect on the mass function of SMBH can be found in Volonteri (2007).
Summarizing, we can define: (i) a normalized local SMBH mass density ξ0 ≡ ρBH,0/(4.3 × 10 5 M⊙Mpc −3 ); (ii) the normalized SMBH mass density at z = zi (where zi is our current observational horizon, i.e. the maximum redshift at which AGN luminosity functions are well known) ξi ≡ ρBH(z = zi)/ρBH,0 and (iii) the normalized mass density of SMBH ejected from galactic nuclei due to the gravitational wave recoil after mergers, ξ lost ≡ ρ BH,lost /ρBH,0. Then, our calculations show that the average radiative efficiency is constrained to be Finally, we can calculate explicitly how much mass have SMBH accumulated in different modes. Specifically, we find that between about 73% and 80% (depending on the choice of the luminosity functions) of the relic black hole mass density has been accumulated in radiatively efficient modes (either HR or HK), while only 20% to 27% in the low radiatively efficient (LK) mode, even if most of the time in the life of a SMBH is likely spent in this latter mode (Cao 2007;Hopkins, Narayan & Hernquist 2006). The accretion rate density evolution for both the entire population and for the SMBH in the LK mode is shown in Figure 4.

The SMBH mass function evolution
So far, we have only discussed integral constraints, i.e. those obtained from the study of the evolution of the integral of the mass (or luminosity) function. However, the solution of eq. (2) provides us with the full evolution of the SMBH mass function, and we proceed now to discuss it in detail.
The SMBH mass function evolution is shown in Figure 5 for our fiducial case of ǫ rad = 0.07. In the different panels, each corresponding to a different redshift, the dashed lines show the low redshift (z = 0.1) mass function as a reference, while the shaded areas represent the calculated mass functions with their own uncertainties. The so-called "downsizing" is apparent in the first four or five panels (i.e. up to redshift ∼ 1.5 − 2): we clearly see how in this time interval the mass function grows much more substantially at the low-, rather than at the high-mass end. This is to be expected, in the light of the results of Heckman et al. (2004). There, on the basis of a sample of about 23,000 SDSS type II AGN for which both an estimate of the black hole mass (via a measurement of the velocity dispersion of the host) and of the intrinsic accretion rate (via measures of the narrow line luminosity) were available, it was shown that it is only the small mass black holes that are growing rapidly at z ∼ 0.2. Similar trends, indicating higher Eddington ratios (and smaller growth times, see below) have been also found recently in the SDSS quasar sample (Netzer & Trakhtenbrot 2007), dominated by un-obscured, luminous AGN at higher redshift. We stress, however, that our results, based on a much more statistically complete census of the AGN population (thanks to the hard X-ray selection), is much less prone than previous ones to selection biases.
Differently from what previously found in M04, we note here that by the epoch of the peak of the AGN/QSO activity (z ∼ 2), the high-end of the mass function decreases significantly, too, in coincidence with the observed drop of the number density of bright QSOs at those redshifts (Fan et al. 2001;Wall et al. 2005;Hasinger, Miyaji & Schmidt 2005;Richards et al. 2006;Silverman et al. 2008;Brusa et al. 2008).

The accretion rate distribution function: downsizing and its reversal
As we have stressed a few times already, our calculations do not assume any particular distribution of accretion rates for AGN, rather derive it consistently from the observed mass and luminosity functions (see § 3). It is thus interesting to show the evolution of the (Eddington-scaled) accretion rate function, here defined as: where Φ M,Ṁ (M,Ṁ ) is given by eq. (5), the relationship betweeṅ M and λ is simply given by λ = ǫ radṀ c 2 /L Edd , and λ is related to the accretion rate via eq. (1). We show the Eddington ratio distribution function evolution in Figure 6. The various panels make evident the change in the shape of the distribution, which results from the combined evolution of the X-ray luminosity function (a form of luminosity dependent density evolution is noticeable in the evolution of Φ λ as well) and of the mass function itself. Overall, the trend is for a progressive flattening of the accretion rate distribution function with increasing redshift, i.e. of a an increasing relative importance of highly accreting objects, as it is reasonable to expect on physical grounds if the SMBH population grows from small mass high-redshift seeds. The evolution of the Eddington ratio distribution function shown in Figure 6 does not contain any information about the typical accretion rate as a function of black hole mass. A useful way to present this information is that of showing the specific growth rate as a function of SMBH mass. We have computed this quantity by taking the ratio of the black hole mass to the average accretion rate Ṁ (M ) . Such a ratio M/ Ṁ defines a timescale, the socalled growth time, or mass doubling time, as it measures the time it would take a black hole of mass M to double its mass if accreting at the currently measured rate ofṀ . The redshift evolution of the growth time (in years) as a function of black hole mass is shown in Figure 7. As a reference, we show in each panel the age of the Universe as a horizontal dashed line. Black holes with growth times longer than the age of the Universe are not experiencing a major growth phase, which must have necessarily happened in the past. On the contrary, objects with growth times shorter than that are actively growing.
The z = 0.1 growth time distribution is in good agreement, both in shape and normalization, with that first measured in the above mentioned work of Heckman et al. (2004). The positive slope of the average growth time vs. M is one of the most direct evidence of a SMBH "downsizing", i.e. of the fact that in the local universe only small mass black holes are actively growing. Interestingly enough, a similar trend has been since unveiled for star forming galaxies (Feulner et al. 2005;Noeske et al. 2007;Pérez-González et al. 2008), where the ratio of galaxy mass M gal to star formation rate is studied as a function of M gal . It is clear that a direct comparison of the SMBH growth times so measured with those of the stars in galaxies of various morphologies/ages/masses may hold the key for a thorough understanding of the SMBHgalaxies co-evolution, and we aim at presenting such a comparison in a following paper.
The redshift evolution of the growth times distribution can be used to identify the epochs when black holes of different masses grow the larger fraction of their mass. Figure 7 shows that as we approach the peak of the black hole accretion rate density (z ∼ 1.5 − 2, see fig. 2), we witness the rapid growth of high mass black holes, too. The growth time distribution as a function of mass becomes flat there, and sinks below the corresponding age of the Universe. This is the period of rapid, widespread growth of the entire SMBH population. Correspondingly, above this redshift, the "downsizing" trend seems to disappear.
Finally, we show in Figure 8 the redshift evolution of the radiated (bolometric) energy density split into various black hole mass bins. On the left, we divide objects according to their mass at each redshift. The low redshift downsizing trend is here apparent in the rapid decline of the radiated energy density of the most massive black holes (log M > 9), which becomes less dramatic when less massive holes are considered. While in the local Universe (i.e. for z < ∼ 0.3) the larger contributors to the total radiative energy density are the numerous, small, rapidly accreting SMBH (with log M < 7), above redshift one most of the radiative energy is produced by larger black holes, typically with masses between 10 8 and 10 9 solar masses. As we move to higher redshift, the numbers of black holes with log M > 8 declines, and the radiative energy density is dominated again by smaller black holes with 8 > log M > 7. The most massive black holes, those with log M > 9, never dominate the radiative energy density injection, essentially because of their small numbers, being always in the exponentially declining part of the mass function, well above its knee (see Figure 5). As we discussed in the Introduction, cosmological studies of SMBH evolution offer the advantage, with respect to those of galaxies, that SMBH do not undergo any morphological transformation as they grow. Thus, once enough information is gathered to solve the continuity equation (2), we can trace back the evolution of any black hole (or of any population of black holes of any given mass at a particular redshift). Indeed, it was by comparing the growth histories of black holes of different masses today that Marconi et al. (2004) first noticed what they called the "antihierarchical" evolution of the SMBH population (see also M04; Shankar et al. 2008b). Obviously, we can do the same here. On the right hand side of Figure 8, we show the redshift evolution of the radiative (bolometric) energy density split into bins of different black hole masses today, i.e. at z = 0.
From this we see that the progenitors of local SMBH with masses in the range 8 < log M < 9 dominate the AGN radiative energy output between z ∼ 0.3 and z ∼ 3, when the progenitors of the largest black holes today (those with log M > 9) started to become at least equally important.

Selection effects: mass function of "active" black holes
From our discussion so far, it should be clear that a complete picture of the growth and evolution of the SMBH population cannot be achieved without appreciating the full extent of the wide distribution of both masses and accretion rates, a point often overlooked in either semi-analytic or numerical models. As a consequence, the very concept of a clear distinction between "active" and "inactive" black holes can be misleading, once we consider within the same scheme black holes accreting close to the Eddington limit and those with an Eddington-scaled accretion rate as small as, say, 10 −4 .
Thus, we advocate here more operational definitions of "active" black hole, the easiest and most straightforward of which are those based on nuclear flux limits. As an illustrative example, we plot in Figure 9 the mass functions of SMBH whose observed Xray flux exceeds some fixed limits. In particular, alongside the total mass function, reproduced here from Fig. 5 as a term of comparison, we plot the mass functions of all SMBH with an observed (i.e. absorbed) 2-10 keV flux above 10 −13 (red), 10 −14.3 (blue) and 10 −15.35 (green) erg/s/cm 2 . The latter two limits have been chosen to match the observational limits of recent deep (Chandra Deep Field South, CDFS; Giacconi et al. 2002) and medium-deep (XMM-COSMOS, XCosmos; Hasinger et al. ) surveys. Due to the relative steepness of the accretion rate functions at low redshifts, the mass functions of active black holes so defined are extremely sensitive to the flux limit defining what "active" means.
This fact implies that, when trying to determine the SMBH mass function from a survey that selects AGN above a given flux limit, one should carefully account for the bias introduced by the fact that SMBH of any mass always have a broad distribution of accretion rates (Babić et al. 2007). For obvious reasons of simplic- ity and coherence, we discuss here hard X-ray flux limits, as we use the hard X-ray luminosity function to unveil the SMBH evolution, and the X-ray background to constrain the relative numbers of obscured and unobscured sources. However, it is clear that, provided robust bolometric corrections are available, similar calculations can be converted into flux limits in any band. For example, Greene and Ho (2007) have recently presented the mass function of active black holes, identified as broad Hα line emitters (in the redshift range 0.1 < ∼ z < ∼ 0.3) in the SDSS database. Comparison with the soft X-ray selected AGN luminosity function of Hasinger et al. (2005) showed that their sample approach completeness at luminosity corresponding to a (rough) X-ray flux limit of 10 −13 erg/s/cm 2 . Indeed, the mass function we show at low z in Fig. 9 for such a flux limit is in quite good agreement with that reported in Greene and Ho (2007), once the fact that our active mass functions do include both obscured and unobscured objects is accounted for.
In Figure 10 we show the ratio of the various "active" mass functions to the total SMBH mass function. These curves can be simply interpreted as "selection functions" for the black hole mass of AGN surveys of corresponding hard X-rays flux limits. However, as far as we can claim that, at all redshifts considered here, there is a one-to-one correspondence between SMBH and galaxies mass (as it was implicitly assumed to derive the z = 0 mass function), then the various curves of Figure 10 could be also interpreted as "AGN fractions" as a function of SMBH (and galaxy) mass. Once again, due to the main tenet of our method of calculation, namely the fact that growing supermassive black holes display a very broad accretion rate distribution at any given mass, the very concept of "AGN fraction" is strongly dependent on the sensitivity limit that defines what an AGN is. In general, the AGN fraction tends to be highest at . The right hand plot shows the same quantities for objects divided according to their instantaneous mass, while the right hand one is for objected partitioned according to their final (z = 0) mass. the high mass end, and declines at low masses. A deep X-ray survey like the CDFS has a quite flat 'mass sensitivity' (or mass fraction) at z ∼ 1, while the shallower XMM-COSMOS at the same redshift starts declining already for log M < 8. At higher redshift even the deepest surveys are highly incomplete (in a SMBH mass sense), and the incompleteness rises steeply with decreasing mass. In the case of spectroscopic studies (see e.g. Decarli et al. 2007), where the definition of an active galactic nucleus relies on the capability of spatially and spectrally resolve AGN-excited emission lines, the AGN fraction will be a more complicated function of SMBH (and galaxy) mass than shown in Figure 10.

THE EVOLVING KINETIC LUMINOSITY FUNCTION OF AGN
So far, we have focused our attention on the growth of the SMBH population under the assumption that, whatever the mode of accretion of a black hole, it s possible to establish a unique correspondence between radiated (bolometric) luminosity and rate of change of black hole mass (see section 2). A choice of local black hole mass function, an expression for the bolometric correction, and a choice of a value for the accretion efficiency was all that was needed to extract the useful information on the SMBH growth from the observed X-ray luminosity function. However, the growth of supermassive black holes through mass accretion is accompanied by the release of enormous amounts of energy which, if it is not advected directly into the hole (see e.g. Narayan & Yi 1995), can not only be radiated away, but also disposed of in kinetic form through powerful, collimated outflows or jets, as observed, for example, in Radio Galaxies. To reiterate the point made before ( § 2), it is only by assessing the relative importance of the radiative and kinetic energy output in either "low kinetic" (LK), "high radiative" (HR) and "high kinetic" (HK) modes Figure 10. Ratio of SMBH mass function above a given flux limit to the total mass function. The various shaded areas correspond to different observed 2-10 keV flux limits, the same as those used in figure 9, i.e. 10 −13 (red), 10 −14.3 (blue) and 10 −15.35 (green) erg/s/cm 2 . that a full picture of the effects of SMBH growth on the environment will possibly emerge. This is important, as we expect that radiative and kinetic feedback differ not only on physical grounds, in terms of coupling efficiency with the ambient gas, but also on evolutionary grounds. The anti-hierarchical growth of the SMBH population discussed in sections 5.2-5.4 necessarily implies a delay between the epoch of Quasar (HR) dominance and that of more sedate LK mode, such that kinetic energy feedback plays an increasingly important role in the epoch of cluster formation and virialization.
In this section, we present a method to complement the information on the cosmological growth of the SMBH population described in the previous sections with an analysis of the evolving kinetic luminosity function of AGN. As in Merloni (2004), we will make use of the radio luminosity function of AGN as the main observable for our study. However, differently from the choice made there, we will focus our attention on the intrinsic radio core luminosity function. The main reason for this choice is twofold.
On the one hand, the observed relation between BH mass, radio and X-ray luminosity (the fundamental plane of active black holes, MHD03) that defines the physical state of LK mode objects is based on the observed (5 GHz) radio core emission, not on the extended one. Also, Merloni and Heinz (2007) have shown that the intrinsic radio core luminosity of an AGN jet is a good indicator of its total kinetic power, once the effects of relativistic beaming are taken into account. On the other hand, the physical model for the radiative emission in the core of a relativistic jet is relatively well established (Blandford & Koenigl 1979;Heinz & Sunyaev 2003), as opposed to the complicated physics of radio emission and particle acceleration in the large scale lobes of a radio jet, where environmental effects do necessarily play an important role. Obviously, the price to pay for the use of radio core emission as a tracer of the jet power, is that relativistic beaming may severely affect the derived parameters, and must be taken into account. In the classical model by Blandford & Koenigl (1979), the flat spectrum radio synchrotron emission of a compact jet core is produced by a superposition of self-absorbed synchrotron spectra, each from a different region in the jet. The model predicts a dependence of the monochromatic radio luminosity Lν on jet kinetic power W kin of the form Lν ∝ W 17/12 kin . More generally, Heinz & Sunyaev (2003) showed that any scale invariant jet model producing a power-law synchrotron spectrum with index αν must obey the relation Lν ∝ W (17+8αν )/12 kin M −αν . By studying a sample of nearby low-luminosity radio AGN for which the kinetic power had been measured,  have shown that, once beaming is statistically taken into account, a tight relation between the 5 GHz intrinsic radio core luminosity (LR) and the jet kinetic power holds, in the form: where W0 ≃ 1.6 × 10 36 and L0 = 10 30 erg s −1 , and the intrinsic scatter is about 0.35 dex. In §6.1, we will use this relation to construct the kinetic jet luminosity function from the observed flat spectrum radio luminosity function ΦL(Lν ) (abbreviated as FSLF below). Similar, but complementary studies of the AGN kinetic luminosity function evolution based on steep spectrum sources has been recently carried out (Körding, Jester & Fender 2008;Shankar et al. 2008a). The results of both these works are in reasonable agreement with those we present here.

Beaming effects and the intrinsic luminosity function of radio jet cores
As pointed out in Heinz, Merloni & Schwab (2007), we can use eq. (16) and the observed FSLF to derive the underlying kinetic luminosity function ΦW of radio jet cores: Without loss of generality, and consistent with most observational studies, we will use a broken power-law to describe φR(LR): In the above equation (16), however, it is the intrinsic radio core luminosity that enters the determination of the kinetic power. Since jets are relativistic, the observed luminosity functions, φ ′ are always affected by Doppler boosting. Assuming all object have the same bulk Lorentz factor γ = (1 − β 2 ) −1/2 , the observed luminosity L of a relativistic jet is related to the emitted one, L, via the relation L = δ p L, where δ, the kinematic Doppler factor of the jet is defined by δ = [γ(1 − βcosθ)] −1 . The exponent p depends on the assumption about the jet structure, and is here taken to be p = 2 as appropriate for continuous jets, rather than discrete ejections (for which p = 3, see Urry & Padovani 1995). If the radio sources have all a pair of oppositely directed jets, and the jets are randomly oriented in the sky, the cosine of the angle between the velocity vector and the line of sight, cos θ have a uniform probability distribution, and the observed sources will be those lying at angles between θmin and θmax, corresponding to Doppler factors of δmin and δmax, respectively (see Urry and Shafer 1984). We assume here that θmax = 0 • and δmin = 1/γ (i.e. neglecting all strongly de-beamed objects). Urry and Padovani (1991) have shown how, under the above assumptions, an intrinsic double power-law LF is transformed by the effects of relativistic beaming. It can be shown that, for astrophysically relevant values of a1, a2 and γ the most general situation is that in which the beamed LF is a triple broken powerlaw (Heinz, Merloni & Schwab 2007). The faint-and bright-end slopes are unchanged, while around the knee an intermediate-slope −(p + 1)/p = −3/2 is found.
The problem of inverting the Urry and Padovani (1991) solution, i.e. of deriving the intrinsic shape of the radio loud AGN LF given the observed one and a few assumption about p and γ is greatly simplified if the intermediate part of the beamed luminosity function is neglected. In practice, in such a case the intrinsic, debeamed, FSLF φ can be obtained from the observed one simply by the following transformations of the LF knee, while keeping both end slopes fixed. In detail, following the procedure described in Urry and Padovani (1991), the LF normalization will be corrected by a factor while the break luminosity transforms as: where Figure 11. Kinetic luminosity function of AGN jets (here plotted as W kin × φ W ) at four different redshifts. In each panel the green area between the black solid lines represents the total kinetic luminosity function, with the uncertainty reflecting the difference in the luminosity functions between Dunlop & Peacock (1990) and De Zotti et al. (2005). The blue area between blue dashed lines represents the contribution from sources in the LK mode, while the dotted lines mark the contribution from sources in the HK state. In the Left plot (FS) the jet core luminosity function includes only flat spectrum sources, while in the right plot (FS+SSmax) a maximal contribution from hidden cores in steep spectrum sources is included (see text for details).
In this way, for any given luminosity function and mean Lorentz factor, it is possible to derive the intrinsic kinetic luminosity function of the AGN jet, by simply applying eq. (17) to the de-beamed radio luminosity function.
As opposed to the case of X-ray selected luminosity functions, the available constraints on the radio cores luminosity functions are not very stringent. Here we will use two different determinations, translated into our Cosmology, and computed at 5GHz: the flat spectrum radio luminosity function from Dunlop & Peacock (1990) (their HIGH-z parametrization of their table C4), and that from De Zotti et al. (2005), given by the sum of flat spectrum radio quasar and BL Lac objects as in eqs. (1)-(3) and table 1 of their paper. As we will show in the following, and as already discussed in Heinz, Merloni and Schwab (2007), the total kinetic power released by AGN is sensitive to the faint-end slope of the kinetic luminosity function. For consistency with observations of local radio luminosity functions down to very low luminosities (Nagar, Falcke & Wilson 2005), we have fixed the z = 0 faint end slope to a1 = 1.85 for both the FSLF of Dunlop & Peacock (1990) and the BL Lac population of De Zotti et al. (2005). Because a1 −1 is (marginally) larger than 0.81, the slope of the W kin -LR relation (16), the total kinetic power is determined by the lower cut-off of our radio luminosity function, chosen here to correspond to the radio luminosity above which the total number of sources is equal to the total number of SMBH with mass above 5 × 10 5 M⊙ (see § 4.4). However, the difference between the total kinetic power obtained in this way and that one would get for a1 < 1.81 is smaller than the uncertainty introduced by the intrinsic scatter in the relation (16). Therefore, for the ease of computation, we have allowed a very mild evolution (flattening) of such a faint-end slope between z = 0 and z = 2, where we have fixed a1 = 1.75. Indeed, there are independent lines of evidence pointing towards a flattening of the faint-end of the AGN radio luminosity function at high redshift (Cirasuolo, Magliocchetti & Celotti 2005), similarly to the well known evolution observed in the X-ray band. Future better constraints on the high redshift evolution of the faint radio AGN population will surely help tightening our constraints on the overall kinetic power density. The observed luminosity functions are then de-beamed following the above described prescription, assuming p = 2. For the mean jet Lorentz factor, we adopt the most likely value based on the statistical analysis of core radio luminosity of a small sample of 15 AGN with measured kinetic power (see Merloni and Heinz 2007 for details) 4 : γ = 7. Also, for the sake of simplicity, and lacking any clear observational evidence, we do not account here for the possibility of varying Lorentz factors with redshift.
So far, we have proceeded under the assumption that the observed FSLF are complete, in the sense that they account for the totality of jet cores. In fact, we have to discuss also the possible contribution from flat spectrum cores which are 'hidden' below a more powerful, extended steep-spectrum component. By definition, the FSLF contains all black holes except those included in the steep-spectrum luminosity function (SSLF), which are dominated by optically thin synchrotron emission. The scaling relation from eq. (16) does not hold for extended, steep-spectrum sources (but see Bîrzan et al. 2004;Best et al. 2006). We can, however, derive an upper limit to the contribution from steep spectrum sources. Assuming a typical optically thin synchrotron spectral index of 0.65, any underlying flat spectrum component would have to fall below ∼20% of the observed 5GHz luminosity, otherwise the source would become too flat to qualify as a steep-spectrum source. Using the steep-spectrum luminosity functions (SSLF) from Dunlop & Peacock (1990; high-z parametrization of their table C4) and De Zotti et al. (2005), we then calculate a maximal core luminosity function comprising both flat spectrum sources and those 'hidden' cores for which we take a contribution equal to 10% of the SSLF. This can be considered as an extreme case in which the objects with powerful (optically thin) radio lobes have a core-to-lobe flux ratio (the so-called core prominence, Bridle et al. 1994) of 0.1 (for a discussion on the observational constraints on the intrinsic core prominence, see Laing et al. 1999, Jackson andWall 1999).
The kinetic luminosity functions calculated from the observed radio cores luminosity function with the help of eq. (16) are shown in Figure 11, where we plot the luminosity function ΦW multiplied by the kinetic power W kin in a representation that highlights the luminosity range where most of the power is released. For the sake of clarity we have chosen to keep separate the results obtained from adopting only the FSLF (left panel [FS]) from those obtained adding to the FSLF a contribution of 10% of the observed SSLF (right plot [FS+SSmax]). As discussed above, we regard this latter as an effective upper limit to the AGN kinetic luminosity function. Figure 11 also shows the decomposition of the total kinetic luminosity functions into the two kinetically dominated modes of accretion, as discussed in § 2. For both cases (FS and FS+SSmax), at z < ∼ 1, the kinetic luminosity functions is dominated by low luminosity (LK) objects, with bright, radio loud QSOs (HK) only contributing to the high power tail. At, and above, the peak of the BHAR density evolution (z ∼ 1.5 − 2.5, see Fig. 2) HK objects start dominating the total kinetic power output of the growing black holes population, even if the relative contribution of LK sources is always substantial.

Kinetic energy density evolution and kinetic efficiency
Integrating the kinetic luminosity functions obtained following the procedure described in the previous section, we arrive to an estimate of the total kinetic energy density as a function of redshift. Figure 12 shows the redshift evolution of the total integrated kinetic energy density (divided by ηc 2 , so that it can be displayed in units of M⊙ yr −1 Mpc −3 ). Once more, we show separately the results of the calculation done including flat spectrum sources only (left panel) and one in which a maximal contribution from hidden cores in steep spectrum sources is added (right panel). The total jet kinetic power density (shown as grey shaded areas) is split into the contributions from LK and HK modes of growth. Because the average Eddington ratio increases towards high redshift, so it does the relative contribution of HK objects (i.e. radio loud QSOs, or powerful FR II), irrespective of the choice of core radio luminosity function. As a term of comparison, the plot shows also the evolution of the total black hole accretion rate density (BHAR), which we discussed at length in section 5.1. Also shown is the estimated total kinetic power output of SNe II, as computed directly from the best fit to the star formation rate density evolution (Hopkins and Beacom 2006). The total kinetic power released by low luminosity (LK) AGN at z < ∼ 0.3, where it is by far the dominant contributor to the total jet kinetic power, is comparable to that released by type II supernovae, which, however, become about a factor of few (up to about an order of magnitude) higher at the peak of star formation activity.
Given the stark differences in the redshift evolution of lowand high-power sources, and the fact that the relative importance of kinetic power output increases as the Eddington ratio decreases (see § 2), we can conclude that the most important feature of the evolution of the kinetic energy output of growing black holes, as compared to the radiative power output (see Figure 8) is the very mild redshift evolution. In this respect, it is interesting to mention that Cavaliere & Lapi (2008) have recently shown how the problem of the 'missing baryons' in groups and cluster of galaxies can be understood in terms of a 'dual' model for AGN feedback (kinetic at low z and radiative at high z) that is in very good agreement with the overall picture of SMBH growth and energy output that emerges from the work presented here.
A fundamental question we are seeking the answer for is the accurate determination of the kinetic energy production efficiency of growing black holes. This bears obvious important consequences for AGN feedback models, but has so far been treated as a free parameter. Our calculations allow for a more precise determination of such parameter. Let us first compute the total integrated (mass weighted) average kinetic efficiency as: These numbers are consistent with those found in Heinz, Merloni and Schwab (2007), once we consider the mean Lorentz factor of γ = 7 used to calculate the intrinsic FSLF, and the steeper slope of the LR-W kin correlation (16). On the other hand, they are somewhat smaller (less than a factor of 2) than those found in both Körding, Jester & Fender (2008) and in Shankar et al. (2008a), where, however, only steep-spectrum sources were used, and the total kinetic power was derived from the correlation between extended radio luminosity and kinetic power of Willott et al. (1999).
Combining the results discussed in section 5.1 on the average radiative efficiency with those shown in eq. (21) we conclude that SMBH during their growth from z ∼ 5 till now convert about 15-25 times more rest mass energy into radiation than into kinetic power, with the exact number depending on the poorly known details of the intrinsic jet cores luminosity function, as well as on our assumptions about the beaming corrections to be made. We note here that numerical simulations of AGN feedback in merging galaxies  fix the feedback efficiency to about 5 × 10 −3 in order to reproduce the M − σ relation. In light of what we have shown, it could be argued that this feedback might be almost entirely provided by kinetically dominated modes of SMBH growth, as long as the coupling efficiency between the kinetic power provided by AGN and the interstellar and intra-cluster medium is equal to one.
We can also compute directly the kinetic efficiency as a function of redshift and SMBH mass, which we show in Fig. 13. In each of the two panels, corresponding to the FS and FS+SSmax cases, the horizontal solid lines show the mass weighted average from eq. (21). The various curves describe the main properties of kinetic feedback as we observe it. For each of the chosen mass ranges, the kinetic efficiency has a minimum when black holes of that mass experience their fastest growth: this is a different way of restating the conclusion that most of the growth of a SMBH happens during radiatively efficient phases of accretion. However, when the mass  fig. 4), while red (dashed lines) and blue (dot-dashed lines) bands represent the relative contribution to the total kinetic power from LK and HK sources, respectively. The green solid line shows the estimated total kinetic power output from SNe II, as inferred by Hopkins and Beacom (2006) based on the measured star formation rate density. The left plot shows the case in which only flat spectrum radio AGN are include in the intrinsic core radio luminosity function (FS), while the right hand plot shows the case in which a maximal contribution from 'hidden' flat spectrum cores in steep-spectrum sources is included. Figure 13. Redshift evolution of the kinetic efficiency ǫ kin . SMBH of different masses are here plotted separately, with a color coding analogous to that of Fig. 8. The right hand plot shows the results obtained including flat spectrum radio sources only, while the right hand one shows the results obtained including also a maximal contribution from the hidden cores of steep spectrum radio AGN. In each plot, horizontal black solid lines mark the mass weighted average values for the kinetic efficiencies, with the dashed lines representing the uncertainties from the particular choice of radio and X-ray luminosity functions.
increases, SMBH are more and more likely to enter the LK mode (see section 2), which increases their kinetic efficiency. More massive holes enter this phase earlier, and by z = 0 they have reached the highest kinetic efficiency of ∼ 2 ÷ 5 × 10 −2 , as it is clearly seen in Figure 13. This is a natural consequence of the observed anti-hierarchical growth of the SMBH population, and of the chosen physical model for the accretion mode of low-Eddington ratio objects.

DISCUSSION AND CONCLUSIONS
We have presented a comprehensive review of our current knowledge of the evolution of AGN and of the associated growth of the supermassive black holes population.
Similar to the case of X-ray background synthesis models, where accurate determinations of the XRB intensity and spectral shape, coupled with the resolution of this radiation into individual sources, allow very sensitive tests of how the AGN luminosity and obscuration evolve with redshift, we have argued that accurate determinations of the local SMBH mass density and of the AGN (bolometric) luminosity functions, coupled with accretion models that specify how the observed AGN radiation translates into a black hole growth rate, allow sensitive tests of how the SMBH population (its mass function) evolves with redshift. By analogy, we have named this exercises 'AGN synthesis modelling'. In performing it, we have taken advantage of the fact that the cosmological evolution of SMBH is markedly simpler than that of their host galaxies, as individual black hole masses can only grow with time, and SMBH do not transform into something else as they grow. Moreover, by identifying active AGN phases with phases of black holes growth, we can follow the evolution of the population by solving a simple continuity equation (2), where the mass function of SMBH at any given time can be used to predict that at any other time, provided the distribution of accretion rates as a function of black hole mass is known (see section 3).
Such a straightforward approach is fully justified as long as SMBH coalescences, which are likely to occur in the nuclei of merged galaxies, do not alter significantly the shape of the mass functions. It is well established, as we have also shown here, that the local black hole mass density is fully consistent with the mass accumulated onto SMBH in AGN episodes for reasonable values of the accretion efficiency (see § 5.1). Mergers only redistribute mass so they do not affect any integral constrain, such as the Soltan argument, and variations thereof. Until now, evaluating the relative importance of SMBH mergers in the mass function evolution has proved itself very difficult (Yu & Tremaine 2002), due to the large uncertainties in both merging rates of galaxies (either theoretically or observationally, see e.g. Bundy, Treu & Ellis 2007;Guo and White 2008;Hopkins et al. 2008 and references therein) and in the physical processes involved in the merger of the two nuclear black holes (see Colpi et al. 2008 for a recent review). Shankar et al. (2008b) have attempted to estimate such effect considering only equal mass merger, at a rate most likely far in excess of the true one, and concluded that the effect of SMBH mergers on the local mass function is overall smaller than the current uncertainties in the AGN bolometric luminosity function itself, and may be relevant, if at all, only at the very high mass end of the distribution. A more detailed comparison, taking into account the full mass and redshift dependencies of the merger rates, is however necessary before a firm conclusion can be drawn.
In order to carry out our calculation, we assumed that black holes accrete in just three distinct physical states, or "modes": a radiatively inefficient, kinetically dominated mode at low Eddington ratios (LK, the so-called "radio mode" of the recent literature), and two modes at high Eddington ratios: a purely radiative one (radio quiet, HR), and a kinetic (radio loud, HK) mode, with the former outnumbering the latter by about a factor of 10. Such a classification is based on our current knowledge of state transitions in stellar mass black hole X-ray binaries as well as on a substantial body of works on scaling relations in nearby SMBH (see § 2). It allows a relatively simple mapping of the observed luminosities (radio cores, X-ray and/or bolometric) into the physical quantities related to any growing black hole: its accretion rate and the released kinetic power. With the aid of such a classification, we have solved the continuity equation for the black hole mass function using the locally determined one as a boundary condition, and the hard Xray luminosity function as tracer of AGN growth rate distribution, supplemented with luminosity-dependent bolometric corrections (Marconi et al. 2004) and absorbing column density distributions. By construction, our adopted intrinsic hard X-ray luminosity functions satisfy the XRB constrain (as well as the sources number counts, and many others), as they adopt the most recent synthesis model as a input (see Gilli et al. 2007 for details). One remaining uncertainty with these models lies in the actual luminosity and redshift distribution of heavily obscured (Compton thick) AGN. However, according to the adopted column density distribution, the contribution to the overall mass growth from Compton thick AGN amounts to just about 20%. Firmer conclusions on this issue will need more accurate determination of number density and typical luminosity of these heavily absorbed (and thus invisible in the 2-10 keV band) sources at redshift between 1 and 2. In fact, preliminary results based on mid-infrared selection criteria (Daddi et al. 2007;Fiore et al. 2008) seem to show quite a good agreement with the predictions of the XRB synthesis models adopted here.
The main results of our study are summarized below.
• We have shown how using standard Soltan (1982) type of arguments, i.e. comparing the local mass density to the integrated mass growth in AGN phases, very tight constraints can be put on the average radiative efficiency of the accretion process: 0.065 ξ 0 (1+ξ lost ) < ∼ ǫ rad < ∼ 0.070 ξ 0 (1−ξ i +ξ lost ) , where ξ0 is the local mass density in units of 4.3 ×10 5 M⊙ Mpc −3 , while ξi and ξ lost are the mass density of z ∼ 5 and "wandering" SMBH, respectively (also in units of the local mass density). Figure 14 shows the constraints derived from eq. (14) in the ξ lost , a * panel, where a * is the mass weighted average spin parameter of the SMBH calculated inverting the classical GR η(a) relation (Shapiro 2005) and taking into account the assumed relationship of eq. (1) between η and ǫ rad . It is interesting to note that both the amount (numbers and masses) of black holes effectively ejected from galactic nuclei due to gravitational wave recoil after a merger and the average radiative efficiency of accreting black holes depend on the spin distribution of evolving SMBH (see e.g. Volonteri 2007, Berti and. Thus, eq. (14) couples implicitly the spin distribution of accreting black holes (through ξ lost and ǫ rad ) and the properties of the seed black hole population, whose density must be reflected in the z ∼ 5 mass density ξi. At face value, our results indicate that, for ξ0 = 1 ≫ ξi, SMBH must have on average a rather low spin, similar to what predicted by models of AGN fuelling that take into account the self-gravity of an accretion disc and result in accretion events which are small and randomly-oriented, thus favoring effective spin-down of the hole (King, Pringle & Hofmann 2008). Figure 14. Joined constraints on the fractional mass density of "wandering" SMBH, ξ lost , ans on the mass weighted spin of SMBH, a * , calculated for ξ 0 = 1. The dark blue area between solid lines is the allowed area for ξ i = 0, while the areas between the lower solid line and the upper dashed (dot-dashed) line represent the allowed area for ξ i = 0.1 (ξ i = 0.3).
• We confirm previous results and clearly demonstrate that, at least for z < ∼ 1.5, SMBH mass function evolves anti-hierarchically, i.e. the most massive holes grew earlier and faster than less massive ones ( § 5.3). By looking at the distribution of SMBH growth times as a function of mass and redshift, for the first time we find hints of a reversal of such a downsizing behaviour at redshifts close to and above the peak of the black hole accretion rate density (z ≈ 2), where we witness the epoch of rapid growth of almost the entire SMBH population. Our results bring the study of AGN 'downsizing' from a phenomenological level, in which one simply describes the observed behaviour of the luminosity functions (Hasinger, Miyaji & Schmidt 2005;Bongiorno et al. 2007), to a more physical one, in which such an observed behaviour is explained in terms of evolution of physical quantities (mass and accretion rate), which could be more directly related to the parallel evolution of the host galaxies population.
• Differently from most previous analytic, semi-analytic and numerical models of black holes growth, we do not assume any specific distribution of Eddington ratios, rather we determine it empirically by coupling the mass and luminosity functions and the set of fundamental relations between observables in the three accretion modes. SMBH always show a very broad accretion rate distribution, and we have highlighted the profound consequences of this fact for our understanding of observed AGN fractions in galaxies, as well as on the empirical determination of SMBH mass functions from large surveys ( § 5.4).
• We have presented the cosmological evolution of the kinetic luminosity function of AGN, based on the evolution of their flat spectrum radio luminosity function and on the empirical correlation between intrinsic radio core luminosity and kinetic power found in Merloni and Heinz (2007). As opposed to the mass growth evolution, the kinetic luminosity function so derived is not very tightly constrained due to poor observational information on the true (in-trinsic) radio core luminosity function at high redshift and to the uncertain distribution of the AGN jet bulk Lorentz factors needed to correct for relativistic beaming of their core emission. Nevertheless, we were able to measure the overall efficiency of SMBH in converting accreted rest mass energy into kinetic power, the "kinetic efficiency" ǫ kin , which ranges between 3 and 5 ×10 −3 , depending on the choice of the radio core luminosity function. This value is somewhat smaller (but less than a factor of 2) than those found independently in both Körding, Jester & Fender (2008) and in Shankar et al. (2008a), where only steep-spectrum sources were used, and the total kinetic power was derived from the correlation between extended radio luminosity and kinetic power of Willott et al. (1999). Nevertheless, the redshift distribution of the kinetic power density released by growing SMBH is in good agreement with that derived by Körding et al. (2008), lending further support to the idea that current measures of the kinetic AGN power are reaching quite robust conclusions.
• Our capability to resolve the mass and accretion rate functions allows us to separate the evolution of both growth rate and kinetic energy density into different mass bins and into the various modes of accretion. In doing so, we found that most (73-80%) of the local mass density was accumulated during radiatively efficient modes (either HR or HK), with the remaining 20-27% in the radiatively inefficient LK mode. The kinetic power density at low redshift is completely dominated by low luminosity AGN, while the contribution from radio loud QSOs (mode HK) becomes significant at z ∼ 2 (see section 6.2). The measured ǫ kin varies strongly with SMBH mass and redshift, being maximal for very massive holes at late times, a property in agreement with what required for the AGN feedback by many recent galaxy formation models.
The richness of details we have been able to unveil for the cosmological evolution of supermassive black holes demonstrates that times are ripe for comprehensive unified approaches to the multiwavelength AGN phenomenology. At the same time, our results should serve as a stimulus for semi-analytic and numerical modelers of structure formation in the Universe to consider more detailed physical models for the evolution of the black hole population.