# Tracing the cosmic growth of supermassive black holes to z similar to 3 with Herschel*

I. Delvecchio, C. Gruppioni, F. Pozzi, S. Berta, G. Zamorani, Alessandro

Cimatti, D. Lutz, D. Scott, C. Vignali, G. Cresci, et al.

► To cite this version:

I. Delvecchio, C. Gruppioni, F. Pozzi, S. Berta, G. Zamorani, et al.. Tracing the cosmic growth of supermassive black holes to z similar to 3 with Herschel*. Monthly Notices of the Royal Astronomical Society, 2014, 439 (3), pp.2736-2754. 10.1093/mnras/stu130. hal-01441011

HAL Id: hal-01441011

https://hal.science/hal-01441011

Submitted on 26 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tracing the cosmic growth of supermassive black holes to $z \sim 3$
with Herschel*

I. Delvecchio,¹† C. Gruppioni,² F. Pozzi,¹ S. Berta,³ G. Zamorani,² A. Cimatti,¹

D. Lutz,³ D. Scott,⁴ C. Vignali,¹ G. Cresci,⁵ A. Feltre,⁶,⁷ A. Cooray,⁸,⁹ M. Vaccari,¹⁰

J. Fritz,¹¹ E. Le Floc'h,¹² B. Magnelli,¹³ P. Popesso,³ S. Oliver,¹⁴ J. Bock,⁹,¹⁵

M. Carollo,¹⁶ T. Contini,¹⁷ O. Le Févre,¹⁸ S. Lilly,¹⁶ V. Mainieri,⁷ A. Renzini¹⁹
and M. Scodeggio²⁰

¹Dipartimento di Fisica e Astronomia, Università di Bologna, via Ranzani 1, I-40127 Bologna, Italy

²INAF – Osservatorio Astronomico di Bologna, via Ranzani 1, I-40127 Bologna, Italy

³Max-Planck-Institut för Extraterrestrische Physik (MPE), Postfach 1312, D-85741 Garching, Germany

⁴Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada

⁵INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

⁶Dipartimento di Fisica e Astronomia, Università di Padova, vicolo Osservatorio 3, I-35122 Padova, Italy

⁷ESO, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany

⁸Department of Physics & Astronomy, University of California, Irvine, CA 92697, USA

⁹Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

¹⁰Astrophysics Group, Department of Physics, University of Western Cape, Bellville 7535, Cape Town, South Africa

¹¹Sterrenkundig Observatorium, Vakgroep Fysica en Sterrenkunde Universeit, Gent, Krijgslaan 281-S9, B-9000 Gent, Belgium

¹²CEA-Saclay, Service d'Astrophysique, F-91191 Gif-sur-Yvette, France

¹³Argelander Institute for Astronomy, Bonn University, Auf dem Hügel 71, D-53121 Bonn, Germany

¹⁴Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK

¹⁵California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

¹⁶Institute of Astronomy, Swiss Federal Institute of Technology (ETH Hönggerberg), CH-8093 Zürich, Switzerland

¹⁷Institut de Recherche en Astrophysique et Planétologie, CNRS, Université de Toulouse, 14 avenue E. Belin, F-31400 Toulouse, France

¹⁸Laboratoire d'Astrophysique de Marseille, CNRS-Université de Provence, rue Frédéric Joliot-Curie 38, F-13388 Marseille Cedex 13, France

¹⁹INAF – Osservatorio Astronomico di Padova, vicolo dell'Osservatorio 5, I-35122 Padova, Italy

²⁰INAF – IASF Milano, via Bassini 15, I-20133 Milano, Italy

## ABSTRACT

We study a sample of Herschel selected galaxies within the Great Observatories Origins Deep Survey-South and the Cosmic Evolution Survey fields in the framework of the Photodetector Array Camera and Spectrometer (PACS) Evolutionary Probe project. Starting from the rich multiwavelength photometric data sets available in both fields, we perform a broad-band spectral energy distribution decomposition to disentangle the possible active galactic nucleus (AGN) contribution from that related to the host galaxy. We find that 37 per cent of the Herschel-selected sample shows signatures of nuclear activity at the 99 per cent confidence level. The probability of revealing AGN activity increases for bright ($L_{1-1000} > 10^{11} L_\odot$) star-forming galaxies at $z > 0.3$, becoming about 80 per cent for the brightest ($L_{1-1000} > 10^{12} L_\odot$) infrared (IR) galaxies at $z \ge 1$. Finally, we reconstruct the AGN bolometric luminosity function and the supermassive black hole growth rate across cosmic time up to $z \sim 3$ from a far-IR perspective. This work shows general agreement with most of the panchromatic estimates from the literature, with the global black hole growth peaking at $z \sim 2$ and reproducing the observed local black hole mass density with consistent values of the radiative efficiency $\epsilon_{\text{rad}}$ ($\sim 0.07$).

**Key words:** galaxies: evolution – galaxies: nuclei – infrared: galaxies.

<footnote>*Herschel* is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>

<footnote>†E-mail: ivan.delvecchio@unibo.it</footnote>
# 1 INTRODUCTION

For more than five decades, since the discovery of quasars, there has been an increasing interest in understanding such extremely bright objects. Though the detailed mechanisms that are the source of the powerful nuclear activity are still debated (Rees 1978; Lodato & Natarajan 2007; Devecchi & Volonteri 2009; Volonteri & Begelman 2010; Ball et al. 2011; Natarajan 2011), it is believed that it is mainly due to mass accretion on to a supermassive black hole (SMBH; Salpeter 1964; Lynden-Bell 1969; Shakura & Sunyaev 1973; Soltan 1982; Rees 1984), allowing them to be revealed as active galactic nuclei (AGN). In the last few years, several studies have highlighted the importance of AGN in the framework of galaxy formation and evolution.

Indeed, many pioneering studies have shown the presence of a correlation between the black hole (BH) mass and other physical properties of its host: galaxy bulge stellar mass, luminosity, velocity dispersion (Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Marconi & Hunt 2003; Gültekin et al. 2009) and dark matter halo mass (Ferrarese 2002).

Observations have suggested a tight AGN/starburst connection through the 'downsizing' scenario, according to which most massive BHs accrete earlier, faster and with more extreme accretion rates than their lower mass counterparts (Ueda et al. 2003; Hasinger, Miyaji & Schmidt 2005; Le Floc'h et al. 2005; Pérez-González et al. 2005; Richards et al. 2006). A similar anti-hierarchical growth seems to be mirrored in the evolution of the cosmic star formation density (SFD; Cowie et al. 1996, 1997; Lilly et al. 1996; Madau et al. 1996; Bell et al. 2005; Juneau et al. 2005; Bundy et al. 2006; Wall, Pope & Scott 2008; Gruppioni et al. 2013), hinting at a deep interplay among the AGN and its host.

Furthermore, semi-analytical models have shown that stellar processes on their own struggle to reproduce the currently observed galaxy population properties, without adding a further source of energy (Benson et al. 2003; Bower et al. 2006; Schawinski et al. 2006; Gabor et al. 2011). This required injection of heat is believed to come from the central active BH, which seems to affect the surrounding galactic environment. Thorough theoretical studies have well established that nuclear accretion processes could impact the star formation history (SFH) of the whole galaxy, making unavailable a fraction of the cold gas required to drive star formation, either by removing or by heating it up (Springel et al. 2005; Croton et al. 2006; Hopkins et al. 2006; Menci et al. 2006). It is clear that a comprehensive study of the AGN accretion history is crucial to shed light on the evolution of galaxies down to the present epoch.

So far, different accretion scenarios have been proposed: major mergers of gas-rich galaxies (Sanders et al. 1988; Di Matteo, Springel & Hernquist 2005; Hopkins et al. 2006; Menci et al. 2008; Hopkins & Quataert 2010); gradually varying processes (Kormendy & Kennicutt 2004), which in turn include internal secular evolution (i.e. large-scale disc instabilities and bars) and external secular evolution (galaxy interactions) and inflow of recycled cold gas from old stellar populations (Ciotti & Ostriker 2007; Ciotti, Ostriker & Proga 2010). Though an accurate determination of the AGN fraction related to each single scenario has not yet been derived, it is currently believed that most luminous AGN ($L_{AGN} > 10^{44}$ erg s$^{-1}$) are triggered by major merger events (Netzer 2009; Santini et al. 2012), while their fainter counterparts are fuelled through secular processes (Lutz et al. 2010; Shao et al. 2010; Rosario et al. 2012).

Recent attempts have been made to trace the SMBH growth over cosmic time (e.g. Marconi et al. 2004; Shankar et al. 2004;

Hasinger et al. 2005; Hopkins et al. 2006; Hopkins, Richards & Hernquist 2007; Merloni & Heinz 2008; Aird et al. 2010; Fanidakis et al. 2012; Mullaney et al. 2012), supporting the connection between black hole accretion rate density (BHAD) and SFD up to $z \simeq 1-1.5$. In particular, Mullaney et al. (2012) suggest that the ratio SFD/BHAD $\sim 1-2 \times 10^3$ in the redshift range $0.5 < z < 2.5$. In agreement with observations, semi-analytical models and hydrodynamical simulations currently hint at an overall decrease of the BHAD from $z \simeq 2$ down to the present day. At higher redshift the situation is more uncertain, since direct diagnostics for BH mass estimates (e.g. reverberation, time variability, etc.) are no longer available, and indirect measurements need to be extrapolated up to these redshifts.

To date the AGN obscured growth is largely uncertain beyond $z \simeq 1-2$, since even the deepest X-ray surveys suffer from a decreasing completeness with increasing redshift and/or increasing obscuration (Brandt & Hasinger 2005). While the dependence of $N_H$ (hydrogen column density) on the intrinsic AGN luminosity is currently well established (Lawrence 1991; Steffen et al. 2003; Ueda et al. 2003; Barger & Cowie 2005; Simpson 2005; Sazonov et al. 2007; Lusso et al. 2013), the evolution with redshift is still debated: on the one hand, a slight evolution with redshift is suggested by some analyses (La Franca et al. 2005; Ballantyne, Everett & Murray 2006; Treister & Urry 2006; Della Ceca et al. 2008: Hasinger 2008; Treister, Urry & Virani 2009; Iwasawa et al. 2012); on the other hand, the space density of the obscured ($\log[N_H/cm^{-2}] > 22$) AGN population seems to be consistent with a no evolution scenario (Ueda et al. 2003; Akylas et al. 2006; Gilli et al. 2010; Vito et al. 2013), at least for low-luminosity ($L_X < 10^{44}$ erg s$^{-1}$) X-ray AGN.

AGN are ubiquitous X-ray emitters (Elvis et al. 1978) and X-ray surveys are only weakly affected by the X-ray host galaxy light. Nevertheless, recent works (Laird et al. 2010; Schawinski et al. 2012; Treister et al. 2012) claim that the host galaxies of high-redshift ($z > 1-2$) AGN could have experienced intense bursts of star formation (e.g. sub-millimetre galaxies and/or ultra-luminous infrared galaxies, ULIRGs), likely hiding heavily obscured AGN with intrinsically high bolometric luminosities, but enshrouded in dense, cold, optically thick gas clouds (see Alexander et al. 2005).

Infrared (IR) observations are only mildly affected by obscuration and can detect warm dust signatures for both unobscured and obscured AGN, giving rise to a complementary perspective to the X-ray one. This is the reason why IR observations represent a crucial tool to investigate and constrain the obscured growth across cosmic time. Recent *Herschel*-based studies are shedding light on the mid-to-far-IR properties of star-forming galaxies (Hatziminaoglou et al. 2010; Hwang et al. 2010; Rowan-Robinson et al. 2010; Elbaz et al. 2011; Sajina et al. 2012; Smith et al. 2012; Magdis et al. 2013; Symeonidis et al. 2013), both in the presence and in the absence of a significant AGN emission component. The importance of IR observations in chasing obscured and unobscured AGN has also been proved through mid-infrared (MIR) photometry, which allows us to develop interesting diagnostics to detect AGN features (Lacy et al. 2004; Stern et al. 2005; Pope et al. 2008; Donley et al. 2012; Juneau et al. 2013; Lacy et al. 2013). However, these techniques are still to be tested on high-redshift ($z > 1$) sizeable samples (Juneau et al. 2011; Trump et al. 2013) and are sensitive to comparatively large AGN-to-host flux ratios in the MIR domain (Barmby et al. 2006; Brusa et al. 2010; Donley et al. 2012). Hence, possible signatures of nuclear accretion for both intrinsically low-luminosity and heavily absorbed AGN might be diluted by the host galaxy light. These complications highlight the need for a complete assessment of the amount of star formation in AGN host galaxies.
The main purpose of this paper is to provide an estimate of the BHAD from an IR survey through a multiwavelength analysis of galaxy spectral energy distributions (SEDs). The Herschel space observatory (Pilbratt et al. 2010) gives an unique opportunity to constrain the far-infrared (FIR) peak of star-forming galaxy SEDs, which can be used as a proxy for their star formation rate (SFR), down to unprecedented flux densities (≈few mJy), reaching to $z > 3$.

The FIR peak is mainly due to optical/UV stellar light reprocessed by dust. However, an additional contribution can arise from warmer nuclear dust heated by material accreting on to the SMBH rather than stellar processes. Spitzer data, together with Herschel capabilities, are paramount to reach a complete assessment of the dust content at different temperatures stored across the IR domain. Herschel does play a key role in investigating the star formation/AGN connection down to intrinsically low-luminosity sources without the need to extrapolate previous stellar+dust emission templates to rest-frame FIR wavelengths. This is the reason why we exploit observations carried out within the PACS Evolutionary Probe (PEP;<fnref n="1" /> Lutz et al. 2011) and the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012) in two different regions: the Great Observatories Origins Deep Survey-South (GOODS-S; Elbaz et al. 2011), the deepest, pencil beam, survey observed by Herschel, and the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007), one of the widest among the PEP fields. Such a combination allows us to sample different luminosity regimes and to span a wide range of redshifts, while preserving a sizeable number of sources.

