The Cosmic History of Black Hole Growth from Deep Multiwavelength Surveys

Significant progress has been made in the last few years on understanding how supermassive black holes form and grow. In this paper, we begin by reviewing the spectral signatures of Active Galactic Nuclei (AGN) ranging from radio to hard X-ray wavelengths. We then describe the most commonly used methods to find these sources, including optical/UV, radio, infrared and X-ray emission and optical emission lines. We then describe the main observational properties of the obscured and unobscured AGN population. Finally, we summarize the cosmic history of black hole accretion, i.e., when in the history of the Universe supermassive black holes were getting most of their mass. We finish with a summary of open questions and a description of planned and future observatories that are going to help answer them.


Introduction
Astrophysical black holes come in a wide range of masses, from ∼ > 3 M for stellar mass black holes [1] to ∼10 10 M for so-called supermassive black holes [2,3]. The best evidence for the existence of a supermassive black hole can be found in the center of the Milky Way galaxy, where from dynamical studies, the mass of the Sgr A * source was established to be ∼4.4×10 6 M [4,5].
Evidence for the existence of supermassive black holes has also been found in other massive nearby galaxies [6], mostly from resolved stellar and gas kinematics. For active galaxies, it has been possible to use the technique known as reverberation mapping [7,8,9]. From these observations, a clear correlation has been established between the mass of the central black hole and properties of the host galaxy such as stellar mass in the spheroidal component [10], luminosity [11], velocity dispersion [12,13] and mass of the dark matter halo [14]. The fact that such correlations exist, even though these components have very different spatial scales, suggests a fundamental relationship between black hole formation and galaxy evolution. Furthermore, it is now well established by simulations [15] that the energy output from the growing central black hole can play a significant role in the star formation history of the host galaxy. In particular, theory suggests that nuclear activity regulates star formation either by removing all the gas [16,17] or by heating it [18]. It is therefore obvious, that a complete study of galaxy evolution requires a comprehensive understanding of black hole growth.
Most current black hole formation models tell us that the first black hole seeds formed at z ∼ > 15. While the exact mechanism for the formation of the first black holes is not currently known, there are several prevailing theories (see the comprehensive reviews by M. Rees [19] and M. Volonteri [20] for more details). One of the most popular possibilities is that the first black hole seeds are the remnants of the first generation of stars, the so-called population III stars, formed out of primordial ultra-low metallicity gas. These black holes formed at z∼20 and have typical masses ∼100-1,000 M . This scenario has problems explaining the very high masses, of ∼10 9 M , estimated for supermassive black holes in z∼6 optically-selected quasars [21]. Alternatively, the first black holes could have formed directly as the result of gas-dynamical processes. It is possible for metal-free gas clouds with T vir ∼ > 10 4 K and suppressed H 2 formation to collapse very efficiently [22], possibly forming massive black hole seeds with M∼10 4 -10 5 M as early as z∼10-15. If instead the UV background is not enough to suppress the formation of H 2 , the gas will fragment and form "normal" stars in a very compact star cluster. In that case, star collisions can lead to the formation of a very massive star, that will then collapse and form a massive black hole seed with mass ∼10 2 -10 4 M [23].
Given the current masses of 10 6−9 M , most black hole growth happens in the Active Galactic Nuclei (AGN) phase [2,24]. With typical bolometric luminosities ∼10 45−48 erg s −1 , AGN are amongst the most luminous emitters in the Universe, particularly at high energies and radio wavelengths. These luminosities are a significant fraction of the Eddington luminosity -the maximum luminosity for spherical accretion beyond which radiation pressure prevents further growth -for a 10 8−9 M central black hole. A significant fraction of the total black hole growth, ∼60% [25], happens in the most luminous AGN, quasars, which are likely triggered by the major merger of two massive galaxies [26]. In an AGN phase, which lasts ∼10 8 years, the central supermassive black hole can gain up to ∼10 7 -10 8 M , so even the most massive galaxies will have only a few of these events over their lifetime. Further black hole growth, mostly in low-luminosity (low Eddington rate) AGN, is likely due to stochastic accretion of cold gas, mostly in spiral galaxies [27].
According to the AGN unification paradigm [28,29], a large fraction of these sources, ∼75% locally, are heavily obscured by optically and geometrically thick axisymmetric material, which explains many of the observed differences among different types of active galaxies. In addition, luminosity [30] and cosmic epoch [31] play a significant role. One constraint on the fraction of obscured AGN and its evolution comes from the spectral shape of the extragalactic X-ray "background" (XRB). Thanks to deep X-ray observations at E ∼ < 10 keV performed by Chandra and XMM-Newton, a very large fraction of the X-ray background, ∼80%, has been resolved into point sources [32], the vast majority of them AGN [33]. Several studies, the first of them ∼20 years ago [34], have used a combination of obscured and unobscured AGN to explain the spectral shape and normalization of the X-ray background with overall good results. The latest AGN population synthesis models [35,36] assume an average ratio of obscured to unobscured AGN of ∼3:1 locally, increasing towards lower luminosities and higher redshifts, as well as a fraction of Compton-thick sources (CT; N H >10 24 cm −2 ) of ∼5-10%, consistent with the value observed at higher energies, E=10-100 keV, of ∼5% by INTEGRAL in the local Universe [35,37], lower by factors of ∼3 than expectations of previous population synthesis models [38,39].
In this paper, we review multiwavelength methods used to trace the growth of SMBHs ( §2), the known properties of unobscured and obscured AGN ( §3 and §4 respectively), the cosmic history of black hole accretion ( §5) and prospects for future observations ( §6). Throughout this paper, we assume a ΛCDM cosmology with h 0 =0.7, Ω m =0. 27 and Ω Λ =0.73, in agreement with the most recent cosmological observations [40].

How to trace SMBH Growth?
One of the main features of the AGN emission is that it covers a very wide range of wavelengths, from radio to Gamma-rays (Fig. 1). While unobscured sources are easily detectable and identified by their UV and soft X-ray continuum and their broad optical emission lines, obscured AGN can only be found at longer, mid-IR, wavelengths or in hard X-rays. Of course, selections at different wavelengths have different biases. For example, while radio surveys are not particularly affected by obscuration, they are most likely to detect radioloud sources, which are only ∼10% of the total AGN population at bright fluxes [41]. However, combining different multiwavelength techniques gets us closer to a complete picture. Below, we briefly describe the most popular AGN selection methods, their advantages and drawbacks.