The paper is structured as follows. The sources in the Herschel-selected sample are described in Section 2, together with a discussion on their counterparts in other wavelengths and their redshift distribution; in Section 3 we describe our approach to decouple the AGN contribution from that related to SF within each global SED and how we evaluate the relative incidence of the nuclear content. In Section 4, we derive the AGN bolometric luminosity function (LF) up to $z \sim 3$, while in Section 5 we present the BHAD evolution over cosmic time, discussing our results and comparing them with previous estimates from the literature. In Section 6, we present our main conclusions.

Throughout this paper, we adopt a Chabrier (2003) initial mass function (IMF) and we assume a flat cosmology with $\Omega_m = 0.27$, $\Omega_\Lambda = 0.73$ and $H_0 = 71 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2 SAMPLE DESCRIPTION

In this work, we exploit large photometric data sets, from the UV to the sub-mm, both in GOODS-S and in COSMOS fields. We take advantage of PACS (Poglitsch et al. 2010) photometry in the FIR (100 and 160 µm; 70 µm data are available for the GOODS-S field only). As a reference sample we use the PEP 160 µm blind catalogue (internal version 1.2 in GOODS-S and version 2.2 in COSMOS). The PEP sample has sources selected above $3\sigma$ detection at 160 µm, corresponding to nominal flux density limits of 2.4 and 10.2 mJy in the GOODS-S and COSMOS fields, respectively. Globally we find 782 and 5105 sources at 160 µm, in the GOODS-S and the COSMOS fields, respectively. The available photometry has been extended up to sub-mm wavelengths (250, 350 and 500 µm) through the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) data of the HerMES (Oliver et al. 2012).

Here we briefly outline the multiband identification procedure, as well as the overall redshift distribution, referring to Berta et al.

(2010, 2011) and Lutz et al. (2011) for further details concerning data reduction, completeness and evaluation of spurious sources and photometric redshift uncertainties.

### 2.1 Multiwavelength identification

PEP-selected sources have been cross-matched with their respective counterparts using a maximum likelihood technique (Sutherland & Saunders 1992; Ciliegi et al. 2001). Starting from each PACS 160 µm detection, we checked for counterparts at shorter wavelengths through a cross-match with detections at 100 µm, 70 µm (in the GOODS-S sample only) and down to 24 µm, still using a maximum likelihood algorithm.

In GOODS-S, we have matched the PACS area in common with the inner region covered by multiwavelength observations in UV, optical and near-IR (NIR). The resulting common area (196 arcmin<fnref n="2" />) is around 65 per cent of the parent sky area covered by PACS. The number of PACS sources within the inner region is 494, which is around 63 per cent of the parent PACS population. This suggests that the surface density of PACS sources remains almost unchanged. In this inner GOODS-S area, the Spitzer-MIPS (Multiband Imaging Photometer for Spitzer) 24 µm selected catalogue (Magnelli et al. 2009) is linked to Spitzer-IRAC (Infrared Array Camera) 3.6 µm positions through prior extraction. MIPS 24 µm detections act as a connecting link between the FIR and optical regimes, since they have been cross-correlated with the Multi-wavelength Southern Infrared Catalog (MUSIC;<fnref n="2" /> Grazian et al. 2006; Santini et al. 2009).

The MUSIC spectral coverage ranges from the far-UV (~3500 Å) up to the MIR (Spitzer-IRAC at 8 µm). As mentioned before, to gather as much information as possible for each source, we consider the Herschel-selected sample within the MUSIC area (196 arcmin<fnref n="2" />). 99 per cent of these sources (494 at 160 µm) have a counterpart in the 24 µm catalogue and in the optical (MUSIC catalogue), above their corresponding flux density limits. Such a high percentage of identified counterparts suggests that 24 µm and optical observations in the GOODS-S field are deeper than the Herschel-PACS ones, with respect to an average MIR-to-FIR galaxy SED.

We have also considered HerMES observations (Oliver et al. 2012), which entirely cover the PEP area with SPIRE maps at 250, 350 and 500 µm. Also the SPIRE counterparts of our PACS-selected sample have been restricted to the common MUSIC area. In HerMES, SPIRE flux densities are taken by following the prior source extraction presented by Roseboom et al. (2012), based on Spitzer-MIPS 24 µm positions as a prior, and cross-matched to the blind PACS catalogue. In the GOODS-S area, we have also considered Spitzer-IRS (Infrared Spectrograph) 16 µm detections (Teplitz et al. 2011), which have been cross-correlated to the available 24 µm positions through a nearest-neighbour match.

In the COSMOS field, the sky area explored by Herschel and observations at shorter wavelengths is around 2 deg<fnref n="2" />. PACS 160 µm sources have been associated with the Spitzer-MIPS 24 µm ones from Le Floc'h et al. (2009) and with the IRAC-selected catalogue from Ilbert et al. (2010), which is already matched to UV, optical and NIR photometry. By removing PEP sources within flagged areas of the optical-to-NIR catalogues, we are left with 4118 sources at 160 µm identified at shorter wavelengths. Around 96.5 per cent of them have an MIPS 24 µm detection in the common area. The

<footnote n="2">²The MUSIC catalogue can be retrieved at http://lbc.oa-roma.inaf.it/goods/goods.php</footnote>

<footnote n="1">¹http://www2011.mpe.mpg.de/ir/Research/PEP/</footnote>
fraction of PACS-selected sources having at least one SPIRE counterpart ($>3\sigma$) reaches 84 per cent in the GOODS-S field (within the MUSIC area) and 86 per cent in the COSMOS field.

In both fields, we checked some properties of PACS sources without a SPIRE counterpart and we ensured that the integrated IR (8–1000 µm) luminosities<fnref n="3" /> and redshifts were consistent within $1\sigma$ uncertainty with those of PACS sources having a SPIRE counterpart.

Overall, we collect photometry at a maximum of 22 (GOODS-S) and 20 (COSMOS) different wavebands, from UV to FIR wavelengths. The numbers of PEP-selected objects fulfilling the optical identification criterion are 494 and 4118, in the GOODS-S and COSMOS, respectively.

## 2.2 Redshift distribution

In GOODS-S, we rely on a large spectroscopic data base,<fnref n="4" /> including more than 3000 sources (Cristiani et al. 2000; Croom et al. 2001; Bunker et al. 2003; Dickinson et al. 2004; Stanway et al. 2004a,b; Strolger et al. 2004; Szokoly et al. 2004; van der Wel et al. 2004; Doherty et al. 2005; Le Fèvre et al. 2005; Mignoli et al. 2005; Vanzella et al. 2008; Popesso et al. 2009; Santini et al. 2009; Balestra et al. 2010; Cooper et al. 2012; Kurk et al. 2013). We used this to enrich the existing spectroscopic data base of the MUSIC catalogue (Grazian et al. 2006; Santini et al. 2009) replacing MUSIC photometric redshifts with the up-to-date spectroscopic ones in the case of ‘secure’ redshift determination (according to the scale adopted by Balestra et al. 2010). In summary, the GOODS-S PEP sample includes 494 sources with 100 per cent redshift completeness (71 per cent spectroscopic and 29 per cent photometric).

A rich redshift data set is also available for COSMOS, including both photometric (Ilbert et al. 2010) and spectroscopic (Lilly et al. 2009; Trump et al. 2009) redshifts.

Recently, Berta et al. (2011) extended the previous photometric redshift data set for most of the remaining PEP sample, obtaining an overall redshift completeness as high as 93 per cent (about 40 per cent spectroscopic). By using the latter data set, we gather a number of 3849 PEP sources with a redshift estimate.

The photometric redshift uncertainties have been evaluated in both fields from previous works. In GOODS-S, Grazian et al. (2006) found an average absolute scatter $\langle|\Delta z/(1+z)|\rangle = 0.045$ over the whole redshift range $0 < z < 6$. In COSMOS, Ilbert et al. (2010) evaluated the scatter on the photo-$z$’s to range from 0.008 (for $I < 22.5$) to 0.053 (for $24 < I < 25$), while the more up-to-date photometric redshifts provided by Berta et al. (2011) have a median absolute deviation (MAD<fnref n="5" />) of around 0.01. The expected fraction of outliers, defined as objects having $\langle|\Delta z/(1+z)|\rangle \ge 0.2$, is as high as ~2 per cent.

The final sample of GOODS-S (GS) and COSMOS (C) sources with redshifts includes 4343 objects (494 for GS and 3849 for C), which are selected at 160 µm. The global redshift distributions of PACS 160 µm selected sources are presented in Fig. 1 (top) for both the PEP fields. As expected, because of the higher flux limit, the COSMOS redshift distribution has a lower median redshift value ($\langle z \rangle \sim 0.7$) than the GOODS-S ($\langle z \rangle \sim 1$), although it also extends to high redshift ($z > 3$). To highlight the range of IR luminosities sampled

<figure><img src="image_2.png" /></figure>

**Figure 1.** Top panel: redshift distribution of the 160 µm PEP sources identified in optical bands, in GOODS-S (red dashed) and COSMOS (blue dashed). Note that the scale of the y-axis is logarithmic. Bottom panel: IR luminosity limit at 160 µm (rest-frame) as a function of redshift. Open circles represent the luminosity limits for GOODS-S (red) and COSMOS (blue) sources. The solid curves show the minimum luminosity (corresponding to the limiting flux) versus redshift for the two fields (GOODS-S in red, COSMOS in blue), obtained for the empirical SED of the source IRAS 20551–4250.

at different redshifts, we show in Fig. 1 (bottom) the minimum 160 µm luminosity versus redshift of the (rest-frame) 160 µm luminosity. Red and blue circles represent GOODS-S and COSMOS sources, respectively. The solid curves represent the minimum luminosity (corresponding to the limiting flux) versus redshift for the two fields (GOODS-S in red, COSMOS in blue). This curve is purely illustrative and is taken from the SED template of IRAS 20551–4250, which is a local star-forming galaxy hosting an AGN (see Polletta et al. 2007).

## 3 FITTING BROAD-BAND SEDs

In this section, we present the SED-fitting analysis performed for the 4343 PACS sources with optical-to-MIR counterparts and either spectroscopic or photometric redshift. Each observed SED has been fitted making use of the MAGPHYS code<fnref n="6" /> (da Cunha, Charlot & Elbaz 2008) and a modified version of this code presented by Berta et al. (2013).

### 3.1 SED fitting with MAGPHYS

MAGPHYS is a public code which uses physically motivated templates to reproduce the observed SEDs from the ultraviolet to sub-mm wavelengths. It has been developed especially to fit star-forming galaxy SEDs, so that it is well suited for our purposes (see e.g. Smith

<footnote n="3">³ Calculated through SED fitting. See Section 3 for details.</footnote>

<footnote n="4">⁴ Publicly available at http://www.eso.org/sci/activities/garching/projects/goods/MasterSpectroscopy.html</footnote>

<footnote n="5">⁵ MAD($x$) = median($|x - \text{median}(x)|$), where $x$ represents the difference between spectroscopic and photometric redshift.</footnote>

<footnote n="6">⁶ MAGPHYS can be retrieved at http://www.iap.fr/magphys/magphys/MAGPHYS.html</footnote>
et al. 2012 for an extensive application of MAGPHYS on Herschel data). Indeed, the basic recipe that this code relies on is a self-consistency between the energy originating from stars, which is partially absorbed by dust, and the IR light due to re-emission by the dust itself. Briefly, the stellar and dust components are computed as follows (see da Cunha et al. 2008 for details).

(i) The Bruzual & Charlot (2003) library of templates is used to compute the integrated light produced by stars in galaxies. These models predict the spectral evolution of stellar populations at $\lambda \le 10$ µm and ages between $10^5$ and $10^{10}$ yr, starting from a Chabrier (2003) IMF and combining parameters representative of different metallicities, SFH and dust contents.

(ii) Each stellar template is modified by the effect of dust extinction, which is parametrized as the cumulative attenuation produced by interstellar medium (ISM) and dusty molecular clouds, according to the angle-averaged model of Charlot & Fall (2000).

(iii) The global dust emission is ascribed to three main constituents of interstellar dust: polycyclic aromatic hydrocarbons (PAHs, in the range 3–12 µm); MIR continuum from hot (130–250 K) dust grains; and thermal emission, both from warm (30–60 K) and cold (15–25 K) dust grains in thermal equilibrium, according to da Cunha et al. (2008) prescriptions.

The analysis of each single observed SED with MAGPHYS starts by building up a large library of templates at the same redshift of the source. Given the set of parameters considered as well as the wide range of values adopted by each of them, in the end 50 000 stellar and 50 000 independent dust emission models are computed for each source. Only about 30 per cent of the IR templates can be combined with each stellar model, because of the underlying assumption that all the absorbed UV and optical stellar radiation has to be accounted for in the integrated IR emission. As a result of such selection, typically $10^9$ different self-consistent (i.e. stellar+dust) model combinations are created to reproduce each observed SED from the UV to the FIR. Then, each model combination is compared to the observed SED, and the stellar+dust model normalization is rigidly scaled to fit the data as well as possible. Moreover, a Bayesian approach is implemented to build a marginalized likelihood distribution for each parameter, which is traced by the distribution of values of such a parameter across the whole library of models. This statistical analysis accounts for the degeneracy level in the parameter space, which occurs when a relatively small set of data points is compared with such a large range of models. The resulting probability distribution function (PDF) allows the user to obtain reliable confidence ranges for parameter estimates.

## 3.2 Fitting with MAGPHYS+AGN

MAGPHYS does not consider a potential AGN contribution, but only star formation processes are involved in building the SEDs. Since we are interested in probing the possible AGN contribution to the observed SEDs, we make use of a modified version of the code (Berta et al. 2013), which accounts also for a possible nuclear heating source.

The code by Berta et al. (2013) includes the AGN component and differs from the original MAGPHYS in the model normalization approach. The stellar+dust normalization is not scaled to fit the data, but is allowed to float within a range of normalization values (the same range adopted by the original version of MAGPHYS), to allow for the additional AGN normalization when fitting the observed SED. This approach results in a three-component simultaneous fit, using an overall $\chi^2$ minimization.

The AGN library of templates (Feltre et al. 2012; see also Fritz, Franceschini & Hatziminaoglou 2006) includes around 2400 models, each one analysed for 10 different lines of sight $\theta$, uniformly distributed in $\theta$ from edge-on to face-on. This set of templates is widely used in the literature, since it is physically motivated and suitable in reproducing a wide variety of AGN features, such as hot circumnuclear dust in MIR dominated SEDs (e.g. Mrk3; see Mullaney et al. 2011) and the silicate features at 9.7 and 18 µm (see Feltre et al. 2012).

Each nuclear template includes both the contribution from the central engine (i.e. accretion disc) and the reprocessed emission from a smooth surrounding dusty torus.

The central heating source is modelled with a broken power law $\lambda L_{\lambda} \propto \lambda^{\alpha}$, with slopes $\alpha_1 = 1.2$, $\alpha_2 = 0$ and $\alpha_3 = -1$, in the wavelength ranges $\lambda_1 = 0.001-0.03$ µm (Hubeny et al. 2000), $\lambda_2 = 0.03-0.125$ µm (Zheng et al. 1997) and $\lambda_3 = 0.125-20$ µm (Hatziminaoglou et al. 2008), respectively. At longer wavelengths, the emission scales as a Planck spectrum in the Rayleigh–Jeans regime. The implementation of the circumnuclear dust component relies on the numerical solution of the radiative transfer equation for a smooth dusty structure.

Each torus model is univocally identified by a set of six different parameters: the outer-to-inner radius ratio $R_{\text{out}}/R_{\text{inn}}$; the opening angle $\Theta$; the equatorial optical depth at 9.7 µm $\tau_{0.7}$; the radial and height slopes of the density profile ($\beta, \gamma$); and the viewing angle $\theta$ of the line of sight.

During the AGN fitting procedure, 100 different AGN templates are randomly selected for each source and for each stellar+dust realization. In this way, it is possible to optimize the random sampling by exploring a wider range of the AGN parameter space for each source.

By solving the radiative transfer equation for a smooth dusty structure irradiated by the accretion disc, a bolometric correction (BC) value is available for each AGN template. Such information allows us to calculate the AGN bolometric luminosity ($L_{\text{bol,AGN}}$) directly from the IR (1–1000 µm rest-frame) luminosity of the best-fitting AGN model, according to the following expression:

$$k_{\text{bol}} = \frac{L_{\text{bol,AGN}}}{L_{1-1000}^{\text{agn}}}. \qquad (1)$$

Our $L_{\text{bol,AGN}}$ estimates are derived without double counting the dust-reprocessed contribution arising from the accretion disc, since the radiative transfer recipe performs a self-consistent re-distribution of the input energy at all wavelengths, according to the adopted set of templates.

### 3.2.1 Reduced library of templates

The set of possible model combinations included in the original version of MAGPHYS is very large (around $10^9$ possible configurations, see Section 3.1), as well as the set of physical parameters that the stellar+dust emission models rely on. The range of values spanned by each parameter is physically reasonable, but the computational time taken by the code strongly increases when the AGN component is added to the SED-fitting procedure. For these reasons, we take advantage of a reduced grid of templates, as described by Berta et al. (2013).

While the full grid of MAGPHYS templates includes 50 000 × 50 000 model configurations, the reduced grid has been built up through a random selection of 1000 stellar templates for each source and 1000 dust emission templates for each stellar template. This selection technique permits to save a factor of 2500 in time
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison between the values obtained with the full (x-axes) and reduced (y-axes) MAGPHYS library of templates for three different integrated parameters: stellar mass, IR (8–1000 µm) luminosity due to star formation and SFR. The full grid uses all the available templates of galaxy SEDs encoded in MAGPHYS, whereas the reduced grid accounts for a small fraction of them resulting from a random selection (see the text for details). Red triangles denote the best-fitting values of each parameter as taken from the GOODS-S sample. The black solid line shows the bisector, while the dashed lines represent the ±1σ confidence regions. The resulting 1σ uncertainties related to these parameters are 0.13, 0.11 and 0.17 dex, respectively.</figcaption></figure>

(2 h per object, instead of 5000), giving a total of ~3 × 10⁵ self-consistent templates for each source (after the 30 per cent cut required by the energy balance argument). Such selection uniformly spans the parameter space covered by the original library, without introducing any significant bias in the adopted grid. In Fig. 2, we compare the values of some integrated physical quantities (SFR, IR luminosity, etc.) obtained using the full and the reduced grid of MAGPHYS. The plots show that mostly the integrated physical parameters are consistent among the full and reduced grids.

As described in Section 3.2, each AGN template is defined by the combination of six physical parameters. We limit the torus library of templates by restricting the multidimensional volume of the parameter space. In particular, we removed the largest optical depth value at 9.7 µm (see Pier & Krolik 1992) and the most extended geometries ($R_{\text{out}}/R_{\text{inn}} > 100$), as no evidence for such large structures has yet been found (Jaffe et al. 2004; Tristram et al. 2007, 2009).

After this cut, the torus grid includes 504 templates, each one computed at 10 different lines of sight. As mentioned before, a random selection of 100/5040 AGN configurations is performed for each source and stellar+dust model combination. This reduced library spans several values of the ratio $R_{\text{out}}/R_{\text{inn}}$ (10, 30, 60 and 100), the opening angle $\Theta$ (40°, 100° and 140°), the equatorial optical depth $\tau_{9.7}$ (0.1, 0.3, 0.6, 1, 2, 3 and 6) and the slopes of the density profile (-1, -0.5 and 0 for $\beta$; 0 and 6 for $\gamma$). The resulting library turns out to be similar to the one used by Pozzi et al. (2012) and Hatziminaoglou et al. (2008), which we refer to for a detailed description of the degeneracy level drawn by torus parameters.

Since the only crucial parameter in the following analysis is the AGN bolometric correction ($k_{\text{bol}}$), used to convert the total IR AGN emission into AGN bolometric luminosity, we show in Fig. 3 how much our restrictions modify the shape of the distribution of the BC. In particular, the reduced grid of AGN templates reaches $k_{\text{bol}}$ values around 60, while the full grid samples up to $k_{\text{bol}} \sim 150$. The largest values are mostly representative of the most obscured AGN ($\tau_{9.7} = 10$). Despite the existing cut in $k_{\text{bol}}$, we point out that the missed values (i.e. $k_{\text{bol}} > 60$) constitute a negligible number of AGN templates with respect to the full grid (~0.3 per cent).

In addition, we checked the stability of our $k_{\text{bol}}$ estimates with increasing number of galaxy templates. We have run the MAGPHYS+AGN SED fitting on a subsample of 50 Herschel sources randomly selected from the GOODS-S field. We compared the $k_{\text{bol}}$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** BC (1–1000 µm, rest-frame) distributions to turn the IR AGN luminosity into bolometric luminosity. The red area traces the full grid of the library by Feltre et al. (2012) and Fritz et al. (2006). The green area represents the reduced library that is adopted in this paper.</figcaption></figure>

estimates taken from the run with 1000 × 1000 MAGPHYS templates with those obtained from the run with 5000 × 5000 templates. We found no systematics between the two SED-fitting analyses, with a 1σ scatter equal to 0.07 dex. Though this test does not account for the full grid of available galaxy templates, we can at least say that the AGN BC does not change significantly with increasing number of galaxy templates.

### 3.2.2 AGN bolometric luminosity and bolometric corrections

For every observed SED, a PDF is created for each parameter to trace its confidence range. As mentioned before, among the physical parameters returned by the code, we have focused on the AGN bolometric luminosity only (computed by accounting for the rest-frame 1–1000 µm AGN emission and adopting its related BC).

We stress that each $L_{\text{bol,AGN}}$ value is not properly a bolometric luminosity, since the input (i.e. accretion) energy of the central engine is limited to the range $10^{-3}-10^3$ µm (i.e. $10^{-1}-10^{-7}$ keV). This means that the X-ray emission is supposed to be negligible in terms of input energy. Such an assumption relies on the fact that
<figure><img src="image_5.png" /><figcaption>Figure 4. Distributions of AGN BC. The regions traced by solid lines represent various trends of $k_{bol}$ computed in different bands by Hopkins et al. (2007). Black circles (with ±1σ uncertainties) come from the sample analysed in this paper, only including sources with relevant signatures of nuclear activity (see Section 3.3 for details).</figcaption></figure>

large $k_{bol}$ values have been found in hard X-rays, of the order of 20–30 (Marconi et al. 2004; Risaliti & Elvis 2004; Hopkins et al. 2007; Pozzi et al. 2007; Vasudevan & Fabian 2009; Lusso et al. 2012), increasing as a function of $L_{bol,AGN}$.

In Fig. 4, we show a comparison between distributions of BC computed in various bands. On the one hand, Hopkins et al. (2007) derived their $k_{bol}$ values by assuming different amounts of obscuration applied to an average, unobscured, type 1 AGN SED (Richards et al. 2006). They found $k_{bol}$ to depend on the input AGN luminosity, consistent with other works (e.g. Marconi et al. 2004; Lusso et al. 2012). On the other hand, we derive the distribution of $k_{bol}$ by integrating the best-fitting AGN model of each source (if it fulfils the criterion imposed by the Fisher test, as explained in Section 3.3) over the whole IR domain (1–1000 µm, rest-frame). Each black circle represents the median $k_{bol}$ value on different luminosity bins (=0.3 dex wide). Error bars correspond to ±1σ uncertainties and are centred on the median value of the respective bin. We find that $k_{bol}$ does not change significantly with increasing AGN bolometric luminosity $L_{bol,AGN}$. This result could suggest that the fraction of the input energy coming from the central engine which is absorbed by the obscuring structure and re-emitted in the IR does not change significantly as a function of $L_{bol,AGN}$.

The average BC derived from the IR are much lower than those derived from X-ray measurements (see Fig. 4). Indeed, the distribution of 1–1000 µm BC has a median value of ~3.8 over the whole range of $L_{bol,AGN}$ covered by our sample, ranging on average from 2 to 5.

### 3.3 Testing the AGN incidence

In order to evaluate the impact of an additional AGN component on a ‘pure’ star-forming galaxy SED, for each source we run both MAGPHYS and MAGPHYS+AGN. As mentioned in Section 3.2.1, a reduced grid has been adopted in both fitting procedures. First, we run the code version including the AGN, and during the random selection of stellar and dust templates, their template identification numbers are registered. Once the SED fitting with MAGPHYS+AGN is

accomplished, we run the MAGPHYS code for each source by using the same set of stellar+dust templates previously extracted for the same source. This implementation allows us to test if the fit provided by similar star formation templates improves when accounting for an additional nuclear emission component.

In Fig. 5, three different examples of best fits are shown as a result of this SED decomposition. The integrated emission has been decomposed into three emission contributions, according to the recipe previously described. The expected stellar light (green line) is partially absorbed by a mixture of dust grains and re-distributed in the MIR/FIR domain (blue line). The AGN emission is denoted by the red curve. The respective $\chi^2$ values are then used to estimate the significance of the improvement of the fit occurring when a nuclear component is added.

Around 98 per cent of the PEP sample fitted with the MAGPHYS+AGN code includes an AGN component in the best-fitting model. However, the addition of the AGN component does not always significantly improve the best fit. For instance, the last source shown in Fig. 5 turns out to exhibit negligible signatures of AGN activity over the entire wavelength range. To pinpoint the significance of the AGN component on the global SED, we apply a Fisher test (F-test) to the two SED-fitting runs, as presented by Bevington & Robinson (2003):

$$F_{\text{test}} = \frac{\chi^2_{\text{noagn}} - \chi^2_{\text{agn}}}{\chi^2_{\text{v,agn}}} \geq F_{\text{threshold}}(\text{CL}) \quad (2)$$

with $\chi^2_{\text{noagn}}$ and $\chi^2_{\text{agn}}$ being the $\chi^2$ values of the best fit without and with the AGN, respectively, whereas $\chi^2_{\text{v,agn}}$ is the reduced $\chi^2$ value referred to the best fit with an AGN component. If this ratio exceeds a given threshold, then the source of interest is included in the following analysis, otherwise it is ruled out and a null AGN contribution is assumed for its specific SED. The threshold value depends on the number of degrees of freedom as well as on the required confidence level (CL).