Optical/UV Continuum
Rest-frame optical/UV selection of AGN, in particular for high-luminosity unobscured quasars, is particularly efficient because the spectral shapes of stars and quasars at these wavelengths produce very different broadband colors due to the presence of the "big blue bump" [43] in quasar spectra from ∼100Å to ∼1µm (Fig. 2). This emission is often attributed to the thermal radiation with temperatures ∼30,000K originating in the accretion disk [44]. This unique spectral shape has been used in the past to identify quasars with great success by optical surveys such as the Palomar-Green (PG) survey [45], the 2 degree field QSO redshift survey [46] or more recently the Sloan Digital Sky Survey (SDSS; [47]), which has found now more than 1 million quasars [48]. However, samples selected in the optical are far from complete, as emission at these wavelengths is strongly affected by reddening or extinction, and most AGN are obscured along our line-of-sight. Furthermore, for lower luminosity sources, the optical light from the host galaxy can easily outshine the nuclear emission. This is particularly important for ground-based observations and high-redshift sources, for which it is very hard to separate the non-thermal and stellar components spatially.

Radio
Historically, identification of AGN based on their radio emission has been very important. In fact, the first discovered quasar, 3C 273, was originally classified as a radio source [50]. In spite of this, radio selection can be very problematic. In radio-loud sources (defined as f 5GHz /f B >10 [51]) radio emission is associated with a strong, non-thermal, component, probably originating in a beamed collimated relativistic jet [52]. In radioquiet sources, which are typically ∼3 orders of magnitude fainter at these wavelengths [53], the radio emission corresponds to the long-wavelength tail of the far-infrared dust emission. As a consequence, radio-selected samples are necessarily biased towards radio-loud sources, which represent only ∼10% of the overall AGN population. Figure 2: Composite rest-frame optical/UV spectrum for the optically-selected quasars in the Sloan Digital Sky Survey, from the work of Vanden Berk et al. [49]. The dashed and dotted lines show power-law fits to the continuum emission.

Optical Emission Lines
As first reported by Baldwin et al. [54], the photoionizing spectrum of a power-law continuum source, such as an AGN, produces very different emission line intensity ratios when compared with that of typical star-forming regions (mostly due to O and B stars). Hence, emission lines can be used to identify the presence of AGN even in galaxies in which the optical/continuum does not show any direct AGN signature, due to obscuration and/or low luminosity. Because the AGN ionizing emission reaches material even a few kilo-parsecs away from the nuclear region, this selection technique is less sensitive to circumnuclear obscuration, and thus provides a more complete AGN view when compared with, for example, optical/UV continuum selection. This technique was used successfully in the SDSS [55,56] to extend the low-redshift AGN sample to lower luminosities. Emission line ratios and diagnostic regions can be seen in Fig. 3. Emission-line selection can also be used at higher redshifts, as shown by the DEEP2 galaxy redshift survey, which selected a sample of 247 AGN at z∼1 from optical spectroscopy using the DEIMOS spectrograph at the Keck observatory [57].
While this is an efficient AGN selection technique, optical spectroscopy is very expensive in telescope time, and is only feasible for relatively bright emission line regions. This selection may be incomplete at the low luminosity end, if the host galaxy can outshine the high-ionization emission lines. It is currently very difficult to extend this selection beyond z∼1, as the relevant emission lines move to observed-frame near-IR wavelengths, where current-generation spectrographs are significantly affected by atmospheric emission, do not cover wide field of views and have limited multi-object capabilities.

X-rays
As was found more than 30 years ago, AGN are ubiquitous X-ray emitters [59]. Their X-ray emission extends from ∼0.1 keV to ∼300 keV and is attributed to inverse-Compton scattering due to high-energy electrons in a hot corona, surrounding the accretion disk. The high-energy cutoff at ∼100-300 keV is presumably due to a cutoff in the energy distribution of the electrons in the hot corona. AGN are typically ∼1-5 orders of magnitude  [56]. Divisions between star-forming galaxies and AGN are shown by the dashed [55] and solid [58] lines; bona-fide AGN sit at the upper right in these distributions.
more luminous in X-rays than normal galaxies, which makes them the dominant extragalactic population at these wavelengths. Most AGN are obscured by photoelectric absorption by gas along the line of sight, which preferentially affects emission at the lower energies. This is often parametrized by the amount of neutral hydrogen column density in the line of sight. Fig. 4 shows typical AGN X-ray spectra including the power-law component and high energy cutoff, for different levels of photoelectric absorption. Deep X-ray surveys with Chandra and XMM-Newton have found the largest AGN densities, ∼7,000 deg −2 [60], ∼10-20 times higher than even the deepest optical surveys. Still, X-ray selected AGN samples are still biased against the most obscured sources. In fact, even the deepest Chandra surveys can miss more than half of the total AGN due to a combination of obscuration and/or low luminosity [61]. Observing at higher energies helps to alleviate the effects of obscuration. Wide-area surveys performed by the Swift [62] and INTEGRAL [63] satellites have detected a large number of AGN in the local Universe at E=10-100 keV, where all but the most obscured, Compton-thick, sources emit strongly. Unfortunately, due to their relatively high flux limits (∼2-3 order of magnitudes shallower than the deeper Chandra observations), surveys at these energies are so far limited to low redshifts only.

Infrared
Intrinsic AGN emission is not particularly strong at near/mid-IR wavelengths. Radiation coming from the accretion disk, often characterized as a power-law, while very strong at UV and optical wavelengths decreases rapidly beyond ∼1 µm [64]. However, as was originally found by IRAS [53] and latter confirmed by ISO [65] and Spitzer [66], AGN are luminous IR sources. This is commonly attributed to re-emission of absorbed energy by dust. This radiation can be found starting at ∼2-3 µm, which corresponds to the dust sublimation temperature, about ∼1,000-2,000 K [67], and extends to ∼100 µm, at the tail of the black body spectrum for a ∼100-1,000 K distribution [68], where the emission from the host galaxy typically starts to dominate. Typical AGN infrared luminosities are 10 44−46 erg s −1 and thus they represent a significant fraction, ∼30% on average [69], of the bolometric luminosity.
One clear advantage of infrared AGN selection is that this re-emission is mostly isotropic, and hence both obscured and unobscured sources have similar detection probabilities. Therefore, it provides a complementary approach to the most common selection techniques described above, which are less sensitive to obscured sources. However, star-forming galaxies are very luminous at these wavelengths as well, so that the host galaxy can easily outshine the central emission, in particular for low-luminosity sources, hence yielding a low efficiency in detecting AGN [70]. Furthermore, the selection function in infrared studies is more complicated, as the probability of detecting the AGN depends on the properties of the host galaxy, such as the amount of dusty star-formation, which is in principle independent of the nuclear luminosity.