When carrying out this analysis at a 99 per cent CL, in the GOODS-S field we find that the AGN detection rate is 48 per cent, while in COSMOS it is around 36 per cent. By evaluating the binomial confidence limits following Gehrels (1986) and Gerke et al. (2007), such a difference turns out to be significant at the >3σ level. To make a more coherent comparison between these percentages, we checked the incidence of some possible selection effects. We find that the AGN detection rate depends on both the observed 160 µm flux density and the presence of the 16 µm detection, which is missing in COSMOS. As a sanity check, we run MAGPHYS and MAGPHYS+AGN on the entire GOODS-S sample without using the 16 µm data during the SED fitting. To compare the resulting AGN detection rate with that obtained in COSMOS, we selected only sources above the COSMOS 160 µm flux density limit ($S_{160} > 10.2$ mJy). The percentage of AGN hosts in GOODS-S obtained in this case decreases to 40 ± 4 per cent and is now consistent with the COSMOS one within 1σ uncertainty. The overall fraction of >99 per cent significant AGN is 37 per cent, while it increases to 45 and 55 per cent at the >95 and >90 per cent CL, respectively. We label as AGN only those sources satisfying the F-test at the >99 per cent CL, as this threshold is fairly conservative and reliable in terms of AGN contribution to the observed SED. In Fig. 5, the upper and middle panels represent sources with a >99 per cent significant AGN contribution, while that shown in the bottom panel has a less significant AGN contribution, hence classified as a ‘purely star-forming galaxy’.

<footnote n="7">⁷ The terms ‘accretion’ and ‘bolometric’ AGN luminosity adopted throughout the paper are assumed to have the same physical meaning.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Three examples of SED decomposition performed without AGN (left-hand panels) and with AGN (central panels). The corresponding PDF of the AGN bolometric luminosity is shown in the right-hand panels. The red circles represent the observed data points. The green line shows the integrated unextincted emission originating from stars. The blue solid line represents the star formation contribution for dust absorption, partially re-distributed across the MIR/FIR range in a self-consistent way. The red line reproduces the AGN contribution and incorporates both the accretion disc and the torus emission, according to the adopted library of templates. The sources shown in the upper and middle panels are classified as AGN at ≥99 per cent confidence, while one example of source with less significant AGN contribution is shown in the bottom panel.</figcaption></figure>

We have verified that the star-forming galaxy SEDs adopted in the code were physically motivated. For instance, we found that our best-fitting galaxy SEDs well reproduce the ‘IR8’ relation, presented by Elbaz et al. (2011), within ±1σ uncertainty. This finding also removes possible systematics related to the treatment of PAH emission component with respect to the total IR luminosity. We also checked the potential degeneracy between AGN and starburst emission components when fitting the rest-frame MIR (3 < λ/μm < 8) data points. Indeed, we note that the ISM emission in the MAGPHYS star-forming SEDs includes dust and PAH emission down to 3 μm (see da Cunha et al. 2008). An excess above the pure stellar continuum will thus not be automatically ascribed to AGN. Besides, given the redshift distribution of our sources in Fig. 1, the rest MIR continuum is sampled at least by two filters in almost all observed SEDs (four IRAC bands, MIPS 24 μm and IRS 16 μm in GOODS-S). Such a rich spectral coverage allows us to better constrain the rest MIR emission through the available templates.

### 3.3.1 AGN detection rate and evolution with $L_{1-1000}$

The sources fulfilling the F-test at the >99 per cent CL are 1607/4343 (37 per cent). Our finding settles in an open debate concerning the fraction of AGN which are detected through independent selection methods. In the following, we check the consistency between our results and the most recent estimates inferred from the analysis of different IR galaxy samples.

In the local Universe, Yuan, Kewley & Sanders (2010) found the percentage of sources hosting an optically identified AGN to depend on the total IR luminosity (8–1000 μm).<fnref n="8" /> This fraction ranges from 20–30 per cent up to ~60 per cent in local luminous ($L_{IR} \sim 10^{11} L_{\odot}$) and ultra-luminous ($L_{IR} \sim 10^{12} L_{\odot}$) infrared galaxies (LIRGs and

<footnote n="8">⁸ For historical reasons, throughout the paper the expression ‘IR luminosity’ (or $L_{IR}$) will refer to the luminosity integrated over the spectral range 8–1000 μm. If calculated in a different range (e.g. 1–1000 μm), it will be explicit.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** AGN detection rate in the 160 µm PEP-selected sample, as a function of the total IR (1–1000 µm) luminosity and redshift. In each redshift bin, the sample has been split into similarly populated luminosity bins. Error bars along the y-axis correspond to $1\sigma$ confidence limits, while downward arrows show $2\sigma$ upper limits. The dashed line shows a linear fit to the data of the second redshift bin to highlight the evolution of the AGN detection rate with redshift.</figcaption></figure>

ULIRGs), respectively. Consistent results have been obtained by Imanishi et al. (2010) and Alonso-Herrero et al. (2012) through Akari-IRC (InfraRed Camera) and Spitzer-IRS spectroscopy of local LIRGs and ULIRGs. These authors independently found that the AGN detection rate increases with increasing $L_{IR}$, from 25 per cent (LIRGs) up to 70 per cent (brightest ULIRGs), strengthening the previous results of Lutz et al. (1998) based on Infrared Space Observatory data (from ~15 per cent at $L_{IR} < 2 \times 10^{12} L_{\odot}$ up to 50 per cent for brighter ULIRGs). The spectral decomposition has also been useful to effectively identify AGN: through detection of IR fine-structure lines, PAH equivalent widths and MIR continuum, several works (e.g. Sajina et al. 2009, 2012; Nardini et al. 2010) claimed the presence of a significant AGN in 60–70 per cent of local and distant ULIRGs, as well as evident signatures of coexisting star formation.

At higher redshifts, the situation is more controversial. Olsen et al. (2013) analysed a mass-selected sample at $z \sim 2$ within the Chandra Deep Field-South and ended up with a high fraction (43–65 per cent) of star-forming galaxies likely hosting an AGN, according to the X-ray classification taken from Xue et al. (2011). Pozzi et al. (2012) studied a sample of 24 ULIRGs at $z \sim 2$ (see Fadda et al. 2010) within the GOODS-S field and selected to have faint 24 µm flux densities ($0.14 < S_{24\mu m}/mJy < 0.55$). They ended up with a smaller fraction (~35 per cent) of ULIRGs featuring signatures of AGN activity and claimed that their IR luminosity emission is dominated by starburst processes.

In this work, the sample has been split into six redshift bins ($0.1 \le z < 0.3$; $0.3 \le z < 0.7$; $0.7 \le z < 1.2$; $1.2 \le z < 1.8$; $1.8 \le z < 2.5$ and $2.5 \le z \le 3.8$) and into different IR luminosities. We check for the evolution of the AGN detection rate with total (AGN+starburst) 1–1000 µm luminosity $L_{1-1000}$ at different redshifts, as shown in Fig. 6. Error bars along the y-axis show the $\pm 1\sigma$ uncertainty (following Gehrels 1986 and Gerke et al. 2007), while the error bars along the x-axis set the bin width. Downward arrows represent $2\sigma$ upper limits to the expected AGN detection rate when only ‘inactive’ galaxies are detected in that bin. We find an overall increase of the AGN detection rate as a function of $L_{1-1000}$, except

in the lowest redshift bin. At $z < 0.3$, the AGN hosts are mainly ‘normal’ ($L_{1-1000} \le 10^{11} L_{\odot}$) star-forming galaxies. At these redshifts, the average AGN detection rate is ~5–10 per cent and we do not find a clear trend with 1–1000 µm luminosity. This is likely due to the fact that the comoving volume covered by COSMOS is not large enough to detect sources with $L_{1-1000} \ge 10^{12} L_{\odot}$.

At higher redshift ($0.3 < z < 3.8$), the AGN detection rate does increase with increasing $L_{1-1000}$, which is consistent with results from previous studies. This finding might suggest that the probability of finding AGN in FIR-selected galaxies increases as a function of $L_{1-1000}$. Besides the trend with $L_{1-1000}$, we point out the evidence for redshift evolution (see Fig. 6). Indeed, the AGN detection rate of $L_{1-1000} \sim 10^{11} L_{\odot}$ does increase from 5 to 30 per cent up to $z \sim 1$; galaxies with $10^{11} \le L_{1-1000}/L_{\odot} \le 10^{12}$ host AGN activity with a probability rising from 25 per cent at $z \sim 0.5$ to 60 per cent at $z \sim 2$. Brighter ($L_{1-1000} \ge 10^{12} L_{\odot}$) IR galaxies show an AGN detection rate rising from 55 per cent at $z \sim 1$ to 70–80 per cent at $z \sim 2–3$.

We also checked the potential dependence of this finding on selection effects. Indeed, our selection at 160 µm is certainly sensitive to warmer dust with increasing redshift, which might favour the detection of galaxies hosting AGN activity in more distant sources. However, given the redshift distribution of our Herschel-selected sample in the GOODS-S and COSMOS fields, we detected 35 objects only at $3 < z < 3.8$, where our selection wavelength spans the spectral range $35 < \lambda < 45$ µm rest frame. Moreover, the AGN templates adopted in this work (from Feltre et al. 2012; Fritz et al. 2006) typically peak around 15–20 µm (rest-frame), going down at longer wavelengths. We examined for each source the fractional energy contribution at 160 µm (observed-frame) due to AGN emission, as returned from the MAGPHYS+AGN SED-fitting analysis. We found the selection wavelength being poorly related to AGN activity, since its energy contribution to the total energy emitted at 160 µm is negligible (<1 per cent) with respect to that ascribed to star formation, apart from a few exceptional cases (about 20 sources) that are generally characterized by an AGN-dominated SED in the MIR. This check makes us reasonably confident that the selection at 160 µm does not affect the overall increase of the AGN detection rate with redshift.

### 3.3.2 Comparison with MIR colour–colour AGN selection

In order to test the robustness of the method that we apply to estimate the AGN contribution in FIR-selected galaxies, we compare our findings with those obtained from an independent AGN diagnostic, specifically the MIR colour–colour selection. This criterion was first introduced by Lacy et al. (2004, 2007) by making use of Spitzer-IRAC colours to isolate the AGN candidates, both obscured and unobscured. However, as already stated by Lacy et al. (2007), a sample of so-called AGN candidates might suffer from a significant contamination from ‘pure’ galaxy SEDs with similar MIR colours. More recently, Donley et al. (2012) revisited the previous method to provide a more reliable AGN selection criterion, although it was less complete at low X-ray luminosities ($L_{[0.5–8]} < 10^{44} \text{ erg s}^{-1}$). In Fig. 7, we compare our classification with that derived from these MIR colour–colour diagnostics. We report the comparison with Donley et al. (2012) in the upper panels and with Lacy et al. (2007) in the bottom ones. The left- and right-hand panels refer to our AGN and galaxy samples, respectively. Red asterisks mark the population of sources ‘inside’ the corresponding wedge, whereas blue asterisks represent the sources in the PEP sample not satisfying the colour–colour criterion.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison with AGN colour–colour selection criteria. Upper panels: comparison with Donley et al. (2012). Bottom panels: comparison with Lacy et al. (2007). The left- and right-hand panels refer to PEP sources classified as AGN and galaxies, respectively. Sources satisfying the corresponding criteria are marked with red asterisks, otherwise with blue asterisks.</figcaption></figure>

As shown in Fig. 7, some PEP sources, although inside the black solid wedge of Donley et al. (2012), are labelled as ‘outside’. This is why, in addition to the requirement of being inside the wedge, the IRAC colours of each source have to fulfil further conditions on IRAC fluxes at the same time to satisfy the Donley criteria.<fnref n="9" />

The colour–colour wedge defined by the Donley criteria (black solid line) includes ~17 per cent of our AGN sample and only ~0.8 per cent of our galaxy sample. According to our classification, we confirm that Donley et al. (2012) selection is very efficient in avoiding contamination from purely star-forming galaxies. The overlapping fraction of AGN obviously increases when considering the wider Lacy et al. (2007) region: we find that 73 per cent of our PEP AGN and 23 per cent of our PEP galaxies are inside the wedge.

We checked the average integrated properties of PEP AGN inside or outside both colour–colour AGN selection criteria. As expected, objects satisfying the Donley criterion are mostly AGN-dominated systems in the MIR (both obscured and unobscured) and populate the highest luminosity tail of different luminosity distributions. Indeed, typical values for the 1–1000 µm AGN luminosity reach ~10¹² L⊙, as well as for the 1–1000 µm luminosity due to star formation. Moving towards the area delimited by the Lacy et al. (2007) criteria and outside the Donley et al. (2012) selection, the typical properties become less and less extreme, with Seyfert-like

SEDs (i.e. dominated by the host galaxy light in the IR) being more common than quasi-stellar object (QSO)-like ones.

In summary, we find the colour–colour criterion by Lacy et al. (2007) to be in reasonably good agreement with our classification and to represent an acceptable compromise between reliability and completeness. In addition, the Lacy wedge effectively rules out sources with 1–1000 µm SF luminosity $L_{1-1000} < 10^{11}$ L⊙, where the AGN contribution to the IR becomes negligible (~ few per cent) and also the probability of picking up AGN is generally lower (see Section 3.3.1 and Fig. 6).

PEP AGN as classified through SED decomposition occur in different regions with respect to PEP galaxies in the MIR colour–colour diagram. This, together with the percentage of our PEP AGN and galaxies fitting with the Donley and Lacy criteria discussed above, further supports the robustness of our classification.

### 3.3.3 Specific SFR versus redshift

Our SED-fitting analysis allows us to get both stellar mass ($M_*$) and SFR for each source. Therefore, to derive the specific star formation rate (sSFR) for the entire Herschel-selected sample and to infer its evolution with redshift.