Unobscured AGN
Because quasars are the most luminous, and thus easily detectable, members of the AGN family, the luminosity function of optical quasars has been well determined for years. In particular, it was found that the number of quasars evolve strongly [71] and peak at z 2 [72]. This evolution has been modeled either as a pure luminosity evolution (PLE), in which the characteristic luminosity changes with redshift while the shape of the luminosity function remains the same [73], or a pure density evolution (PDE), so that only the normalization of the luminosity function depends on redshift [71]; however, it was quickly discovered that at least in the PG quasar survey, the shape of luminosity function also evolves with redshift and thus neither a PLE nor a PDE provides a good description [45].
In Fig. 5 we show one of the latest measurements of the quasar luminosity function at 0.4<z<2.6, from Croom et al. [74]. They conclude that a luminosity-dependent density evolution provides a better fit to the optical quasar luminosity function, than either a PLE or PDE. Similar conclusions were reached by studying X-ray selected sources using soft X-ray (0.5-2 keV) observations [75] or hard X-ray (2-10 keV) data [76,77] which includes both obscured and unobscured AGN. However, taking advantage of the large number of sources in their sample, ∼10,000, Croom et al. [74] concluded that the best fit to the observed luminosity function is obtained by using a model based on a luminosity evolution + density evolution (LEDE). The most important difference between a LEDE and a PLE fit is a change in amplitude and bright-end slope at high redshift.
An important conclusion obtained from observations of the quasar luminosity function is the evidence for "cosmic downsizing" [81], i.e., that the most massive black holes get most of their mass at high redshift, while at low redshift only low mass black holes are still growing. This is observed both in optical [74] and hard X-ray luminosity functions [76,81], thus indicating that this result is independent of obscuration. Recent deep optical surveys such as the Great Observatories Origins Deep Survey (GOODS; [82]), the Cosmic Evolution Survey (COSMOS; [83]), the NOAO Deep Wide Field Survey (NDWFS) and the Deep Lens Survey (DLS;  [78] blue), the 2dF-SDSS LRG and QSO survey (2SLAQ; [79]; green) and combining the latter with the general SDSS QSO sample ( [80]; red).The dotted line shows the measured luminosity function in the 1.53<z<1.81 range for reference. The strong evolution of the QSO luminosity function is clearly seen in this figure. This evolution is best fitted by a LEDE model, as described in the text.
[84]) have produced significant advances in extending the quasar luminosity function to higher redshifts, z>3. Fig. 6 shows the quasar luminosity function at z∼4 and the redshift dependence of the quasar spatial density [82,84]. While the presence of downsizing is clear up to z∼2.5, at higher redshifts it is less convincing, most likely due to poor statistics and incompleteness. As argued by Glikman et al. [85], the slope of the faint end of the quasar luminosity function is critical in determining the contribution of these sources to the ionization of the intergalactic medium. Based on current measurements, quasars contribute ∼60% of the ionizing photons at z∼4 and thus are the dominant source at this redshift.
At even higher redshifts, z∼5-6, current deep surveys do not cover enough area to detect a significant number of sources. However, wide-area survey such as SDSS [87] and the Canada-France High-z Quasar Survey [88] have been able to find a sizable sample, ∼40 high-luminosity quasars, at these high redshifts. According to these samples, there is a large decrease in the number density of high-redshift quasars, when compared to z∼2, suggesting that the peak of the quasar activity is at z∼2.5 [80]. The decline towards higher redshifts is given by 10 −0.47z [87] from z=3 to z=6. A similar trend is also observed for high luminosity sources in X-ray surveys [89], which are also dominated by unobscured sources. This indicates that due to their very low spatial density, unobscured quasars do not contribute significantly to the early hydrogen re-ionization of the intergalactic medium at z∼6 [87,88], in contrast to the situation at z∼4.

Obscured Accretion
The space density and evolution of the unobscured AGN population has been well studied, mostly from optical and soft X-ray surveys. However, we know that a large fraction of the SMBH growth happens in heavily obscured systems. Observations of the nearest AGN suggest that the local ratio of obscured to unobscured sources is ∼4:1 [90]. A similarly high fraction of obscured AGN has been used to explain the spectrum and normalization of the extragalactic XRB, as shown by the latest AGN population synthesis models [35, 38, 39, [85]. The filled red circles show the measurements presented on that work, while the open red circles were reported previously by the same group [84] and show the change with increased spectroscopic completeness. The triangles were obtained from the SDSS quasar sample at z=4.25 [80]. The blue squares are the space densities of z∼4 QSOs from Ikeda et al. [83] and the dot-dashed line is their best-fit double power law. The lower right-hand legend lists the best-fit parameters to a double power-law (solid line) shaded region represents the 1σ uncertainties. Dashed and dotted lines show the z∼3 QLF [86], representing different fits to the observed QLF. Right panel: Quasar space density as a function of redshift from the work of Ikeda et al. [83]. Dotted lines used the combined 2SLAQ, SDSS, SWIRE, NDFWS and DLS samples, while the dashed lines combine the COMOS and 2SLAQ sources. While AGN downsizing is clearly visible at z<2.5, at higher redshifts the situation is more uncertain. Figure 7: Observed spectrum of the extragalactic XRB from HEAO-1 [92], Chandra [32], XMM [93], INTE-GRAL [94] and Swift [95] data. The thick black solid line shows the population synthesis model for the XRB spectrum of Treister et al. [35]. Red, blue and thin black solid lines show the contribution to this model from unobscured, obscured Compton thin and CT AGN respectively. 91]. The XRB gives an integral constraint to the AGN population and its evolution; the most recent deep surveys show that ∼90% of the observed 2-8 keV XRB radiation can be attributed to resolved AGN [32]. In Fig. 7, we show the latest AGN population synthesis models for the XRB which uses a local ratio of obscured to unosbcured AGN of ∼3:1, plus a luminosity and redshift dependence, as described below [35]. The largest uncertainty in such models stems from the normalization mismatch of the data in the 1-10 keV energy range.
A possible dependence of the fraction of obscured AGN on luminosity was first suggested nearly 20 years ago [30], and confirmed since by hard X-ray surveys [76,81,96]. A possible physical explanation, is the socalled "receding torus," in which the size of the inner opening angle depends on luminosity [30,97]. More recent observations found a luminosity dependence of the ratio of mid-IR to bolometric flux for unobscured AGN, consistent with this idea [69]. Alternatively, it has been proposed that the observed dependence of the obscured fraction on luminosity could be explained either by the effects of photo-ionization on the X-ray obscuring matter [98] or by the Eddington limit on a clumpy torus [99]. In Fig. 8 we show the observed fraction of obscured AGN as a function of luminosity obtained by combining data from deep Chandra X-ray surveys [100]. A consistent result is observed for AGN selected in a hard X-ray (E=14-195 keV) survey, as shown in the right panel of Fig. 8, indicating that this trend is not due to selection biases.
In the left panel of Fig. 8 we compare the observed dependence of the fraction of obscured sources on luminosity with the expectations for different geometrical parameters of the obscuring material. If the height of the torus is roughly independent of luminosity, the change in covering fraction is due to a change in inner radius (the original "receding torus" model), hence a rough L −1/2 dependence for the contrast should be expected [30,101]. If the effects of radiation pressure are incorporated, in the case of a clumpy torus, a L −1/4 dependence is expected. As can be seen in Fig. 8, a L −1/2 dependence is too steep compared to observed data. This implies that the height of the obscuring material cannot be independent of the source luminosity and provides evidence for a radiation-limited structure.
The dependence of the fraction of obscured AGN on redshift is more controversial. While some studies [31,91,102,103,104] found a small increase in the fraction of obscured AGN at higher redshifts, other results suggest that this fraction is constant [76,105]. These discrepancies can be understood due to a combination of small samples and the use of N H to classify AGN, which produces a well-know redshift bias [105]. The fraction  [31]. Black circles with solid error bars show the fraction obtained combining these two samples. The dependence in the AGN population synthesis model of Gilli et al. [39] for the intrinsic and observed fractions of obscured AGN are shown by the dashed and dotted red lines; the dependence used by Treister & Urry [31] is shown by the solid blue line. The dashed magenta line shows the expected dependence for a radiation-limited torus [99], while the dotted green line shows the expectation for the original "receding torus" [30], both normalized to the observed value in the 10 42−43 erg s −1 bin. Right panel: Same as left panel but using an AGN sample at z∼0 selected in hard X-rays from Swift/BAT observations [37]. The fact that the same luminosity dependence is observed in both samples indicates that it is not due to selection effects.
of obscured AGN as a function of redshift for a large, ∼2,000 sources, X-ray selected sample [100], using optical emission lines to separate obscured and unobscured AGN, is shown in Fig. 9. It increases significantly with redshift, roughly as (1+z) α , with α=0.3-0.5 (thin dashed lines, bottom panel, Fig. 9; best fit, α 0.4, thick dashed line). This value of α does not change significantly if a different host galaxy evolution is assumed, and it is consistent with the value of 0.3 reported by other studies [91,102,104].
Since star-forming galaxies may be expected to have more dust, the increase in the relative fraction of obscured AGN at high redshift may be due to an increase in the contribution to obscuration by galactic dust. By combining hard X-ray and mid-infrared observations, a similar ratio of hard X-ray to mid-infrared flux for obscured and unobscured AGN has been found [106], contrary to the predictions of the simplest AGN unification paradigm, in which the obscuration comes from the dusty torus and therefore the mid-infrared emission is reduced due to self-absorption. This result can be explained if the obscuration comes from a much more extended region, i.e., kiloparsec, galactic scales rather than a compact parsec-scale torus. Furthermore, signatures for extended absorbing regions have been detected in nearby galaxies like NGC 1068 [107] and NGC 4151 [108]. Heavy absorption at kiloparsec-scales has routinely been found in ultra-luminous infrared galaxies (ULIRGs), which suffer a very strong evolution [109]. Hence, it seems likely that the change in the relative fraction of obscured AGN could be related to galactic-scale absorption, in particular since some ULIRGs also contain an obscured AGN (e.g., Arp 220; [110]).
Below, we review in more detail our current knowledge of the obscured AGN population at three different cosmic epochs: z 0, z=1-3 and z>6.