Previous studies have highlighted the presence of a tight correlation between stellar mass and SFR in the local Universe (e.g. Peng et al. 2010, 2012). These parameters trace the well-known ‘main-sequence’ (MS) of star-forming galaxies. Such a relation has been extended to $z \sim 1$ (Elbaz et al. 2007) and up to $z \sim 2$ (Daddi et al. 2007; Rodighiero et al. 2011) and $z \sim 3$ (Magdis et al. 2010). In particular, Rodighiero et al. (2011) argued the presence of a bi-modal regime in star formation between galaxies lying on the MS and those on the off-sequence (i.e. $>0.6$ dex above the MS), with the former fuelled through steady star formation events, whereas the latter are experiencing enhanced and starbursting star formation.

In this subsection, we compare the mean evolutionary trend of the sSFR with redshift with independent findings from the literature. The stellar mass is inferred from the best-fitting template of the MAGPHYS+AGN SED fitting. The SFR has been derived by converting the IR (rest-frame 8–1000 µm) luminosity through a standard Kennicutt (1998) relation, rescaled to a Chabrier IMF. In the case of ‘active’ galaxy, the AGN contribution has been subtracted from the total IR luminosity, otherwise the IR emission is supposed to be due to star formation only.

In Fig. 8, we show the overall evolution of the sSFR with redshift for sources with $10 < M_* < 11$. The joint (GOODS-S and COSMOS) sample has been split into six different redshift bins, as previously done in Fig. 6. Black dots represent each source of our sample, while magenta triangles show the median sSFR in each redshift bin. Error bars cover the $\pm 1\sigma$ population in each z-bin along the vertical axis and set the bin width along the horizontal axis. In this plot, only sources above the nominal SFR limit have been considered. The cut in SFR at a given redshift bin corresponds to the minimum SFR value spanned by the SFR–z distribution at the upper bound of the same z-bin. This cut ensures to get a complete sample in SFR in each z-bin.

For comparison, in Fig. 8 we report other independent estimates of the sSFR as a function of redshift. Rodighiero et al. (2010) studied a Herschel-selected sample in the GOODS-N field, while Rodighiero et al. (2011) extended the previous analysis with Herschel data to the COSMOS and the GOODS-S field. Daddi et al. (2007) and Dunne et al. (2009) analysed a sample of k-selected galaxies, in the GOODS-S and in the Ultra Deep Survey,

<footnote n="9">⁹ Given the monochromatic fluxes ($S_\lambda$) at 3.6, 4.5, 5.8 and 8.0 µm, the criteria established by Donley et al. (2012) set the following additional conditions: $S_{3.6} < S_{4.5} < S_{5.8} < S_{8.0}$ at the same time, which make the black solid wedge a multidimensional region.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8. sSFR as a function of redshift in the stellar mass range $10 < M_* < 11$. Black dots represent the sSFR obtained through our SED decomposition for the joint (GOODS-S and COSMOS) sample. Magenta triangles set the median sSFR in each redshift bin. Horizontal error bars set the bin width; vertical error bars enclose the ±1σ population of sources in each redshift bin. Various estimates from the literature, complemented with their respective error bars, are reported for comparison in the same stellar mass range.</figcaption></figure>

respectively. Elbaz et al. (2007) carried on an IRAC-based work (3.6 and 4.5 µm selection) in the GOODS-S. Magdis et al. (2010) investigated a sample of IRAC-detected Lyman break galaxies in the Hubble Deep Field-North.

We find our trend to be in fair agreement with the other independent estimates, within 1σ uncertainty. However, a more detailed analysis of the AGN contribution in the $M_*$-SFR plane, for both detected and undetected X-ray sources, is postponed to a forthcoming paper.

# 4 AGN BOLOMETRIC LF

In this section, we investigate the evolution of the AGN bolometric LF over the whole joint GOODS-S and COSMOS sample. Globally we collect 4343 objects, but just the fraction showing a substantial AGN contribution has been considered for the following analysis (37 per cent at the 99 per cent CL). The sample of interest has been split into six different redshift bins (as done in Fig. 6) and 10 luminosity bins, from $10^8$ up to $10^{15}$ L⊙, with a logarithmic bin width equal to 0.7.

## 4.1 Method

The computation of the LF relies on the non-parametric $1/V_{max}$ method introduced by Schmidt (1968), where $V_{max}$ represents the maximum comoving volume that allows a source of a given luminosity to be observable above the flux limit of the survey. Since we consider sources coming from two different samples (GOODS-S and COSMOS), we use the method of Avni & Bahcall (1980), who extended the previous $1/V_{max}$ method to coherently join and analyse simultaneously different samples. For each source of our sample, we compute the $V_{max}$ as follows:

$$V_{max} = \int_{z_{min}}^{z_{max}} \frac{dV}{dz} \Omega(z) dz \quad (3)$$

with $z_{min}$ being the lower boundary of a given redshift bin and $z_{max}$ the minimum between the upper boundary of the redshift bin and the maximum redshift value available for the source of interest to be detected. Here $\Omega(z)$ represents the effective area of visibility and is computed as the fraction of the geometrical projected sky area $\Omega_{geom}$, where the sensitivity of the instrument is enough to let the source be observable, according to the formula:

$$\Omega(z) = \Omega_{geom} f_c(z), \quad (4)$$

where $f_c(z)$ is the flux completeness correction obtained from simulations by Berta et al. (2010, 2011). For a given luminosity and redshift bin, we computed the LF as follows:

$$\Phi(L, z) = \frac{1}{\Delta \log L} \sum_{i=1}^{n} \frac{1}{V_{\text{max},i}}, \quad (5)$$

where $\Delta \log L$ is equal to 0.7, as mentioned before. The COSMOS field has a 93 per cent redshift completeness, regardless of PACS flux densities, so that we multiply $\Phi(L, z)$ in COSMOS by a factor of 1.07, whereas the GOODS-S field is 100 per cent complete in redshift and therefore no correction is needed. Moreover, we account for the fraction of the COSMOS sample without an MIPS 24 µm counterpart inside the MIPS area. Indeed, in the context of a multicomponent SED fitting, the lack of the 24 µm detection prevents the possible torus component from being constrained, resulting in low-reliability best fits. This corresponds to ~3.5 per cent of the PEP sample, roughly independent of redshift and AGN luminosity. We take into account their contribution by uniformly increasing $\Phi(L, z)$ by a factor of 3.5 per cent.

## 4.2 AGN bolometric LF and its uncertainties

In Fig. 9, we show our accretion LF in different redshift bins. Error bars correspond to ±1σ Poissonian uncertainties (Marshall 1985), while in the case of one single object in a given luminosity bin, we provide a 90 per cent uncertainty following the recipe of Gehrels (1986). We computed separate LFs for GOODS-S and COSMOS (red and blue areas, respectively), and the total LF for the joint (GOODS-S + COSMOS) sample (black circles). As expected, in all redshift bins the GOODS-S field extends down to lower luminosities, whereas COSMOS typically samples higher luminosities.

The AGN bolometric LF shown here has been already corrected for incompleteness in accretion luminosity. We briefly describe our method as follows and refer the reader to Appendix A for a more detailed explanation. Our approach to account for the incompleteness in accretion luminosity follows that presented by Fontana et al. (2004), except for the parameters involved in such an analysis: accretion power rather than stellar mass, and 160 µm flux $S_{160}$ rather than K-band flux. Since a non-negligible fraction of active galaxies might be missed by our FIR-based selection, it is necessary to quantify such missed AGN population. Starting from an FIR-selected sample, the incompleteness in accretion luminosity might be evaluated by looking at the distribution of flux ratio between accretion flux $S_{accr}$ (defined as the bolometric flux corresponding to a given $L_{bol,AGN}$ and redshift) and $S_{160}$, traced by our AGN population. By assuming that this distribution does not change either with redshift or with Herschel flux density, it is possible to shift that down in $S_{accr}$ and to virtually sample relatively high values of $L_{bol,AGN}$ in weakly star-forming galaxies, not observable by Herschel. We iterate this shift down in $S_{accr}$ as long as the expected number of missed AGN is equal to the number of observed AGN. The latter step identifies a threshold in accretion flux $S_{accr,lim}$, corresponding to correction
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Accretion LF in different redshift bins for the PEP sample. Red and blue shaded areas trace the ±1σ Poissonian uncertainties separately for GOODS-S and COSMOS observed data points, respectively. The joint (GOODS-S and COSMOS) LF is represented by black circles, with ±1σ uncertainties (black vertical bars) related to the joint sample. The green dashed areas show the ±1σ range of ΦMC(L, z) obtained through Monte Carlo simulations (see the text for details).</figcaption></figure>

for incompleteness by a factor of 2. The curve which parametrizes the incompleteness is used to compute the maximum comoving volume where each source with $S_{accr} \ge S_{accr, lim}$, either detected or undetected by Herschel, is expected to be placed. We carried out this analysis separately for GOODS-S and COSMOS fields. As mentioned before, the AGN bolometric LF shown in Fig. 9 already incorporates the correction for incompleteness. This means that our LF is supposed to account for all active galaxies with $S_{accr} \ge S_{accr, lim}$, either above or below the Herschel detection limit. This allows us to remove the observational bias due to our FIR selection.

In the luminosity bins which are populated by both GOODS-S and COSMOS AGN at the same time, the two accretion LFs are in reasonably good agreement with each other. In the first redshift bin, we find a large gap (~1 dex) in the connection point between the two surveys, but the very low statistics (only two objects per field) makes the two LF broadly consistent within 2σ.

A Monte Carlo (MC) simulation has been performed to account for all the uncertainties in $L_{bol, AGN}$ associated with the objects in the same redshift and luminosity bin. For each object belonging to a given ($L_{bol, AGN}, z$) bin, we randomly extract an individual $L_{bol, AGN}$ value from the whole PDF distribution, using the intrinsic shape of the PDF as the weight. Then, using these random values of $L_{bol, AGN}$, we compute the LF for this simulated sample. By iterating the same MC simulation 100 times, we characterize a range of $\Phi_{MC}(L, z)$ defining the uncertainty region associated with each LF data point $\Phi(L, z)$. In Fig. 9, we show our accretion LF (black circles) with its related ±1σ confidence region (green dashed areas). The size of the uncertainty region is typically comparable to the Poissonian error bars, except for the less populated bins (large Poissonian error bars).

The limits drawn by MC simulations have been used to provide an uncertainty range for the integrated LF (see Section 5).

## 4.3 Evolution of the AGN bolometric LF

We fit our set of data points making use of a modified version of the Schechter function (Saunders et al. 1990), as follows:

$$ \Phi(L)d \log L = \Phi^* \left(\frac{L}{L^*}\right)^{1-\alpha} \exp \left[ -\frac{1}{2\sigma^2} \log_{10} \left(1 + \frac{L}{L^*}\right) \right] d \log L \quad (6) $$

which behaves as a power law for $L \ll L^*$ and as a Gaussian in $\log L$ for $L \ll L^*$. The couple of parameters ($\Phi^*, L^*$) represent the normalization and luminosity of the knee of the distribution, respectively. Two other parameters ($\alpha, \sigma$) are set to shape the low- and high-luminosity tails of the best-fitting function.