Obscured AGN in the Local Universe
Nearby AGN are found in the so-called Seyfert galaxies [111], which are know to host low luminosity and/or obscured nuclear activity [90,112]. These growing supermassive black holes have been identified because of ; gray circles) and combining both samples (black circles with solid error bars). The expected observed fraction for an intrinsic fraction of 3:1 obscured to unobscured AGN, accounting for optical and X-ray selection effects, is shown by the black solid line. As can be seen, while the observed fraction of obscured AGN declines toward higher redshifts, if the X-ray and optical selection effects and the luminosity dependence of the obscured AGN fraction are taken into account, this decline should be even stronger. Bottom panel: Inferred fraction of obscured AGN relative to an intrinsically constant fraction after correcting for selection effect and including the luminosity dependence of the obscured AGN fraction. Symbols are the same as for the upper panel. The corrected fraction of obscured AGN increases with redshift following a power-law of the form (1+z) α with α=0.4±0.1.
their high-ionization optical emission lines and in some cases their blue UV/optical continuum. The first nearby AGN catalogues, produced ∼40 years ago [113], contained ∼200 quasars. Roughly speaking, ∼5-15% of the galaxies near the Milky Way contain an active nucleus [90], and ∼75% of these active galaxies are obscured. In fact, two of the three nearest AGN are Compton-thick (NGC 4945 and the Circinus galaxy [114]). Hence, optical surveys are not particularly efficient in unveiling this accretion, while observations at other wavelengths, in particular in the infrared [115], and hard X-rays, are more complete.
Surveys at hard X-ray energies, E>15 keV, have been very successful in providing the most complete AGN samples in the local Universe. As long at the neutral hydrogen column density is lower than ∼10 24 cm −2 , the direct AGN emission is mostly unaffected at these energies. Current observations at E>10 keV with the International Gamma-Ray Astrophysics Laboratory (INTEGRAL; [116]) and Swift [117] satellites are available only at relatively high fluxes, and hence low redshifts, z<0.05.
Using the IBIS coded-mask telescope [118], INTEGRAL surveyed ∼80% of the sky down to a flux of 5 mCrab in the 17-60 keV. Krivonos et al. [119] report the properties of 130 AGN detected in these all-sky observations. A large number of unidentified sources remain in this full INTEGRAL catalog (48) but only seven are found at high galactic latitude (|b|>5 o ), and thus are likely of extragalactic origin. Five of the 130 known AGN are Compton thick. Using similar observations from the all-sky Swift/BAT survey, a catalog of 103 AGN [62] contains five AGN with estimated N H greater than 10 24 cm −2 . However, we caution that some of these N H measurements were obtained by fitting a single absorbed power-law to the X-ray spectrum, while heavily absorbed AGN have more complex spectra [120,121], so the N H estimates are likely to be lower limits. Figure 10 shows the cumulative number counts of AGN, with CT sources shown separately, as a function of hard X-ray flux. In order to avoid the necessity of specifying a standard spectrum to convert fluxes to different energy bands, we show the INTEGRAL and Swift sources separately, but note that there is good agreement (within ∼40%) in the normalization between the two distributions if a standard band conversion is assumed. At these high fluxes the slope of the log N-log S is Euclidean, implying a uniform spatial distribution, as expected given the low redshifts of these sources. The number of CT AGN found by these surveys is surprisingly low, compared to the sample of known CT AGN in the local Universe, most likely due to the effects of obscuration even at these high energies [35,37,122]. A study of optically-selected local Seyfert 2 galaxies with hard X-ray information [112] found 12 CT AGN in a total of 45 Seyfert galaxies. Three were detected by Chandra and/or XMM, while the rest are mostly reflection-dominated sources, too faint to be detected by either INTEGRAL or Swift even though they are nearby, moderate luminosity AGN. This suggests that even hard X-ray surveys miss quite a few Compton thick AGN.
The observed fraction of CT AGN in the INTEGRAL and Swift/BAT hard-X-ray selected samples is low, ∼5%. A very similar and consistent value, 4.6%, was recently obtained from a sample of 307 objects detected in the three-years all-sky Swift/BAT survey [37]. This was initially surprising, since previous AGN population synthesis models that can explain the XRB used much higher CT fractions of ∼15-20% [38,39], i.e., factors of 3-4 higher. We now know that even observations at these high energies can be affected by obscuration, if the column density is high enough. For example, ∼50% of the source flux in the 15-55 keV range can be lost if logN H >24.5 [123]. As pointed out by Malizia et al. [122], and as can be clearly seen in Fig. 11, the INTEGRAL all-sky observations, which have a similar flux limit as the Swift/BAT images, show a steep decline in the number of obscured sources, from ∼80% at z<0.015, to ∼20-30% at higher redshifts. Also, all the CT AGN in this sample were found at z<0.015. Hence, these authors concluded that the INTEGRAL observations are affected by obscuration at larger distances, and that the intrinsic fraction of CT sources is ∼25%, as observed at z<0.015. However, it is worth mentioning that these additional sources, because of their very high column densities, do not contribute significantly to the XRB -although they certainly contribute to black hole growth.
The cumulative contribution of CT AGN to the XRB, as a function of redshift, determined from population synthesis models, is shown in Fig. 12. As can be seen, the total contribution of CT AGN to the XRB is ∼9%, and about 50% of it comes from sources at z<0.7. Similarly, only ∼2% of the XRB is provided by CT AGN at z>1.4, while CT AGN at z>2 contribute only ∼ < 1% to the XRB. Conversely, the 5% uncertainty in the absolute measurement of the XRB intensity translates into an uncertainty of a factor of ∼5 in the number of CT AGN at  [38], which at these high fluxes has a Euclidean slope. The dashed lines mark the Euclidean slope normalized to the number of Swift and INTEGRAL CT AGN (5% of the total). The gray lower limits show the previously-known transmissiondominated AGN with hard X-ray observations, not detected in the INTEGRAL or Swift surveys. These are lower limits since they were selected from pointed observations and are thus highly incomplete.  [122]. All of these sources have X-ray luminosities lower than 10 46 erg s −1 . The remarkably strong decline in this fraction at z∼0.015 is likely due to selection effects and not intrinsic to the AGN population. That the fraction of CT AGN in the first bin is ∼20%, in contrast to the ∼5% found overall suggests current hard X-ray surveys are not sensitive enough to observe CT AGN beyond z∼0.01. This will improve dramatically with NuSTAR [124].   keV Swift/BAT band as a function of redshift from population synthesis models [35]. As shown by the vertical dashed lines, 50%, 80% and 90% of the total CT AGN contribution come from sources at z<0.7, 1.4 and 2.0, respectively. Only ∼1% of the total XRB intensity comes from CT AGN at z>2. The data point at z∼0 corresponds to the contribution to the XRB by the CT AGN detected by Swift/BAT, while the data points at high redshift were obtained from the CT AGN in the CDF-S [125]. Solid error bars correspond to transmission-dominated sources only, while the data points with dashed error bars include all the sources in the sample.
z>2. Hence, the number of CT AGN at high redshift is largely unconstrained by the XRB.

Obscured AGN at Intermediate Redshifts (1 ∼ < z ∼ < 3)
As shown in the previous section, the number of heavily-obscured AGN at intermediate redshifts, z ∼ > 1, is largely unconstrained by the XRB due to model degeneracies, or by current X-ray surveys at E>10 keV, which do not have the required sensitivity. NuSTAR will change this situation dramatically. However, for now we are forced to use alternative methods to determine the amount of black hole growth occurring in these sources. We explore here two of these techniques, which have been particularly successful: X-ray stacking and mid-IR AGN selection.
Using the deepest available X-ray observations obtained with Chandra, as in the example shown in Fig. 13, a number of moderately-obscured AGN have been identified at z>1. The observed-frame hard X-ray band (2-8 keV band) at these redshifts covers higher rest-frame energies, thus making them less affected by obscuration. For example, the vast majority of the sources the GOODS fields with high X-ray to optical flux ratios are obscured at z∼2 [126,127,128]. Furthermore, from X-ray spectral analysis, ∼30 CT AGN candidates have been identified in the CDF-S [125] and CDF-N [129]. However, it is clear that X-ray selection remains highly incomplete for obscured sources at these redshifts [61].
Because much of the energy absorbed at optical to X-ray wavelengths is later re-emitted in the mid to far-IR, it is expected that AGN, in particular the most obscured ones, should be very bright mid-IR sources [66]. Sources having mid-IR excesses, relative to their rest-frame optical and UV emission, have been identified as potential CT AGN candidates at z∼2 [130,131,132,133,134]. However, because of the strong connection between vigorous star formation and AGN activity in the most luminous infrared sources [26], the relative contribution of these two processes remains uncertain and controversial [135,136,137]. Significant progress has been made by virtue of deep Spitzer observations, in particular using the 24 µm band. At z∼1-2, this There are ∼760 individual sources in ∼450 arcmin 2 . This is a smooth enhanced image corresponding to the full (0.5-8 keV) X-ray band. Image and data obtained from http://cxc.harvard.edu/cdo/cdfs.html. emission corresponds to rest-frame wavelengths of ∼10 µm, where the contrast between AGN and star formation is largest. In order to look for high-luminosity obscured AGN missed by X-ray observations, Fiore et al. [131,138] defined the "mid-IR excess" region as f 24 /f R >1000 and R-K>4.5 (Vega), as shown in Fig. 14. As argued by several authors, the fraction of CT AGN in these infrared-excess samples is very high, >70% [131,134,138].  The solid histogram shows the distribution for the sources not detected in X-rays, while the dashed hatched histogram considers only the X-ray detected sources. A KS test shows that it is perfectly likely (∼16%) that these two distributions were drawn from the same parent distribution. The dotted histogram shows the slightly lower redshift distribution (divided by 1,000) for all the sources with a 24 µm detection and a measured photometric redshift in the ECDF-S.
Because sources in the mid-IR excess region are, by definition, very faint at optical wavelengths, it has been very difficult to use optical spectroscopy to measure accurate redshifts. Instead, most surveys have had to rely on (hopefully) accurate photometric redshifts [139,140]. The distribution of photometric redshifts for the sources in the mid-IR excess region in the ECDF-S is shown in Fig. 15; most mid-IR excess sources, have 1<z<3. While the majority of these sources are not detected in X-rays, a significant signal is found in X-ray stacks [131,134,138]. As shown in Fig. 16, the strong stacked detection at E>5 keV clearly indicates the presence of a large number of heavily-obscured AGN in this infrared-excess sub-sample. Specifically, Treister et al. [134] reported that heavily-obscured AGN were ∼80-90% of the mid-IR-excess sources in the ECDF-S. A similarly high fraction of ∼80% was found in the CDF-S [131] and other fields [138]. Optical spectral fitting of these sources indicates evidence for substantial young stellar populations, younger than 100 Myrs [134]. This suggests that these sources are simultaneously experiencing significant star-formation and heavilyobscured AGN activity. The best-fit stellar masses for ECDF-S infrared-excess sources range between 10 9 and 10 12 M with a median stellar mass of ∼10 11 M [134]. Hence, in general these are very massive galaxies.
In order to study in more detail the evolution of the CT AGN space density, in Fig. 17 we present the existing measurements of the CT AGN space density as a function of redshift. Reasonable agreement, in particular at z<1, is found between both observed values and existing hard X-ray luminosity functions and evolution [77]. However, at z>1.7 and high luminosity, a clear discrepancy is found. Treister et al. [35], concluded that this difference of a factor of 2-3 could be due either to incompleteness in the Swift/BAT and INTEGRAL CT AGN samples at z=0 used to fix the luminosity function normalization (because reflection-dominated AGN are missed) or to contamination by other types of sources in the observed values at high redshifts. However, after adding the measurements obtained using the infrared-excess sources in the ECDF-S, it appears not only that the systematic difference is still present but perhaps more importantly that there is a strong increase in the number of CT AGN from z 1.7 to 2.4. This is not described by any existing luminosity function. It is unlikely that this evolution is due to selection effects, as results from different fully independent surveys and selection techniques are combined in Fig. 17, namely X-ray selected sources [125], 24 µm-selected sources Figure 16: Stacked background-subtracted Chandra counts as a function of rest-frame energy, for sources with f 24 /f R >1000 and R-K>4.5 in the 4 Msec CDF-S field (filled circles; [134]). The cyan dashed lines (stars) shows the simulated spectra for the high-mass X-ray binary (HMXB) population normalized using the relation between star-formation rate and X-ray luminosity [141]. The blue dashed lines (open squares) show simulated thermal spectra corresponding to a black body with kT =0.7 keV. An absorbed AGN spectrum, given by a power-law with Γ=1.9 and a fixed N H =10 24 cm −2 , is shown by the red dashed lines (open circles). In addition, a scattered AGN component, characterized by a 1% reflection of the underlying unobscured power-law, is shown by the green dashed lines (open triangles). The resulting summed spectrum (black solid lines) is in very good agreement with the observed counts. The strong detection in the stacked spectrum at E>5 keV, confirms the presence of a significant number of heavily-obscured AGN in these IR-excess objects [134]. [134,138], and a sample of CT AGN found using mid-IR spectroscopy [133]. This result can be interpreted in the context of galaxy evolution models [142], where quasar activity is driven by galaxy mergers and the supermassive black hole is initially completely surrounded by dust, before radiation pressure removes it and a "classical" unobscured quasar is visible [25].
Several groups [143] have found that the fraction of galaxies containing an AGN is a strong function of their IR luminosity. In Fig. 18 we present the stacked spectra for the sources in the CDF-S, grouped in bins of IR luminosity [144]. We can see by comparing these spectra that the relative emission at E 5 keV, where we expect the AGN emission to dominate even for heavily-obscured sources, changes with IR luminosity. In other words, there is a clear trend, with stronger high energy X-ray emission at increasing IR luminosity. The spectra shown in Fig. 18 cannot be directly interpreted, as the detector-plus-telescope response information is lost after the conversion to rest-frame energy and stacking. Hence, simulations assuming different intrinsic X-ray spectra have to be used in order to constrain the nature of the sources dominating the co-added signal.
The observed stacked spectral shape cannot be explained by any plausible starburst spectrum. An AGN component dominating at E>5 keV, is required. The average intrinsic rest-frame 2-10 keV AGN luminosity needed to explain the observed spectrum, assuming that every source in the sample contains an AGN of the same luminosity, is 6×10 42 erg s −1 for sources with L IR >10 11 L , 3×10 42 erg s −1 for sources with L IR >5×10 10 L , 5×10 41 erg s −1 for 5×10 10 L >L IR >10 10 L and 7×10 41 erg s −1 for L IR >10 10 L . All of these are (intrinsically) very low-luminosity AGN; even if there is a range, it is extremely unlikely to include high-luminosity quasars like those discussed in previous stacking papers. This is not too surprising, actually, because the surveyed volume (even to high redshift) is small, so rare objects like high-luminosity Figure 17: Space density of Compton thick AGN as a function of redshift, as published by Treister et al. [134]. Filled triangles show the ECDF-S results [134]. Squares: X-ray selected sources in the CDF-S [125]. Star: Measurement obtained using mid-IR spectroscopy [133]. Pentagons: Values obtained using mid-IR excess sources in COSMOS [138]. Solid lines show the expected space density of Compton thick AGN from the luminosity function of Yencho et al. [77], with the overall normalization fixed to the results of the INTEGRAL and Swift/BAT surveys [35], while the dashed lines show the expectations based on the luminosity function of Della Ceca et al. [103]. Red symbols show measurements and expectations for L X >10 43 erg s −1 sources, while the blue symbols are for L X >10 44 erg s −1 . While for the lower luminosity sources a good agreement is found between observations and expectations, higher luminosity sources at z>1.8 lie well above the luminosity function.
quasars do not appear. If the heavily-obscured AGN in these stacked samples have the same median intrinsic luminosity as the X-ray detected sources with similar IR luminosities, this would indicate that 15% of the galaxies with L IR >10 11 L contain a heavily-obscured AGN. This fraction is ∼10% in the L IR >5×10 10 L and 5×10 10 L >L IR >10 10 L samples. For sources with L IR >10 10 L this fraction is <5%. This extra AGN activity (in addition to the X-ray detected sources) can account for ∼22% of the total black hole accretion. Adding this to the obscured black hole growth in X-ray detected AGN [145], we confirm that most of this growth, ∼70%, is significantly obscured and missed by even the deepest X-ray surveys [25,61].