We make use of a non-linear least-squares fitting routine to find the set of parameters that better reproduce the observed $\Phi(L, z)$ through a modified Schechter function. In Fig. 10, the best-fitting functions have been overplotted on the data. The black line represents our best fit, while the red dot-dashed line shows the quasar bolometric LF as derived by Hopkins et al. (2007, hereafter H07) from X-ray data. They used a wide compilation of observed AGN LFs, from the MIR through optical, soft and hard X-ray data to build up a self-consistent AGN bolometric LF. H07 adopted a distribution of luminosity-dependent BC for each selection waveband, as well as a large baseline of intrinsic (i.e. corrected for obscuration) AGN SED shapes (from Richards et al. 2006; Steffen et al.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Accretion LF (open circles with $1\sigma$ Poissonian uncertainties) derived from PEP data. The black solid line represents our best-fitting LF whereas the red dash-dotted line is the best-fitting curve of the AGN bolometric LF, as derived by H07. In each redshift bin, we report the median redshift value of the AGN population, both for this work and for H07.</figcaption></figure>

<table><caption>**Table 1.** List of best-fitting parameters of the AGN bolometric LF with related $\pm 1\sigma$ uncertainties.</caption><thead><tr><th>z-bin</th><th>α</th><th>σ</th><th>log<sub>10</sub>(L*/L<sub>☉</sub>)</th><th>log<sub>10</sub>Φ*(Mpc<sup>-3</sup> dex<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0.1 ≤ z &lt; 0.3</td><td>1.97 ± 0.25</td><td>0.54 ± 0.31</td><td>10.91 ± 1.96</td><td>-4.35+0.93</td></tr><tr><td>0.3 ≤ z &lt; 0.7</td><td>1.48 ± 0.22</td><td>0.54 ± 0.31</td><td>11.34 ± 0.14</td><td>-4.03 ± 0.31</td></tr><tr><td>0.7 ≤ z &lt; 1.2</td><td>1.48 ± 0.22</td><td>0.54 ± 0.31</td><td>11.96 ± 0.15</td><td>-4.47 ± 0.31</td></tr><tr><td>1.2 ≤ z &lt; 1.8</td><td>1.48 ± 0.22</td><td>0.54 ± 0.31</td><td>12.61 ± 0.24</td><td>-5.00 ± 0.55</td></tr><tr><td>1.8 ≤ z &lt; 2.5</td><td>1.48 ± 0.22</td><td>0.54 ± 0.31</td><td>12.96 ± 0.18</td><td>-5.15 ± 0.44</td></tr><tr><td>2.5 ≤ z ≤ 3.8</td><td>1.48 ± 0.22</td><td>0.54 ± 0.31</td><td>12.92 ± 0.22</td><td>-5.40 ± 0.61</td></tr></tbody></table>

2006). However, the intrinsic SEDs of Richards et al. (2006) also include an additional IR ($\lambda > 1$ µm) contribution arising from the torus reprocessed radiation. This IR excess (with respect to a black-body emission in the Rayleigh-Jeans regime) has to be removed to avoid double counting the input energy of the accretion disc. As already done by Merloni & Heinz (2013), we independently evaluated such an excess to be of the order of 33 per cent. As a consequence, all the AGN bolometric luminosities taken from H07 have been reduced by a factor of ~0.33 (~0.18 dex). We split our set of data points in order to sample similar mean redshifts to H07.

There are four functional parameters controlling the Schechter function, but not more than five luminosity bins are sampled per redshift bin (i.e. there is either 0 or 1 degree of freedom). As a consequence, not all the best-fitting parameters could be sufficiently constrained at the same time. For this reason, we choose to let all the functional parameters be free in the third redshift bin only, and to fix $\alpha$ and $\sigma$ at those resulting values. In the

first redshift bin, the $\alpha$ value was allowed to vary, but $\sigma$ was fixed.

This configuration implies that the shape of the Schechter function is defined in the first and third redshift bins and is allowed to scale rigidly along the $(x,y)$ axes to fit the observed LF in the other redshift bins.

In Table 1, we list all the best-fitting parameters achieved through the least-squares fitting, together with their related $\pm 1\sigma$ errors. Since only one object populates the highest luminosity bin of the first redshift slice, the resulting best fit does not allow us to constrain either $\Phi^*$ or $L^*$ (see Table 1).

Despite the uncertainties characterizing the first z-bin, there is evidence in favour of an evolution of the parameters ($\Phi^*, L^*$) as a function of redshift. As shown in Fig. 11 (top panel), the normalization value has a maximum around $z \sim 0.5$, then it decreases to $z \sim 3$; also H07 do find a similar evolution of $\Phi^*$ with redshift. In the middle panel, the knee of the modified Schechter LF increases with redshift, at least up to $z \sim 2$. Comparing the LF of our work
<figure><img src="image_12.png" /><figcaption>Figure 11. Top panel: evolution with redshift of the normalization value $\Phi^*$, both for this work (black circles) and for H07 (red circles), taken from the respective best-fitting LF. Middle panel: evolution of the knee luminosity of the LF. Error bars correspond to the best-fitting uncertainty in each redshift bin. Bottom panel: number density of different luminosity bins as a function of redshift.</figcaption></figure>

with table 2 of H07,<fnref n="10" /> the knee of the LF is systematically at lower luminosities than what was inferred by H07 in the respective redshift bins, except at $z > 1.5$, where our $L^*$ estimates are consistent with H07. This is likely due to the fact that we are biased against early-type (and/or weakly star-forming) galaxies, which are known to host the most luminous AGN in the local Universe. Letting the number density of early-type galaxies decrease towards higher redshift, one might also expect that fewer and fewer passive galaxies host an active SMBH, as predicted by the major merger evolutionary scenario (e.g. Sanders et al. 1988). This would imply that the bulk of the SMBH growth would be progressively dominated by late-type galaxies, so that at higher redshift our $L^*$ values become consistent with those found by H07. In Fig. 11 (bottom panel), we show the evolution with redshift of the AGN number density, for different luminosity bins. As already pointed out from several works (e.g. Hasinger et al. 2005; Silverman et al. 2008) from X-ray data, we confirm an anti-hierarchical growth for SMBHs harboured in star-forming galaxies selected by Herschel, with the most luminous ones evolving faster and peaking at higher redshift.

# 5 SMBH GROWTH ACROSS COSMIC TIME

Given the estimated LF, we are able to derive the black hole accretion rate density (BHAD or $\Psi_{bhar}$) over cosmic time, from $z \sim 3$ down to the present epoch. This quantity is fundamental for characterizing the effective growth of AGN and is defined by the following expression:

$$ \Psi_{bhar}(z) = \int_{0}^{\infty} \frac{1 - \epsilon_{rad}}{\epsilon_{rad} c^2} L_{bol,AGN} \phi(L_{bol,AGN}) d \log L_{bol,AGN}, \quad (7) $$

where all the ingredients are already available. The only parameter that still needs to be adopted is the radiative efficiency. We use $\epsilon_{rad} = 0.1$ and constant with redshift and intrinsic AGN luminosity (see H07).

## 5.1 Comparison with previous results

In Fig. 12, we report our BHAD estimate compared with previous findings from the literature.

Merloni & Heinz (2008, hereafter MH08) estimated the BHAD from the evolution of the hard X-ray LF taken by Silverman et al. (2008) and accounting for a distribution of luminosity-dependent $N_H$ (hydrogen column density along the line of sight) and a set of X-ray BC from Marconi et al. (2004). By solving the continuity equation for the SMBH mass function (assuming its local value from Shankar, Weinberg & Miralda-Escudé 2009) and combining that with the evolution of the hard X-ray LF, they were able to trace the BHAD also as a function of BH mass and accretion rate. It is worth mentioning that from fig. 4 of MH08 it is shown that most of the local ($z \sim 0.1$) SMBH growth is hidden in radiatively inefficient accreting systems, whose released energy is dominated by kinetic feedback rather than radiative losses. Such ‘silent’ (Eddington ratio $\lambda_{Edd} < 3 \times 10^{-2}$) AGN population probably does not enter our sample and could be responsible for the discrepancy (by a factor of about 1.7) shown in Fig. 12 in the local BHAD between MH08, H07 and our estimate. As mentioned in Section 3.2.2, H07 built up an observed quasar LF by consistently connecting different LFs, each one computed in a single waveband.

To guarantee a coherent comparison between our study and other works taken from the literature, we checked that both MH08 and H07 have accounted for the fraction of missed AGN because of their own selection effects. MH08 estimated the fraction of highly obscured ($N_H > 1.5 \times 10^{24} \text{cm}^{-2}$) AGN from the X-ray background synthesis models by Gilli, Comastri & Hasinger (2007). In particular, they calculated the expected fraction between observed and unabsorbed 2–10 keV X-ray luminosity, by assuming a canonical AGN X-ray spectrum (absorbed power law with photon index $\Gamma = 1.8$; Tozzi et al. 2006), and as a function of obscuration in the range $21 < \log(N_H/\text{cm}^{-2}) < 24$. For $\log(N_H/\text{cm}^{-2}) > 24$, they assumed that the ratio between observed and intrinsic X-ray luminosity was of the order of 2 per cent (Gilli et al. 2007).

H07 estimated the incompleteness fraction with a set of luminosity-dependent $N_H$ distributions. First, they collected a wide set of luminosity-dependent AGN spectral shapes and calculated the BC distributions, in the case of no obscuration. In addition, they implemented three possible models for obscuration. For each parametrization, H07 calculated the expected amount of extinction in X-rays, optical and MIR. Finally, through a convolution between AGN templates, BC distributions and $N_H$ distributions, they were able to predict the incompleteness fraction of AGN in each band.

<footnote n="10">¹⁰ Note that the $L^*$ values derived by H07 have been corrected in this work by a factor of ~0.18 dex.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Black hole accretion rate density estimate from the AGN bolometric LF as a function of redshift (black circles). The red shaded area shows the ±1σ uncertainty region. Previous estimates from different selection wavelengths (from Merloni & Heinz 2008, and H07) are reported for comparison.</figcaption></figure>

Our estimates (black circles) have been derived by integrating the best-fitting curve of the LF, down to $10^8$ L⊙, in each redshift bin. The red shaded area represents the ±1σ uncertainty region. The latter has been computed by accounting for the 100 different integrated $\Phi_{MC}(L, z)$ values and plotting their cumulative distribution at 16 and 84 percentiles in each redshift bin. The BHAD seems to evolve quickly from $z > 3$ to $z \sim 2$, where it shows a peak, and then decreases towards the present epoch. Despite our AGN sample having been selected and analysed independently of the previous ones, the overall trend is consistent with others taken from the literature, as shown in Fig. 5. This is not surprising, as both our estimate and that obtained by H07 and MH08 have been already corrected for their own incompleteness effects and should be in fair agreement as a proof of their mutual consistency. All the estimates are plotted for the canonical case $\epsilon_{rad} = 0.1$. A different value of $\epsilon_{rad}$ simply results in a change of the BHAD normalization.

An important constraint that allows us to test the consistency of the integrated BHAD($z$) is to deal with Soltan's argument, i.e. deriving the local BH mass density and comparing it to the estimate achieved from observations. This represents the overall energy density released by SMBHs during cosmic time, down to $z = 0$. One of the most recent estimates comes from Shankar et al. (2009), as mentioned before. Tuning the integrated BHAD evolution in order to reproduce the observed $\rho_{bh,0}$ is a way to check the range spanned by the radiative efficiency $\epsilon_{rad}$. From the standard relativistic accretion theory, it is known that $\epsilon_{rad} = 0.06-0.20$ (Novikov & Thorne 1973). We integrate the BHAD over the redshift range from $z = 6$ to $10^{-5}$ (but taking $z = \infty$ as the upper boundary would change our result by a factor of $< 1$ per cent). By assuming $\epsilon_{rad} = 0.1$, as shown in Fig. 13, our local BH mass density is

$$ \rho_{bh,0} = 3.1^{+1.0}_{-0.8} \times 10^5 \, \text{M}_\odot \, \text{Mpc}^{-3} \quad (8) $$

which is consistent with previous estimates. Indeed, for a typical matter-to-radiation conversion efficiency $\epsilon_{rad} = 0.1$, Yu & Tremaine (2002) found $\rho_{bh,0} = 2.9 \pm 0.5 \times 10^5 \, \text{Mpc}^{-3}$; Marconi et al.

(2004) reported $\rho_{bh,0} = (4.6^{+1.9}_{-0.9}) \times 10^5 \, \text{Mpc}^{-3}$; Shankar et al. (2009) achieved $\rho_{bh,0} = (4.2^{+1.0}_{-0.0}) \times 10^5 \, \text{Mpc}^{-3}$. All these measurements come from integration of the SMBH mass function by assuming $\epsilon_{rad} = 0.1$.

<figure><img src="image_14.png" /><figcaption>**Figure 13.** SMBH radiative efficiency $\epsilon_{rad}$ as a function of the local BH mass density predicted by this work (black solid line, with ±1σ uncertainty shown by the dashed lines). The vertical red dashed line shows the estimate of the observed local BH mass density by Shankar et al. (2009).</figcaption></figure>

To fully reconcile the SMBH growth derived here with the observed local BH mass density (e.g. the one from Shankar et al. 2009), we need to adopt $\epsilon_{rad} = 0.076^{+0.023}_{-0.018}$, which is in broad agreement with previous measurements. The typical 10 per cent efficiency comes from Soltan (1982), and corroborated by Fabian & Iwasawa (1999). Later on, Merloni, Rudnick & Di Matteo (2004) found $0.04 < \epsilon_{rad} < 0.12$ and H07 report $\rho_{bh,0} = 4.81^{+1.24}_{-0.99} \times 10^5 \, \text{Mpc}^{-3}$, which results in a slightly higher radiative efficiency ($\epsilon_{rad} \sim 0.1$). Finally, MH08 end up with an intermediate value ($\sim 0.07$), which is in good agreement with our results.

# 6 CONCLUSIONS

We have presented an estimate of the cosmic evolution of SMBH growth, from $z \sim 3$ down to the present epoch. Our estimate of the BHAD across the cosmic time is the first one derived by using an
FIR-selected sample. We selected 4343 PEP sources with 160 µm flux density >3σ detection limit in the GOODS-S and COSMOS fields. The available multiband photometric coverage allowed us to explore the observed SEDs from the UV to the FIR, in particular thanks to the valuable PACS and SPIRE data from the Herschel satellite. Moreover, the availability of Herschel data both in the GOODS-S and in COSMOS fields was essential to explore different and complementary luminosity regimes.

Our analysis relies on a robust broad-band SED decomposition performed to decouple the AGN contribution from the host galaxy content in each global SED. We determined reliable estimates of the AGN accretion luminosity, without any bias against highly obscured AGN and based on smaller BC than those adopted from X-rays (~4 instead of 20–30).

Our main conclusions are as follows.

(i) The percentage of star-forming galaxies in the redshift range $0.1 < z < 3.8$ harbouring an AGN (with a 99 per cent significance) is as high as 37 per cent.

(ii) At $z > 0.3$ we find the AGN detection rate rising as a function of the total (1–1000 µm rest-frame) IR luminosity, suggesting that the AGN activity at these redshifts is more likely to take place in enhanced ($L_{1-1000} \sim 10^{12} L_{\odot}$) star-forming systems, rather than in ‘normal’ ($L_{1-1000} \sim 10^{11} L_{\odot}$) ones. Furthermore, we find also an additional dependence of the AGN detection rate on the redshift, since the probability of picking up AGN in galaxies with a given luminosity increases with redshift, at least up to $z \sim 2$. At $z < 0.3$, the average AGN detection rate reaches 5–10 per cent and does not evolve significantly with IR luminosity.

(iii) We compare our classification based on SED decomposition with that of the Lacy et al. (2007) and Donley et al. (2012) colour–colour AGN selection criteria. As expected, the latter is more reliable and conservative than the former: ~17 per cent of PEP AGN fulfil the Donley et al. (2012) criterion and are characterized by AGN-dominated systems in the MIR, with extreme ($>10^{12} L_{\odot}$) star formation and AGN accretion luminosities. On the other hand, ~73 per cent of PEP AGN lie inside the Lacy et al. (2007) wedge: these sources share intermediate properties between QSO-like SEDs and purely star-forming systems, showing more modest IR luminosities $L_{1-1000} \sim 10^{11-11.5} L_{\odot}$. In addition, we point out that our PEP AGN and PEP galaxies occupy different regions across the MIR colour–colour diagram. Such differences highlight the fact that our SED classification effectively labels as ‘AGN’ sources with near/mid-IR SEDs systematically different from those of purely SF galaxies.

(iv) The best-fitting AGN bolometric LF is in good agreement with that obtained by H07, despite the AGN selection criteria being completely independent one from the other.

(v) We find clear evidence in favour of ‘downsizing’ behaviour of PEP AGN. Indeed, both from the evolution of the best-fitting LF (Fig. 11, top and middle panels) and from the evolution of the number density with redshift (Fig. 11, bottom panel), we see that the bulk of the SMBH growth takes place in more powerful objects at earlier times. The similarity between our findings and those previously suggested by X-ray-based works (e.g. Hasinger et al. 2005; Silverman et al. 2008) may imply that AGN hosts, both from X-rays and from FIR selection, have likely experienced similar evolutionary paths.

(vi) The evolution of the BHAD and the resulting BH radiative efficiency ($\epsilon_{\text{rad}} \sim 0.07$) presented in this paper are in agreement with those obtained from X-ray-based results. This overall consistency suggests that the FIR selection, together with a proper SED

decomposition, allows us to obtain similar integrated luminosity densities with respect to X-ray-selected AGN samples. However, this does not necessarily mean that the two selection methods trace the evolution of the same AGN population.

In this work, we have aimed at constraining the mean evolutionary trend of the SMBH growth as seen by Herschel-selected galaxies. However, our analysis is not deep enough to shed light on the nature of our AGN sample, or on the evolution of the Compton-thick AGN population. This issue would require principally cross-matching with available X-ray data. In a future work, we will investigate more thoroughly the physical properties of the AGN populations by combining IR and X-ray observations, so as to obtain a more complete perspective on the AGN evolution through cosmic time.

## ACKNOWLEDGEMENTS

This paper uses data from Herschel's photometers PACS and SPIRE. PACS has been developed by a consortium of institutes led by MPE (Germany) and including: UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy) and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

ID is deeply grateful to Iary Davidzon, Andrea Negri, Fabio Vito and Alessandro Marconi for inspiring discussions and suggestions. ID also thanks Elisabeta Lusso for useful tests on codes for SED decomposition and Micol Bolzonella for kindly providing a tool to estimate error bars.

## REFERENCES

* Aird J. et al., 2010, MNRAS, 401, 2531

* Akylas A., Georgantopoulos I., Georgakakis A., Kitsionas S., Hatziminaoglou E., 2006, A&A, 459, 693

* Alexander D. M., Bauer F. E., Chapman S. C., Smail I., Blain A. W., Brandt W. N., Ivison R. J., 2005, ApJ, 632, 736

* Alonso-Herrero A., Pereira-Santaella M., Rieke G. H., Rigopoulou D., 2012, ApJ, 744, 2

* Avni Y., Bahcall J. N., 1980, ApJ, 235, 694

* Balestra I. et al., 2010, A&A, 512, A12

* Ball W. H., Tout C. A., Żytkow A. N., Eldridge J. J., 2011, MNRAS, 414, 2751

* Ballantyne D. R., Everett J. E., Murray N., 2006, ApJ, 639, 740

* Barger A. J., Cowie L. L., 2005, ApJ, 635, 115

* Barmby P. et al., 2006, ApJ, 642, 126

* Bell E. F. et al., 2005, ApJ, 625, 23

* Benson A. J., Bower R. G., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2003, ApJ, 599, 38

* Berta S. et al., 2010, A&A, 518, L30

* Berta S. et al., 2011, A&A, 532, A49

* Berta S. et al., 2013, A&A, 551, A100

* Bevington P. R., Robinson D. K., 2003, Data Reduction and Error Analysis for the Physical Sciences. McGraw-Hill, New York
Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Brandt W. N., Hasinger G., 2005, ARA&A, 43, 827

Brusa M. et al., 2010, ApJ, 716, 348

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Bundy K. et al., 2006, ApJ, 651, 120

Bunker A. J., Stanway E. R., Ellis R. S., McMahon R. G., McCarthy P. J., 2003, MNRAS, 342, L47

Chabrier G., 2003, ApJ, 586, L133

Charlot S., Fall S. M., 2000, ApJ, 539, 718

Ciliegi P., Gruppioni C., McMahon R., Rowan-Robinson M., 2001, Ap&SS, 276, 957

Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038

Ciotti L., Ostriker J. P., Proga D., 2010, ApJ, 717, 708

Cooper M. C. et al., 2012, MNRAS, 425, 2116

Cowie L. L., Songaila A., Hu E. M., Cohen J. G., 1996, AJ, 112, 839

Cowie L. L., Hu E. M., Songaila A., Egami E., 1997, ApJ, 481, L9

Cristiani S. et al., 2000, A&A, 359, 489

Croom S. M., Smith R. J., Boyle B. J., Shanks T., Loaring N. S., Miller L., Lewis I. J., 2001, MNRAS, 322, L29

Croton D. J. et al., 2006, MNRAS, 365, 11

da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595

Daddi E. et al., 2007, ApJ, 670, 156

Della Ceca R. et al., 2008, A&A, 487, 119

Devecchi B., Volonteri M., 2009, ApJ, 694, 302

Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604

Dickinson M. et al., 2004, ApJ, 600, L99

Doherty M., Bunker A. J., Ellis R. S., McCarthy P. J., 2005, MNRAS, 361, 525

Donley J. L. et al., 2012, ApJ, 748, 142

Dunne L. et al., 2009, MNRAS, 394, 3

Elbaz D. et al., 2007, A&A, 468, 33

Elbaz D. et al., 2011, A&A, 533, A119

Elvis M., Maccacaro T., Wilson A. S., Ward M. J., Penston M. V., Fosbury R. A. E., Perola G. C., 1978, MNRAS, 183, 129

Fabian A. C., Iwasawa K., 1999, MNRAS, 303, L34

Fadda D. et al., 2010, ApJ, 719, 425

Fanidakis N. et al., 2012, MNRAS, 419, 2797

Feltre A., Hatziminaoglou E., Fritz J., Franceschini A., 2012, MNRAS, 426, 120

Ferrarese L., 2002, ApJ, 578, 90

Ferrarese L., Merritt D., 2000, ApJ, 539, L9

Fontana A. et al., 2004, A&A, 424, 23

Fritz J., Franceschini A., Hatziminaoglou E., 2006, MNRAS, 366, 767

Gabor J. M., Davé R., Oppenheimer B. D., Finlator K., 2011, MNRAS, 417, 2676

Gebhardt K. et al., 2000, ApJ, 539, L13

Gehrels N., 1986, ApJ, 303, 336

Gerke B. F. et al., 2007, MNRAS, 376, 1425

Gilli R., Comastri A., Hasinger G., 2007, A&A, 463, 79

Gilli R., Vignali C., Mignoli M., Iwasawa K., Comastri A., Zamorani G., 2010, A&A, 519, A92

Grazian A. et al., 2006, A&A, 449, 951

Griffin M. J. et al., 2010, A&A, 518, L3

Gruppioni C. et al., 2013, MNRAS, 432, 23

Gültekin K., Cackett E. M., Miller J. M., Di Matteo T., Markoff S., Richstone D. O., 2009, ApJ, 706, 404

Hasinger G., 2008, A&A, 490, 905

Hasinger G., Miyaji T., Schmidt M., 2005, A&A, 441, 417

Hatziminaoglou E. et al., 2008, MNRAS, 386, 1252

Hatziminaoglou E. et al., 2010, A&A, 518, L33

Hopkins P. F., Quataert E., 2010, MNRAS, 407, 1529

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Robertson B., Springel V., 2006, ApJS, 163, 1

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731 (H07)

Hubeny I., Agol E., Blaes O., Krolik J. H., 2000, ApJ, 533, 710

Hwang H. S. et al., 2010, MNRAS, 409, 75

Ilbert O. et al., 2010, ApJ, 709, 644

Imanishi M., Nakagawa T., Shirahata M., Ohyama Y., Onaka T., 2010, ApJ, 721, 1233

Iwasawa K. et al., 2012, A&A, 546, A84

Jaffe W. et al., 2004, Nature, 429, 47

Juneau S. et al., 2005, ApJ, 619, L135

Juneau S., Dickinson M., Alexander D. M., Salim S., 2011, ApJ, 736, 104

Juneau S. et al., 2013, ApJ, 764, 176

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kormendy J., Kennicutt J. R. C., 2004, ARA&A, 42, 603

Kurk J. et al., 2013, A&A, 549, A63

La Franca F. et al., 2005, ApJ, 635, 864

Lacy M. et al., 2004, ApJS, 154, 166

Lacy M., Petric A. O., Sajina A., Canalizo G., Storrie-Lombardi L. J., Armus L., Fadda D., Marleau F. R., 2007, AJ, 133, 186

Lacy M. et al., 2013, preprint (arXiv:1308.4190)

Laird E. S., Nandra K., Pope A., Scott D., 2010, MNRAS, 401, 2763

Lawrence A., 1991, MNRAS, 252, 586

Le Fèvre O. et al., 2005, A&A, 439, 845

Le Floc'h E. et al., 2005, ApJ, 632, 169

Le Floc'h E. et al., 2009, ApJ, 703, 222

Lilly S. J., Le Fèvre O., Hammer F., Crampton D., 1996, ApJ, 460, L1

Lilly S. J. et al., 2009, ApJS, 184, 218

Lodato G., Natarajan P., 2007, MNRAS, 377, L64

Lusso E. et al., 2012, MNRAS, 425, 623

Lusso E. et al., 2013, ApJ, 777, 86

Lutz D., Spoon H. W. W., Rigopoulou D., Moorwood A. F. M., Genzel R., 1998, ApJ, 505, L103

Lutz D. et al., 2010, ApJ, 712, 1287

Lutz D. et al., 2011, A&A, 532, A90

Lynden-Bell D., 1969, Nature, 223, 690

Madau P., Ferguson H. C., Dickinson M. E., Giavalisco M., Steidel C. C., Fruchter A., 1996, MNRAS, 283, 1388

Magdis G. E., Rigopoulou D., Huang J.-S., Fazio G. G., 2010, MNRAS, 401, 1521

Magdis G. E. et al., 2013, A&A, 558, A136

Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D., Frayer D. T., Willmer C. N. A., 2009, A&A, 496, 57

Magorrian J. et al., 1998, AJ, 115, 2285

Marconi A., Hunt L. K., 2003, ApJ, 589, L21

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169

Marshall H. L., 1985, ApJ, 299, 109

Menci N., Fontana A., Giallongo E., Grazian A., Salimbeni S., 2006, ApJ, 647, 753

Menci N., Fiore F., Puccetti S., Cavaliere A., 2008, ApJ, 686, 219

Merloni A., Heinz S., 2008, MNRAS, 388, 1011 (MH08)

Merloni A., Heinz S., 2013, in Oswalt T. D., Keel W. C., eds, Evolution of Active Galactic Nuclei. Springer, Dordrecht, p. 503

Merloni A., Rudnick G., Di Matteo T., 2004, MNRAS, 354, L37

Mignoli M. et al., 2005, A&A, 437, 883

Mullaney J. R., Alexander D. M., Goulding A. D., Hickox R. C., 2011, MNRAS, 414, 1082

Mullaney J. R. et al., 2012, ApJ, 753, L30

Nardini I., Risaliti G., Watabe Y., Salvati M., Sani E., 2010, MNRAS, 405, 2505

Natarajan P., 2011, preprint (arXiv:1105.4902)

Netzer H., 2009, MNRAS, 399, 1907

Novikov I. D., Thorne K. S., 1973, in Dewitt C., Dewitt B. S., eds, Black Holes (Les Astres Occlus). Gordon and Breach, New York, p. 343

Oliver S. J. et al., 2012, MNRAS, 424, 1614

Olsen K. P., Rasmussen J., Toft S., Zirm A. W., 2013, ApJ, 764, 4

Peng Y.-j. et al., 2010, ApJ, 721, 193

Peng Y.-j., Lilly S. J., Renzini A., Carollo M., 2012, ApJ, 757, 4

Pérez-González P. G. et al., 2005, ApJ, 630, 82

Pier E. A., Krolik J. H., 1992, ApJ, 401, 99

Pilbratt G. L. et al., 2010, A&A, 518, L1

Poglitsch A. et al., 2010, A&A, 518, L2

Polletta M. et al., 2007, ApJ, 663, 81
Pope A. et al., 2008, ApJ, 675, 1171
Popesso P. et al., 2009, A&A, 494, 443
Pozzi F. et al., 2007, A&A, 468, 603
Pozzi F. et al., 2012, MNRAS, 423, 1909
Rees M. J., 1978, in Berkhuijsen E. M., Wielebinski R., eds, Proc. IAU Symp. 77, Structure and Properties of Nearby Galaxies. Reidel, Dordrecht, p. 237
Rees M. J., 1984, ARA&A, 22, 471
Richards G. T. et al., 2006, ApJS, 166, 470
Risaliti G., Elvis M., 2004, in Barger A. J., ed., Astrophysics and Space Science Library, Vol. 308, Supermassive Black Holes in the Distant Universe. Kluwer, Dordrecht, p. 187
Rodighiero G. et al., 2010, A&A, 518, L25
Rodighiero G. et al., 2011, ApJ, 739, L40
Rosario D. J. et al., 2012, A&A, 545, A45
Roseboom I. G. et al., 2012, MNRAS, 419, 2758
Rowan-Robinson M. et al., 2010, MNRAS, 409, 2
Sajina A., Spoon H., Yan L., Imanishi M., Fadda D., Elitzur M., 2009, ApJ, 703, 270
Sajina A., Yan L., Fadda D., Dasyra K., Huynh M., 2012, ApJ, 757, 13
Salpeter E. E., 1964, ApJ, 140, 796
Sanders D. B., Soifer B. T., Elias J. H., Madore B. F., Matthews K., Neugebauer G., Scoville N. Z., 1988, ApJ, 325, 74
Santini P. et al., 2009, A&A, 504, 751
Santini P. et al., 2012, A&A, 540, A109
Saunders W., Rowan-Robinson M., Lawrence A., Efstathiou G., Kaiser N., Ellis R. S., Frenk C. S., 1990, MNRAS, 242, 318
Sazonov S., Revnivtsev M., Krivonos R., Churazov E., Sunyaev R., 2007, A&A, 462, 57
Schawinski K. et al., 2006, Nature, 442, 888
Schawinski K., Simmons B. D., Urry C. M., Treister E., Glikman E., 2012, MNRAS, 425, L61
Schmidt M., 1968, ApJ, 151, 393
Scoville N. et al., 2007, ApJS, 172, 1
Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
Shankar F., Salucci P., Granato G. L., De Zotti G., Danese L., 2004, MNRAS, 354, 1020
Shankar F., Weinberg D. H., Miralda-Escudé J., 2009, ApJ, 690, 20
Shao L. et al., 2010, A&A, 518, L26
Silverman J. D. et al., 2008, ApJ, 679, 118
Simpson C., 2005, MNRAS, 360, 565
Smith D. J. B. et al., 2012, MNRAS, 427, 703
Soltan A., 1982, MNRAS, 200, 115
Springel V. et al., 2005, Nature, 435, 629
Stanway E. R. et al., 2004a, ApJ, 604, L13
Stanway E. R., Bunker A. J., McMahon R. G., Ellis R. S., Treu T., McCarthy P. J., 2004b, ApJ, 607, 704
Steffen A. T., Barger A. J., Cowie L. L., Mushotzky R. F., Yang Y., 2003, ApJ, 596, L23
Steffen A. T., Strateva I., Brandt W. N., Alexander D. M., Koekemoer A. M., Lehmer B. D., Schneider D. P., Vignali C., 2006, AJ, 131, 2826
Stern D. et al., 2005, ApJ, 631, 163
Strolger L.-G. et al., 2004, ApJ, 613, 200
Sutherland W., Saunders W., 1992, MNRAS, 259, 413
Symeonidis M. et al., 2013, MNRAS, 431, 2317
Szokoly G. P. et al., 2004, ApJS, 155, 271
Teplitz H. I. et al., 2011, AJ, 141, 1
Tozzi P. et al., 2006, A&A, 451, 457
Treister E., Urry C. M., 2006, ApJ, 652, L79
Treister E., Urry C. M., Virani S., 2009, ApJ, 696, 110
Treister E., Schawinski K., Urry C. M., Simmons B. D., 2012, ApJ, 758, L39
Tremaine S. et al., 2002, ApJ, 574, 740
Tristram K. R. W. et al., 2007, A&A, 474, 837
Tristram K. R. W. et al., 2009, A&A, 502, 67
Trump J. R. et al., 2009, ApJ, 696, 1195
Trump J. R. et al., 2013, ApJ, 763, L6
Ueda Y., Akiyama M., Ohta K., Miyaji T., 2003, ApJ, 598, 886

van der Wel A., Franx M., van Dokkum P. G., Rix H.-W., 2004, ApJ, 601, L5
Vanzella E. et al., 2008, A&A, 478, 83
Vasudevan R. V., Fabian A. C., 2009, MNRAS, 392, 1124
Vito F. et al., 2013, MNRAS, 428, 354
Volonteri M., Begelman M. C., 2010, MNRAS, 409, 1022
Wall J. V., Pope A., Scott D., 2008, MNRAS, 383, 435
Xue Y. Q. et al., 2011, ApJS, 195, 10
Yu Q., Tremaine S., 2002, MNRAS, 335, 965
Yuan T.-T., Kewley L. J., Sanders D. B., 2010, ApJ, 709, 884
Zheng W., Kriss G. A., Telfer R. C., Grimes J. P., Davidsen A. F., 1997, ApJ, 475, 469

## APPENDIX A: INCOMPLETENESS IN ACCRETION LUMINOSITY

As mentioned in Section 4.2, in this appendix we describe in detail our approach to account for the incompleteness in accretion luminosity. We follow a similar approach to that developed by Fontana et al. (2004), who build the stellar mass function correcting for the incompleteness in stellar mass, starting from a K-band-selected sample. The same argument could be in principle generalized, under some assumptions, to deal with the evaluation of the incompleteness affecting some unobservable quantity when selecting sources with any other observable quantity. We make use of a similar algorithm, but referred to the 'accretion flux' $S_{accr}$ as unobservable and to the 160 µm flux as observable. We define $S_{accr} = L_{bol,AGN}/(4\pi D_L^2)$ and compare it with the corresponding 160 µm flux, simply rescaling the accretion power by its luminosity distance factor $D_L$.

As described in Section 2, the reference sample was selected above a given detection threshold (2.4 mJy in the GOODS-S, 10.2 mJy in the COSMOS field). Nevertheless, we are potentially ruling out some objects that are below the observational limit, but likely having a relatively high accretion power (i.e. a larger $L_{bol,AGN}$ value with respect to sources with FIR flux density above the detection limit). In order to quantify the missing fraction of AGN as a function of the observed 160 µm flux density, the following hypotheses have been made.

(i) We assume that the distribution of the flux ratio $S_{accr}/S_{160}$ does not depend on redshift. Indeed, we are not able to characterize such distribution over sizeable samples at any redshift, especially at $z > 2$ where the low statistics prevents us from deriving the ratio with high significance.

(ii) We rely on the basic assumption that the $S_{160}$ completeness function derived by Berta et al. (2010, 2011) for 160 µm sources in GOODS-S preserves the same trend also below the detection limit. Furthermore, the $S_{160}$ flux completeness function derived for COSMOS sources has been extrapolated below 10.2 mJy on the existing one, by assuming the same trend for that derived in GOODS-S.

We fit the logarithmic distribution of flux ratio with a Gaussian function, for GOODS-S and COSMOS samples separately. In support to the assumption (ii), we observed a similar distribution of flux ratio among COSMOS and GOODS-S, which means that most likely the two fluxes scale in a similar way as a function of $S_{160}$. Such a similar trend at lower $S_{160}$ allows us to assume that our distribution is also preserved for lower $S_{160}$ than the GOODS-S detection limit. As expected, the distributions of flux ratio are quite broad (about 0.6 dex), so they do not imply a tight correlation between star formation and AGN accretion activities (see e.g. Rosario et al. 2012). However, the approach developed by Fontana et al. (2004) allows us to take into account the full distribution of flux ratios, even if quite broad.
<figure><img src="image_15.png" /><figcaption>Figure A1. Distribution of GOODS-S (red asterisks) and COSMOS (blue squares) sources classified as AGN, as a function of S<sub>accr</sub> and S<sub>160</sub>. Dotted lines show the flux density limits at 160 µm (S<sub>160, lim</sub>), whereas horizontal dashed lines set the thresholds in accretion flux that we found through this analysis. The red and blue solid lines trace the upper value of the distribution of flux ratio S<sub>accr</sub>/S<sub>160</sub>, as a function of the accretion flux. The red and blue shaded areas define the flux-to-flux regions that we have explored to calibrate our correction for incompleteness in S<sub>accr</sub>. See the text for details.</figcaption></figure>

In order to correct for the incompleteness in S<sub>accr</sub>, it is necessary to evaluate the intrinsic number of AGN that are missed because S<sub>160</sub> is lower than the detection limit S<sub>160, lim</sub>, as a function of the accretion flux. In Fig. A1, we plot the distribution of COSMOS (blue squares) and GOODS-S (red asterisks) sources as a function of both 160 µm and accretion fluxes. The vertical dotted lines set the FIR detection limit of each PEP field. Each coloured solid line marks the upper boundary of the distribution of flux ratios sampled by each value of the accretion flux. For sake of simplicity, hereafter we focus on the red asterisks only (representing GOODS-S). The following reasoning may be also extended to the blue squares (representing COSMOS).

Following Fontana et al. (2004), we set the starting point for our analysis in such a way that the highest flux ratio of the distribution samples the flux density limit S<sub>160, lim</sub> value of the survey (see the downward arrow placed where the solid line crosses the dotted one in Fig. A1). We may reasonably assume that above such a threshold we are 100 per cent complete in accretion flux.

Since we are interested in quantifying the missing AGN fraction as a function of the accretion flux, we follow the downward arrow (see Fig. A1), moving along the solid line down to lower and lower values of S<sub>accr</sub>. According to the hypothesis (ii), we assume that also the distribution shifts down with the same trend. For each step in S<sub>accr</sub>, we need to evaluate the expected (i.e. intrinsic) number of sources N<sub>exp</sub>, which is nothing more than the ‘classic’ log N<sub>exp</sub> − log S<sub>160</sub>. This is achievable by extrapolating the existing S<sub>160</sub> completeness curve down to lower and lower FIR flux densities.

As we are interested in deriving such information as a function of the accretion flux rather than FIR flux density, we convert the

obtained log N<sub>exp</sub> − log S<sub>160</sub> into log N<sub>exp</sub> − log S<sub>accr</sub> through a convolution with the distribution of flux ratios. In other words, we use our distribution to derive the expected number of sources for each input accretion flux. From the comparison between the observed and expected number of sources, it is possible to obtain the detectable AGN fraction f<sub>det</sub> as a function of S<sub>accr</sub>:

$$f_{\text{det}}(S_{\text{accr}}) = \frac{\int_{S_{160,\text{lim}}}^{S_{160,\text{max}}} \frac{dN_{\text{exp}}}{dS_{160}} g(\bar{x}) dS_{160}}{\frac{S_{160,\text{max}}}{S_{160,\text{min}}} \frac{dN_{\text{exp}}}{dS_{160}} g(\bar{x}) dS_{160}} \quad (\text{A1})$$

where (x̄) = S<sub>accr</sub>/S<sub>160</sub> and g(x̄) is the interpolated value of the distribution of flux ratios at [S<sub>accr</sub>, S<sub>160</sub>]. While S<sub>160, lim</sub> is constant, S<sub>160, min</sub> and S<sub>160, max</sub> change step by step in accretion flux, since they represent the minimum and maximum FIR fluxes, respectively, sampled by the distribution of flux ratios for a given accretion flux S<sub>accr</sub>.

The quantity in equation (A1) is by definition ≤1 and scales downwards as long as S<sub>accr</sub> decreases, since the FIR-based selection becomes progressively more incomplete in accretion flux. As done by Fontana et al. (2004), we choose to run the iterations until the detectable fraction matches the lost fraction (i.e. f<sub>det</sub> = 0.5). The latter value identifies the minimum accretion flux S<sub>accr, lim</sub> that a source should have to enter a >50 per cent complete accretion-based sample, regardless of its FIR flux density. In Fig. A1, the shaded areas delimit the flux-to-flux regions that we have explored to calibrate our correction for incompleteness in accretion flux. Such thresholds in accretion flux are equal to 4 × 10<sup>-14</sup> and 5 × 10<sup>-13</sup> erg s<sup>-1</sup> cm<sup>-2</sup> in the GOODS-S and the COSMOS field, respectively (horizontal dashed lines in Fig. A1).

Within our AGN sample, only those sources with accretion fluxes larger than the corresponding threshold have been considered for the AGN bolometric LF. By applying the above mentioned corrections, the predicted number of sources in each L<sub>bol, AGN</sub> bin may increase up to a factor of 2. Nevertheless, the V<sub>max</sub> method (equation 3) was referred to the accessible comoving volume for a source to be detected at 160 µm, but not necessarily to enter at the same time an accretion-based selection. As a consequence, we define an ‘effective’ maximum comoving volume V<sub>max</sub><sup>eff</sup> as

$$V_{\text{max}}^{\text{eff}} = \min[V_{\text{max}}, V_{\text{max,accr}}], \quad (\text{A2})$$

where V<sub>max, accr</sub> is the equivalent of the classic V<sub>max</sub> definition, except for the fact that V<sub>max, accr</sub> is computed as a function of the accretion flux rather than of FIR flux. In other words, we calculate a V<sub>max, accr</sub> value for each object by following equations (3) and (4) and replacing the FIR flux completeness function f<sub>c</sub>(z) with that based on the accretion flux f<sub>det</sub>(S<sub>accr</sub>).

Finally, the effective volume V<sub>max</sub><sup>eff</sup> (equation A2) is computed for sources having S<sub>accr</sub> ≥ S<sub>accr, lim</sub> and S<sub>160</sub> ≥ S<sub>160, lim</sub>, fitting both simultaneously with FIR and accretion-based selections. At the same time, V<sub>max</sub><sup>eff</sup> allows us to correct for the incompleteness in accretion flux, that is to place our Herschel-selected AGN within the appropriate comoving volume where all AGN with S<sub>accr</sub> ≥ S<sub>accr, lim</sub>, either detected or undetected by Herschel, would be observable.

This paper has been typeset from a TeX/LA TeX file prepared by the author.