Obscured AGN at High redshifts, z>3
As mentioned in §3, most measurements of black hole accretion at high redshift, z>3, come from optical observations of unobscured sources. This is not only because obscured sources are obviously fainter at most wavelengths, but also because large areas have to be covered in order to survey a significant volume at high redshifts. As can be seen in Fig. 19, there is a clear decline in the number of luminous quasars at z>2, although the decline is shallower for X-rays compared to optical surveys, since X-ray selection is less biased [77,146]. However, it is important to point out that: (i) these results are limited to the highest luminosity sources, logL X >44.5 erg s −1 , which do not represent the average AGN and do not contribute much to the extragalactic XRB [35], and (ii) only relatively unobscured sources are included. In particular, heavily obscured, Compton-thick, AGN are systematically underrepresented in these surveys. As we will describe, these missing populations can have a significant impact in our understanding of cosmic supermassive black hole growth.
In order to search for the presence of growing supermassive black holes in young galaxies, Treister et  Fig. 16. The resulting summed spectra (black solid lines) are in very good agreement with the observed counts. The strong detection in the stacked spectrum at E>5 keV, in particular at the higher IR luminosities, confirms the presence of a significant number of heavily-obscured AGN in these samples.
al. [152] stacked X-ray images of z>6 galaxy candidates selected based on the optical and near-IR dropout techniques, selected from a sample of 197 galaxies, 151 in the CDF-S and 46 in the CDF-N from Bouwens et al. [153]. Using the 4 Msec Chandra observations of the CDF-S and the 2 Msec data available on the CDF-N, this corresponds to a total exposure time of ∼7×10 8 seconds (∼23 years). Significant detections in both the soft and hard X-ray bands, were obtained. However, these detections have recently been questioned (after the submission of this review) by several authors [154,155,156] due to a possible bias in the background subtraction technique used by Treister et al. While a full discussion is clearly beyond the scope of this paper, we note that a full analysis using an optimal weighting scheme (as in [152]), as well as considering the effects of faint, undetected, sources in the background, remains to be done. At a minimum, the results presented below can be considered as upper limits. The corresponding average rest-frame 2-10 keV luminosity, derived from the observed-frame hard band, is 6.8×10 42 erg s −1 . Since none of these sources were individually detected in X-rays, at least 30% of the galaxies in this sample likely contain an AGN [152]. Furthermore, there is a factor ∼9 difference between the fluxes measured in the observed-frame soft and hard bands. The only explanation for this relatively large flux ratio in the hard to soft bands is very high levels of obscuration. As can be seen in Fig. 20, at z∼6, a minimum column density of N H 10 24 cm −2 , i.e. Compton-thick obscuration, is required. Given that this ratio is observed in a stacked X-ray spectrum, this implies that there are very few sources with significantly lower levels of obscuration, which in turn means that these sources must be nearly Compton-thick along most directions (∼4π obscuration). Similar sources have also been observed in the local Universe [157] but appear to be rare. Furthermore, for z ∼ < 3 we know that the fraction of obscured AGN increases with decreasing luminosity [63,76,102] and increasing redshift [31,91]. Hence, it is not entirely surprising that the sources studied here, given their low luminosities and high redshifts, are heavily obscured. In fact, the discovery of a Compton-thick AGN at z∼5 selected using the drop-out technique has been recently reported [158]. Figure 19: Number density of X-ray selected AGN as a function of redshift for sources with logL 2−10 >44.5 erg s −1 , as published by Brusa et al. [147]. Lines show AGN luminosity functions [39,77,146,148,149] while individual measurements are shown by green circles [150] and black squares [149]. The red star at z=6 was obtained from the optical luminosity function assuming no evolution in α ox . The black shaded area shows the observations of optically-bright SDSS QSOs [151]. Figure 20: Expected ratio of the observed-frame hard to soft flux as a function of obscuring neutral Hydrogen column density (N H ). The black solid line was derived assuming an intrinsic power-law spectrum with slope Γ=1.9 and photoelectric absorption. The gray zone shows the measured ratio for the stack of galaxies at z 6 and the ±1 standard deviation limits. A column density of N H 10 24 cm −2 , i.e. Compton-thick obscuration, is required to explain the observed hard to soft X-ray flux ratio [152].
This relatively large number of growing supermassive black holes at z ∼ > 6 is contrary to the picture obtained from optical observations of high-luminosity quasars at similar distances, which are much rarer [21,151]. This is particularly important in our understanding of the early hydrogen re-ionization, which can be either due to young stars and/or growing supermassive black holes [159]. The results presented by Treister et al. [152] show that while growing supermassive black holes do not contribute much to hydrogen re-ionization, this is not because their numbers drop steeply at z>4 as previously suggested [160], but because large amounts of obscuration found in these sources imply that UV and soft X-rays do not escape.

The Cosmic History of Black Hole Accretion
Direct black hole mass measurements, either through stellar or gas dynamics, are available for only a few nearby galaxies. However, thanks to the tight correlation between mass of the supermassive black hole and other properties such as velocity dispersion and others, it has been possible to estimate the black hole mass function at z 0 [161,162,163,164]. This is commonly done starting from the observed galaxy luminosity or velocity function and assuming either a constant black hole to stellar mass ratio [161] or the M-σ relation [163]. Both the overall shape of the black hole mass function and the integrated black hole mass density, which can only be computed at z 0, can be used to infer properties of the AGN population. This was first used in the so-called "Soltan's argument" [24], which says that the intrinsic bolometric AGN luminosity, L, is directly linked to the amount of mass accreted by the black hole,Ṁ acc : where ε is the accretion efficiency and c is the speed of light. A typical value assumed for the efficiency is ∼10% [24,163].
Recent comparisons of the black hole mass function to the distribution inferred from the observed AGN luminosity indicate that the average efficiency is 8%, the Eddington ratio is ∼50%, and the average lifetime of the visible AGN phase is ∼10 8 years [163,164]. By studying the black hole mass distribution at the high mass end, M>10 9 M , Natarajan & Treister [165] found that the observed number of ultra-massive black holes is significantly lower than the number density inferred from the AGN hard X-ray luminosity function. They concluded that this is evidence for an upper limit to the black hole mass, which can be explained by the presence of a self-regulation mechanism.
The observed black hole mass density at z 0 obtained by integrating the black hole mass function, ranges from 2.9×10 5 [162] [35] obtained a value of 4.5× 10 5 M Mpc −3 , perfectly consistent with the observed value, indicating that at least locally, X-ray detected AGN can account for most or all of the black hole growth.
The black hole mass function can be measured observationally for unobscured, high luminosity AGN at higher redshift, taking advantage of the known correlation between black hole mass and observational quantities such as luminosity and emission line width [167]. These correlations are calibrated using more direct black hole measurements, available for a few, mostly local, sources [8,9]. The large number of unobscured quasars with optical spectroscopy provided by the SDSS and other optical surveys have been very useful in determining the black hole mass function up to high redshifts, [168,169,170], and for lower luminosity sources using deeper surveys such as the AGN and Galaxy Evolution Survey (AGES; [171]) and COSMOS [172]. Using these black hole mass functions as a function of redshift as constraints, recently Natarajan & Volonteri [173] concluded that the observational data are inconsistent with the hypothesis that these black holes are created as the remnants of population III stars. Instead, they argue massive seeding models are required [174].
While a clear picture of the history of black hole growth is emerging, significant uncertainties still remain. In particular, while the spectral shape and intensity of the extragalactic X-ray background have been used to constrain the AGN population, the number of heavily obscured accreting supermassive black holes beyond z∼1 is not properly bounded. Infrared and deep X-ray selection methods have been useful in that sense, but have not provided a final answer, due to confusion with star-forming galaxies in the infrared and the effects of obscuration in X-rays. At higher redshifts, the situation is even more unclear, and only a few, very rare, high luminosity quasars are known. Unless high-redshift AGN luminosity functions are pathological, these extreme sources do not represent the typical growing black holes in the early Universe. As a consequence, and in spite of recent advances [152,173], the formation mechanism for the first black holes in the Universe is still unknown.

Prospects
Scheduled for launch in February 2012, NuSTAR will be the first focusing high energy (E=5-80 keV) X-ray mission, reaching flux limits ∼100 times fainter than INTEGRAL or Swift/BAT observations and comparable to Chandra and XMM-Newton at lower energies. During the first two years of operations, NuSTAR will likely observe, as part of the guaranteed time program, two extragalactic fields: the ECDF-S and the central 1 deg 2 part of COSMOS, for a total of 3.1 Msec each. These deep high-energy observations will enable to obtain a nearly complete AGN survey, including heavily-obscured Compton-thick sources, up to z∼1.5 [175]. A similar mission, ASTRO-H [176], will be launched by Japan in 2014. Both missions will provide angular resolutions ∼ < 1 , which in combination with observations at longer wavelengths will allow for the detection and identification of most growing supermassive black holes at z∼1.
There is little doubt that the Atacama Large Millimeter Array (ALMA) will revolutionize our understanding of galaxy evolution. Sources of mm and sub-mm emission traced by ALMA includes thermal emission of the warm/cold dust, which traces star formation, synchrotron radiation associated with relativistic particles and free-free radiation from HII regions. In particular, CO rotational transition lines have been used to trace the spatial distribution, kinematics, temperature and mass of the molecular gas [177]. The sensitivity of ALMA will allow for the detection of luminous IR galaxies (L IR >10 11 L ), which have been found to often host a heavily-obscured AGN [144], up to z∼10. Furthermore, with ALMA it will be possible to study separately the molecular dust surrounding the central black hole and those in star forming regions in the host galaxy. Due to their limited sensitivity and relatively bad angular resolution, currently-available mm/sub-mm telescopes are not ideal to study star-forming regions even in nearby galaxies. This will dramatically change thanks to ALMA, which will have orders of magnitude better sensitivity and HST-like angular resolution. The first call for ALMA observations was released on March 31, 2011 for observations starting on September 30, 2011. It is expected that the complete array will be in full operation in 2013. The superb spatial resolution and sensitivity of ALMA will allow to uniquely identify the optical/near-IR counterpart of the mm-submm sources. Furthermore, ALMA will directly provide the redshift of the mm-submm sources through the detection of CO rotational transition lines, up to very high redshifts. Combining these new data with existing multiwavelength information will finally allow us to complete the census of supermassive black hole growth since the epoch of cosmic re-ionization.