# The Spitzer mid-infrared AGN survey. II-the demographics and cosmic evolution of the AGN population

## Abstract

We present luminosity functions derived from a spectroscopic survey of AGN selected from Spitzer Space Telescope imaging surveys. Selection in the mid-infrared is significantly less affected by dust obscuration. We can thus compare the luminosity functions of the obscured and unobscured AGN in a more reliable fashion than by using optical or X-ray data alone. We find that the AGN luminosity function can be well described by a broken power-law model in which the break luminosity decreases with redshift. At high redshifts (), we find significantly more AGN at a given bolometric luminosity than found by either optical quasar surveys or hard X-ray surveys. The fraction of obscured AGN decreases rapidly with increasing AGN luminosity, but, at least at high redshifts, appears to remain at % even at bolometric luminosities . The data support a picture in which the obscured and unobscured populations evolve differently, with some evidence that high luminosity obscured quasars peak in space density at a higher redshift than their unobscured counterparts. The amount of accretion energy in the Universe estimated from this work suggests that AGN contribute about 12% to the total radiation intensity of the Universe, and a high radiative accretion efficiency is required to match current estimates of the local mass density in black holes.

quasars:general – galaxies:Seyfert – infrared:galaxies – galaxies:starburst
8

## 1. Introduction

The existence of a large population of obscured AGN has been recognized for as long as AGN have been known. The vast majority of well-studied obscured objects, however, are local, of relatively low bolometric luminosity. The fraction of obscured AGN at the high (quasar) end of the AGN luminosity function, where the population is at cosmological distances has been the subject of debate. Reliable statistics exist only for the % of radio-loud AGN, where luminous radio galaxies are seen to exist in similar numbers to radio-loud quasars when selected on an isotropic property such as low frequency radio emission (e.g. Willott et al. 2000). The situation for radio-quiet objects is less certain, and it is only in the last decade that progress has been made in quantifying the population of dust obscured radio-quiet quasars through advances in X-ray, optical and mid-infrared selection techniques.

Selection of AGN in the mid-infrared allows the discovery of powerful AGN and quasars whose optical and soft X-ray emission is hidden by dust (e.g. Lacy et al. 2004, 2007; Stern et al. 2005; 2012; Martínez-Sansigre et al. 2005; Donley et al. 2012; Eisenhardt et al. 2012). Mid-infrared selection has the powerful characteristic that it permits the selection of samples of AGN containing both moderately obscured and unobscured objects of similar bolometric luminosities, allowing an estimate to be made of the importance of the obscured AGN population to the AGN population as a whole. Mid-infrared selection is also a useful complement to other techniques for finding obscured quasars such as hard X-rays (Norman et al. 2002; Brusa et al. 2010), which tend to be restricted to smaller fields and thus lower-luminosity objects, and optical techniques based on narrow-line emission (Zakamska et al. 2003; Alexandroff et al. 2013), which are restricted to redshift ranges in which strong emission lines fall.

Mid-infared selection has its limitations, for example, Juneau et al. (2013) show that, at and low luminosities (), both X-ray and mid-infrared selection miss a large fraction of AGN that are recovered by the Mass-Excitation technique, an optical diagnostic which compares the [Oiii]/H ratio to the stellar mass (Juneau et al. 2011). These objects are probably of such low luminosity that hot dust emission does not stand out against stellar light in mid-infrared color-color plots. Similar incompletenesses to lower luminosity AGN have been noted by Hickox et al. (2009), Park et al. (2010), Mendez et al. (2013) and Messias et al. (2014). In particular, Messias et al. show that the completeness of the AGN selection technique of Lacy et al. (2007) rises from % compared to X-ray selection at X-ray luminosities to 100% complete at (i.e. 5m luminosities to ). Nevertheless, mid-infrared selection finds many AGN that are missing from X-ray surveys even at low luminosities (e.g. Park et al. 2008), so mid-infrared remains a powerful technique for finding AGN.

Quantifying the abundance of the obscured population is important for two reasons. First, it affects the argument of Soltan (1982), where the total amount of matter accreting onto black holes is compared to the remnant black hole mass density today. The comparison allows us to estimate the mean value of the radiative efficiency of accretion, , which varies from for a non-rotating (Schwartzshild) black hole to for a maximally-rotating Kerr black hole. Previous estimates suggest (e.g. Hopkins, Richards & Hernquist 2007, hereafter H07; Martínez-Sansigre & Taylor 2009; Delvecchio et al. 2014), though an analysis of the X-ray background (Elvis, Risaliti & Zamorani 2002) suggests a higher bound, , consistent with most black holes spinning and with current AGN censuses being significantly incomplete. The second reason is to aid in our understanding of unified schemes of AGN. In the standard orientation-based unified scheme (e.g. Antonucci 1993), the ratio of obscured to unobscured quasars defines the opening angle of the obscuring torus. (In a refinement of this model, Elitzur (2012) suggests that the clumpy nature of the torus means that on an object-by-object basis the angle to the line of sight is only a statistical predictor of obscuration, nevertheless the obscured to unobscured ratio is still related to the torus opening angle in a statistical sense.) The luminosity dependence of the obscured to unobscured ratio (Lawrence 1991; Simpson 2005; Gilli, Comastri & Hasinger 2007; Lusso et al. 2013) has been put forward as evidence of the “receding torus” model in which more luminous objects are able to move the inner radius of the obscuring torus out, reducing the covering factor of the dust. Competing with this theory, in the case of high luminosity and (especially) high redshift quasars is the evolutionary theory, whereby quasars begin life highly obscured as a result of mergers of dusty galaxies and eventually break out from their cocoons of dust and gas due to outflows (Sanders et al. 1988; Hopkins et al. 2008). In this case the ratio of obscured to unobscured objects can be used to estimate the fraction of its lifetime that a quasar exists in the obscured state.

To date, estimates of the fraction of dust-obscured quasars have been affected by either small sample size (Martínez-Sansigre et al. 2005; Lacy et al. 2007), a restricted redshift range based on a requirement that [Oiii] be detected in an optical spectrum (e.g. Reyes et al. 2008), or a very specific infrared color selection (Assef et al. 2014). All these estimates suggest that the obscured population of quasars equals or exceeds the unobscured population in space density, even at high luminosities. To increase the sample size available to us, we have undertaken a larger survey (detailed in Lacy et al. 2013; hereafter Paper I). This survey is drawn from the 54 deg covered by the Spitzer Wide-area Infrared Extragalactic Survey (SWIRE; Lonsdale et al. 2003) and the Spitzer Extragalatic First Look Survey (XFLS; Lacy et al. 2005; Fadda et al. 2006). AGN and quasars were selected on the basis of their mid-infrared properties from several samples flux-limited at 24m, with flux density limits ranging from 0.61-9.6mJy. Further selection based on mid-infrared colors at 3.6-8.0m was then used to remove most non-AGN from the survey and provide a candidate list of 786 objects for which optical and/or near-infrared spectroscopy was obtained. Of these, 672 were confirmed spectroscopically as AGN or quasars. The remainder either showed no firm evidence for an AGN in their spectra, had featureless spectra, or had spectra with only a single weak emission line. In this, the second paper in the series, we compare the fractions of obscured and unobscured quasars in this survey, derive the luminosity functions of the different AGN types, and compare to AGN surveys in other wavebands. We assume a cosmology with , and .

## 2. The evolution of the AGN population as seen in the mid-infrared

### 2.1. Photometric data

We used the optical through near-infrared photometric data tabulated in Table 2 of Paper I, combined with 12m flux densities from the Wide field Infrared Survey Explorer (WISE) (Wright et al. 2010). The 12m flux densities were particularly useful for constraining the mid-infrared SEDs of the AGN between the Spitzer 8m and 24m bands. We matched the objects in Paper I to the WISE all-sky survey using a match radius, recovering matches to 896 out of 963 (93%) of objects. (Some of those not matched had faint detections in WISE 12m from which we could estimate flux densities.) These data were used as the basis of the spectral energy distribution (SED) fits described in the Appendix.

### 2.2. The luminosity-redshift plane

Figure 1 shows the redshift-luminosity plane for our survey. The rest-frame 5m luminosity, was calculated by fitting the spectral energy distribution of each object with either a pure type-1 AGN model, or a combination of an AGN and a host stellar population, as detailed in the Appendix. The presence of infrared emission from the host galaxy may lead to some subtle selection effects, which are also discussed further in the Appendix, though we do not believe they have a significant effect on our derived luminosity function, at least at .

The rest-frame 5m luminosity of the AGN (without a stellar host contribution) was then used to calculate . The 14 samples from which the survey was drawn, with their differing 24m flux density limits and areas, account for the relatively wide luminosity range at a given redshift (a factor of ), in all but the highest redshift bins. For the analysis presented in this paper, we use the “statistical sample” of Paper I, which comprises the 662 objects in a 90% spectroscopically-complete subsample of the entire survey (479 of which are confirmed AGN). The “statistical sample” was defined by exploiting the close correlation of and emission line flux density. Objects in each of the 14 samples were sorted in order of decreasing . The value of at which the survey completeness fell below 90% was then noted, and objects with brighter than this included in the “statistical sample”. This procedure prevents the samples becoming too biassed towards objects which are easy to obtain spectra for (e.g. normal type-1 quasars). As described in Paper I, we classify our AGN as normal, blue type-1 AGN, red type-1 AGN (dust reddened, but still showing broad lines in the rest-frame optical, i.e. with towards the AGN), type-2 AGN (showing narrow lines only, even in the rest-frame optical, and thus with towards the AGN), non-AGN (selected as AGN candidates but showing no signs of AGN activity in their optical spectra), stars, and objects with featureless spectra. Our statistical sample consists of 122 normal type-1 AGN, 86 red type-1’s and 271 type-2’s, 118 non-AGN, 23 stars and 42 with featureless spectra.

### 2.3. The luminosity function

We performed maximum likelihood fits both to the luminosity functions of the AGN population as a whole, and to the individual luminosity functions for the normal (blue) type-1s, type-2s and red type-1s. From an initial estimate of the luminosity function, we noted that, although a power-law provided a visually reasonable fit to the data, there was a small, but significant flattening at low luminosities, consistent with the luminosity function being a broken power-law undergoing luminosity evolution (although our data do not probe below the break at , so the faint-end evolution is unconstrained at high redshifts).

We therefore fit a double power-law function of the form:

 ϕ(L,z)=ϕ∗(L/L∗)γ01+(L/L∗)γ02, (1)

where is the comoving space density of AGN, is the characteristic space density, both in units of comoving Mpc, is the rest-frame luminosity at 5m and is the break luminosity at 5m, both in units of . The evolution in is the usual cubic expression (e.g. H07):

 log10L∗(z)=log10L∗0+k1ϵ+k2ϵ2+k3ϵ3, (2)

where , is a free parameter in the fit, corresponding to the break luminosity at redshift zero, is the faint-end slope and is fixed at . We also fit a model with a fixed faint-end by setting in the left-hand term in the denominator of Equation 1 and allowing the bright-end slope, to vary as

 γ2=γ02+τϵ. (3)

To perform the actual fits, we evaluated the log-likelihood function, , in the usual way (Marshall et al. 1983), by minimizing (using the Powell minimization routine from the scipy.optimize library9):

 S=−2∑i=1,naln[ϕ(Li,zi)]+2∑j=1,nf∫LmaxLmin∫ϕ(L,z)AjdVdzdzdL (4)

where the first term sums over the AGN being considered, and the second sums over the samples, each of area . is the minimum luminosity detectable at redshift in a given sample, and the maximum luminosity (for samples where there is an upper flux density limit). For the k-corrections in this calculation, mean [5.8]-[8.0] and [8.0]-[24] spectral indices for each AGN type were measured and applied as appropriate to the AGN type and the redshift range of the calculation.

The fitting was performed by optimizing the parameters , , , , , and for AGN with and . We fit all three AGN types separately (normal type-1s, type-2s and red type-1s). We also fit the type-2s and red type-1s together to give a luminosity function for the total obscured population. In order to assess the possible effects of incompleteness, we fit to a set including the objects classified optically as non-AGN (some fraction of which are, in fact, AGN, see Paper I) and also containing the featureless spectrum objects which have no redshifts, distributed evenly in comoving volume between , and refer to this as the “Maximal” fit. The best-fit parameters are given in Table 1.

For the purposes of comparison with previous work, we also fit a pure power law of the form

 ϕpl=ϕ∗(L/L∗[z])γ2 (5)

where is allowed to vary according to Equation 3.

To determine whether the luminosity evolution model of Equation 1 was a better fit than the fixed faint-end slope model, or the pure power-law model of Equation 5, we examined the likelihood ratio of the fits. For two fits giving values of and using models with and degrees of freedom, respectively, is distributed as a distribution with degrees of freedom (e.g. Freeman et al. 1998). On this basis, the fixed faint-end slope model is slightly, but not significantly worse than the luminosity evolution model (, with one degree of freedom). The smallness of this difference is perhaps to be expected given the poor constraints on the faint end and its evolution. for the pure power-law versus the broken power law model is 30 (Table 1), and the difference in degrees of freedom is only one, so the broken power-law is formally a much better fit.

To make the binned estimates, we binned the data in five redshift bins (, , , and ) assigning a minimum luminosity to each bin at which the bin was approximately complete (unshaded region in Figure 1). The lower redshift limit was chosen to avoid , where the Lacy et al. (2004) AGN selection can be confused by the 6.2m Polycyclic Aromatic Hydrocarbon (PAH) emission. The bin spacing is in intervals of .

The values of the binned estimates were obtained using the formalism of Hasinger, Miyaji & Schmidt (2005), where the fitted luminosity function of Equation (1) is used to predict the number of sources in each bin, by integrating the luminosity function over the luminosity and redshift ranges for each bin for each individual field in the survey, properly accounting for the upper and lower flux limits in each field. The value of the luminosity function at the bin centre, , is then corrected by the ratio of the number of predicted sources in the bin to the actual number observed, - i.e.:

 ϕijb(Li,zj)=ψmdlnobsnpred.

Although this has the disadvantage of making the binned estimated somewhat model dependent, the advantage of this technique is that it compensates both for the tendency of objects to be found close to the bottom of a luminosity bin, and at the higher redshift end of redshift bins where the luminosity function is increasing rapidly with redshift (i.e. at ).

### 2.4. Evolution of the AGN number density with redshift

Figure 2 shows the evolution of the luminosity function of all spectroscopically-confirmed Spitzer AGN in our statistical sample with and . The left-hand panel, Figure 2(a), shows that the luminosity evolution model seems to fit the data well across most of the range of redshift and luminosities. In the center panel, (b) we show the alternative fixed faint-end slope and simple power-law models, both of which seem to represent the data reasonably well (though the single power-law is, statistically, a worse fit than either of the broken power-law models). On the right-hand panel (c) we show the “Maximal” model compared to the luminosity evolution one. Deeper surveys, with AGN identifications at , are needed to determine the faint end behaviour and evolution, though Juneau et al. (2013) show that the faint end of the luminosity function seems to have evolved little since . The evolution at is relatively poorly constrained, but, as described in Section 2.6, there is evidence for a turnover in number densities of high luminosity AGN at . AGN at the break luminosity contribute most of the luminosity density at a given redshift. Our finding that the break luminosity is a strong function of redshift, decreasing with cosmic time at is consistent with the “downsizing” scenario of galaxy evolution (Cowie et al. 1996), in the sense that, for a given mean Eddington rate, the characteristic mass of the black holes contributing the bulk of the AGN luminosity density is decreasing with cosmic time.

When the optically-classified “non-AGN” and blank spectra with randomly-assigned redshifts are added into the survey to produce the Maximal luminosity function, as shown in the right-hand panel of Figure 2 as the dot-dashed line, we see an increase in the magnitude of the faint-end slope. At high luminosities, above the break, however, we see little more than an increase in the normalization (as expected). The larger numbers of low- AGN is consistent with the population of optically-classified “non-AGN” being in fact dominated by low luminosity, heavily obscured AGN. These conclusions must remain speculative, however, as it is very difficult to rule out the possibility that the some of the “non-AGN” are interlopers.

### 2.5. Comparison with other AGN luminosity functions

In Figure 3, we plot a comparison of our mid-infrared AGN luminosity function with the hard X-ray luminosity functions from LaFranca et al. (2005) and Aird et al. (2010), the mid-infrared luminosity function from Matute et al. (2006; based on surveys using the Infrared Space Observatory) at , and the estimate of the bolometric luminosity function of H07. In Figure 4 we plot the normal (blue) type-1 quasar luminosity function from our survey, and compare to both the combined Sloan Digital Sky Survey (SDSS) and 2-degree Field (2dF) quasar surveys luminosity function (Croom et al. 2009), and the mid-infrared luminosity for type-1 AGN from Brown et al. (2006) at .

For both simplicity and ease of comparison, a single bolometric correction, independent of redshift and luminosity, was assumed for each population (Table 2). For our mid-infrared AGN, we assume a uniform bolometric correction of 10. This number was derived as follows. We began with the bolometric correction from 5m to the total type-1 quasar emission, from Richards et al. (2006). We then corrected for unobscured type-1s to allow for the fact that the mid-infrared emission is likely to be much more isotropic than the optical/UV, and have its origin in reprocessed UV light originally emitted away from our line of sight. We thus derive a “true” bolometric correction for type-1 objects, , where 0.6 is the typical ratio of the UV/optical to the total quasar emission. Considering the type-2s on the other hand, the lower mid-infrared luminosities of type-2 objects at a given radio luminosity (Paper I, section 6) suggests that type-2 objects have higher mid-infrared extinctions than the type-1s. (This is also manifest as a steeper 8-24m spectral index for the type-2s, compared to for the type-1s.) Residual extinction at mid-infrared wavelengths results in an underestimate of the mid-infrared luminosity for the type-2s by a factor . Thus the effective value of the bolometric correction for the type-2 population will be . For simplicity, we average these two corrections to obtain a value of . However, we will return to this issue in more detail in future work, where we will undertake a detailed comparison of type-1 and type-2 AGN SEDs.

At low redshifts (), we expect to be incomplete due to contaminating emission from PAHs of up to % (Paper I, section 2.2), and indeed we find that at we are below most other AGN luminosity functions in our bin, including the mid-infrared one of Matute et al. (2006). At , however, the luminosity function matches those in other wavebands much better, and we see a small excess over the X-ray luminosity function at low luminosities (), and a very large one over the optical quasar luminosity function. Considering the blue type-1 objects in our sample alone, however, Figure 4 shows that we agree very well with prior estimates of the normal type-1 luminosity function at all redshifts.

As is already well-known (Hasinger 2008; Aird et al. 2010), the hard X-ray luminosity function peaks at much lower redshift than either the mid-infrared or the optical quasar luminosity function. This difference between the X-ray and other luminosity functions has often been ascribed to a luminosity-dependent bolometric correction (e.g. Steffen et al. 2006; H07, Han et al. 2012). Any redshift-dependent contribution due to absorption effects has been thought to be low, though it has long been known that X-ray absorbers are commonly seen in individual high redshift quasars (Elvis et al. 1994b). The partial correlation analysis of Steffen et al. (2006) in particular seems to indicate a stronger luminosity than redshift dependence on the ratio between optical and X-ray fluxes, at a given optical or X-ray luminosity (though for X-ray luminosities the result is only significant at 3), which would translate into a difference in the bright-end slope of the luminosity function. However, we find no evidence that a strong luminosity-dependent bolometric correction is required based on our comparison of luminosity functions, and we suspect that the apparently high (13.6) significance of the anticorrelation on optical luminosity is a by-product of luminosity-dependent obscuration effects, i.e. AGN with lower luminosities tend to have higher gas columns in the X-ray, increasing the magnitude of . In Figure 3, we see that, at , the mid-infrared, optical quasar, and hard X-ray luminosity functions are very similar at the high luminosity () end of the luminosity function. Instead of a difference in the slope of the AGN luminosity functions between infrared and X-ray, which is what we would expect if there was a strong luminosity dependence to the bolometric correction, what we see is a fall-off in the hard X-ray AGN population at all luminosities (particularly when the Aird et al. luminosity function is considered, though this luminosity function is thought to be somewhat low at (Kalfountzou et al. 2014)). This suggests that the difference in the luminosity functions is predominantly a redshift effect. The lack of a luminosity-dependent bolometric correction in the X-ray is particularly illustrated when the bolometric luminosity function of H07 (which assumes a luminosity-dependent bolometric correction) is plotted in Figure 3, when we see a significant overprediction of the space density of AGN at .

### 2.6. Evolution by type of AGN

Figure 5 shows the evolution of the mid-infrared AGN population split up by type. The top two panels, (a) and (b), show the luminosity functions at mean redshifts of 0.65 and 1.25, respectively, for the different types of AGN, with the plotted points showing the binned estimates. The black dotted line indicates the fit for all confirmed AGN, with the cyan line that for the Maximal model. Compared to the type-1 luminosity functions (shown in blue and red, for normal and red type-1s, respectively), the type-2 luminosity function (green) has a much steeper faint-end slope. In the lower two panels, (c) and (d), the number density of AGN above luminosities of and , respectively, are plotted against redshift. Up to all the populations shown (all AGN, blue and red type-1s, type-2s and the combination of red type-1s and type-2s, representing the obscured population) show a similar evolution. The blue type-1 population shows the same evolution as optically-selected quasars, with a peak in number density at , but the obscured quasar population (magenta line in panel (d)) shows a peak in number density that is broader and shifted to earlier epochs (). This is consistent with some theoretical predictions (Fanidakis et al. 2012), results from studies of the radio galaxy population (Jarvis et al. 2001), and a study of the infrared luminosity function of quasars from the Sloan Digital Sky Survey (Vardanyan, Weedman & Sargsyan 2014). Our models are fairly poorly constrained beyond the peak, however, with only 13 objects of any type at , so the reality of the difference at high redshifts may be affected by small number statistics. Nevertheless, we have performed some tests to investigate whether the difference in the evolution we see might be statistically significant. We take for a null hypothesis that the evolution of obscured AGN (type-2 and red type-1s) and the unobscured (type-1) population is the same. We then compare the values for fits which fix the evolution of the obscured population to that of the unobscured one, and vice-versa. Fixing the values of the evolution parameters and for the obscured population to those fit for the unobscured population, and allowing the remainder of the parameters to be fit results in a difference in the values of 17, which, with the three degrees of freedom from fixing the parameter values allows us to eliminate the null hypothesis that the evolution is the same at % confidence. Forcing the evolution terms for the unobscured population to be the same as those for the obscured one results in a less significant , % confidence. We also compared to the evolution of all AGN. As expected, fixing the evolution of the obscured AGN to match that of all AGN resulted in no significant difference, but the difference from forcing the type-1 evolution to match that of all AGN resulted in ruling out the null hypothesis at the 95% confidence level.

In summary, we find that the evolution of the obscured and unobscured AGN populations are significantly different, at % confidence using our strongest statistical test. The most obvious difference in the luminosity function evolution is at high redshifts, where the number density of high luminosity obscured quasars appears to peak at a higher redshift than that of the unobscured population. Small number statistics, however, mean this interpretation should be treated with caution.

### 2.7. The obscured fraction as a function of luminosity and redshift

Figure 6 shows histograms of number of objects as a function of redshift and luminosity, split up by type. Figure 7 plots the fraction of obscured objects as a function of both luminosity and redshift with the measured points plotted with error bars and the luminosity function models plotted as dashed lines. These plots show that the fraction of obscured objects varies from % for Seyfert galaxies to % at the highest quasar luminosities. In addition, there is evidence that the nature of the obscuration in the survey seems to change as a function of redshift and/or luminosity. At low luminosities and redshifts, most of the obscured objects are highly obscured, type-2, objects, but at high luminosities and redshifts the obscured population is dominated by less obscured red type-1 objects. The right-hand panel suggests that, despite the change in the nature of the obscuration, the overall obscured fraction is not a strong function of redshift when comparing bins at and , and indeed, as discussed in Section 2.6 below, it is not until that we see a possible difference. Of course selection effects may well be playing a role here as heavily obscured objects become progressively harder to identify at high redshifts even in samples selected in the mid-infrared. As we discuss in Paper I, however, we believe such selection effects are not having a severe effect on our sample, even at .

### 2.8. The AGN luminosity density and the black hole mass density

Figure 8 shows the AGN bolometric luminosity density and corresponding accreted black hole mass density since . Panel (a) shows the luminosity density as a function of redshift, with dusty type-1 objects dominating at high redshifts, normal type-1s at and type-2 objects at . As our best fit luminosity evolution model has a steep faint-end slope, we need to assume a low luminosity cutoff at high luminosities in order to prevent an excessive amount of AGN luminosity density from being estimated. We assume, as a fiducial model, that the fitted form of the luminosity function we have adopted can be extrapolated a factor of three below the current limit of the survey. This results in luminosity densities in line with those we obtain from our alternative fixed faint-end slope model. With this assumption, we find more AGN than previous surveys, and integrating to the present day (e.g. Martínez-Sansigre & Taylor 2009) indicates a total energy, , emitted by all AGN over the history of the Universe of , where the lower bound comes from cutting off the faint-end of the luminosity function at the luminosity corresponding to the flux limit of the deepest sample in our survey (0.61mJy at 24m), and the upper bound from extrapolating the luminosity function to ten times fainter than this limit. To obtain the observed energy density today due to AGN, , we need to divide the integrand by (Soltan 1982), resulting in a value of , or a radiation intensity , within the range of estimated from the X-ray background, and suggesting that % of the total radiation intensity of the Universe is contributed by AGN (Elvis et al. 2002). We therefore need a radiative efficiency of accretion, at the high end of previous estimates to match the local mass density in black holes (, Elvis et al. (2002); , Martínez-Sansigre & Taylor 2009) (Figure 8 (b), (c)). We use the compendium of local black hole mass density estimates of Graham & Driver (2007), with a mean of to estimate a best fit value of .

As a function of cosmic time, the evolution of the black hole mass density shows a very rapid rise with at , then continues to rise, but at a smaller rate. Fifty percent of the current mass density has accreted by . At , we predict that most black hole mass accretion occurs in dusty type-1 AGN, this switches to normal type-1 quasars at , and to type-2 AGN at when powerful quasars become very rare.

## 3. Discussion

We present evidence that, when largely free from the effects of dust obscuration, the AGN luminosity function for actively-accreting black holes seems to be well represented by a double power-law with a break that is a strong function of luminosity, though the form of the faint end remains unconstrained at . Downsizing, in the sense that the characteristic AGN luminosity at which most of the AGN luminosity density is contributed is a strongly increasing function of redshift, is clearly present at . Luminosity functions are often difficult to interpret in terms of physical models; however, at low luminosities, we know that most local AGN are not triggered by major merger events (e.g. Kocevski et al. 2012). At high luminosities, especially in the obscured population, major mergers may indeed play the dominant role (e.g. Urrutia et al. 2008; Treister et al. 2012; Hopkins, Kocevski & Bundy 2013). We therefore speculate that the increase in the break luminosity with redshift is driven by the increased frequency of major mergers of gas-rich galaxies at high redshifts.

A comparison of the mid-infrared luminosity function with those from other wavelengths indicates that we are seeing many more dust-obscured AGN than the normal type-1 AGN population, and more than the hard X-ray population, particularly at high redshifts. We find that the evolution of the obscured and the unobscured AGN population are significantly different, and we see tentative evidence for an earlier peak in the space density of obscured quasars, consistent with that seen in the radio-loud population. Our study then is consistent with the evolutionary explanation for obscured AGN, whereby AGN begin heavily obscured and outflows clear out material from the line of sight to the AGN.

Our new luminosity function has allowed us to improve the constraints on the luminosity density produced by AGN, and on the integral of this over cosmic time. Our overall AGN luminosity density agrees well with the range predicted from the X-ray background (Elvis et al. 2002), suggesting that we have indeed accounted for the vast majority of AGN in the Universe, and that about 12% of the luminosity density of the Universe is contributed by AGN. When compared to the mass density in compact central objects in galaxies today, this luminosity density allows an estimate of the mean radiative efficiency of accretion onto them. This number constitutes one of the best constraints we have on the nature of the compact central objects in galaxies. Although it is usually assumed that these are black holes, there are some theoretical alternatives, such as gravastars, that can deliver comparable radiative efficiencies (Harko, Kovács & Lobo 2009). However, it seems likely that the radiative efficiencies of such objects are (Bambi 2012), comparable to the current limit from the X-ray background (Elvis et al. 2002), but below our best estimate of . Significant systematic uncertainties remain, however. AGN missing from our survey, because they are not confirmable in other wavebands (such as the “non-AGN” objects in this paper) or are missing from the mid-infrared selection criteria (such as low luminosity AGN) will only serve to increase our estimate of the radiative efficiency. However, there remain uncertainties on the local black hole mass density, the appropriate bolometric corrections to apply to both obscured and unobscured objects, and the low luminosity end of the luminosity function at high redshifts that could result in either an increase or a lowering of our estimate for . The bolometric correction uncertainty can be addressed by using a full multiwavelength dataset, including archival Herschel data, to obtain accurate template SEDs for the obscured AGN population. Improved depth can be obtained by using more sophisticated AGN selection techniques based on the complete infrared SED to exploit the full depth of SWIRE to select AGN up to four times fainter than the current limit, and photometric redshifts based on new deep near-infrared surveys in the SWIRE fields (SERVS, Mauduit et al. 2012; VIDEO, Jarvis et al. 2013, and UKIDSS DXS, Lawrence et al. 2007) to obtain redshifts for these AGN, whose emission lines can be too faint to detect even using 8m class telescopes.

We thank the referee whose comments significantly improved the paper. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

## Appendix A Optical through mid-infrared SED fitting

In order to disentangle the contribution of host galaxy light from our AGN fluxes, we performed SED fitting on the photometry described in Section 2.1, using the modeling techniques of Sajina et al. (2006, 2012, see also Lacy et al. 2007b and Hiner et al. 2011). We fit the SEDs at wavelengths from optical through 24m using the photometry described in Section 2.1. For unobscured and red type-1s we fit a pure quasar template from Richards et al. (2006) (with an SMC reddening law applied in the case of the dusty red quasars, and a host galaxy stellar population added in a few low luminosity objects). These appear to work fairly well, apart from the 1-5m region where the AGN SED is known to vary as a function of luminosity (e.g. Gallagher et al. 2007). For type-2s, we used the mid-infrared AGN template of Sajina et al. (2012), and a host galaxy stellar population, modeled as a single stellar population from Bruzual & Charlot (2003) with an age of 0.1,0.64, 1.4 or 5 Gyr picked to best fit the optical SED. 0.64 Gyr was used as the default in cases where the optical/near-IR data were too sparse to distinguish between models (the vast majority of cases), however, 21 were fit with with a 1.4Gyr population, 23 with a 100Myr population, and ten with a 5Gyr population. All the stellar population models had similar behavior in the near-infrared region of the spectrum, however, so the exact choice of stellar age does not affect the AGN flux density at rest-frame 5m by more than a few percent. Examples of our fits are shown in Figure 9. Full SED fits through the far-infrared will be discussed in detail in a future paper.

The SED fitting highlights some ways in which our sample may be subject to subtle selection effects due to the prescence of mid-infrared emission from the host galaxy, particularly at low luminosities and redshifts. Our use of flux limits at 24m ensures that photospheric emission from stars contibutes negligibly (% ) to the flux density at the selection frequency. However, there is the potential for some objects at low redshifts to have their 24m flux densities boosted by emission from warm dust due to star formation. Although this will not contribute significantly to the AGN flux density estimated at 5m, it may lead to objects being included in the sample that would otherwise have fallen below the flux density limit at 24m. However, these same objects are, to a greater extent, selected against in the survey as their 7.7m PAH features will lift them out of the AGN selection region at (Section 2.5). As the slope of the warm dust SED is very steep at wavelengths between m and 24m, the contribution of warm dust emission to the observed 24m flux density at higher redshifts is likely to be less than a few percent.

### Footnotes

1. affiliation: National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, Virginia 22903
2. affiliation: National Optical Astronomy Observatory, 950 North Cherry Avenue, Tucson, AZ 85719
3. affiliation: Department of Physics and Astronomy, Tufts University, 212 College Avenue, Medford, MA 02155
4. affiliation: Gemini Observatory, Northern Operations Center, 670 N. A’ohoku Place, Hilo, HI 96720
5. affiliation: UCO/Lick Observatory, P.O. Box 85, Mount Hamilton, CA 95140
6. affiliation: Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16, 14482, Potsdam, Germany
7. affiliation: Spitzer Science Center, Caltech, Mail Code 220-6, Pasadena, CA 91125; mlacy@ipac.caltech.edu,sajina@ipac.caltech.edu, lisa@ipac.caltech.edu
8. slugcomment:
9. http://scipy.org

### References

1. Aird, J., Nandra, K., Laird, E.S. et al. 2010, MNRAS, 401, 2531 (A10)
2. Alexandroff, R., Strauss, M.A., Greene, J.E., Zakamska, N.L., Ross, N.P. et al. 2013, MNRAS, 435, 3306
3. Antonucci, 1993, ARA&A, 31, 473
4. Assef, R.J., Eisenhardt, P.R.M., Stern, D. et al. 2014, ApJ, submitted (arXiv:1408.1092)
5. Bambi, C. 2011, PhLB, 705, 5
6. Brown, M.J.I., Brand, K., Dey, A. et al. ApJ, 638, 88 (B06)
7. Cowie, L.L., Songaila, A., Hu, E.M. & Cohen, J.G. 1996, AJ, 112, 839
8. Croom, S.M., Richards, G.T., Shanks. T. et al. 2009, MNRAS, 399, 1755
9. Vardanyan, V., Weedman, D. & Sargsyan, L. 2014, ApJ, 790, 88
10. Delvecchio, I., Gruppioni, C., Pozzi, F. et al. 2014, MNRAS, in press
11. Donley, J.L., Koekemoer, A.M., Brusa, M. et al. 2012, ApJ, 748, 142
12. Elitzur, M. 2012, ApJ, 747, L33
13. Elvis, M., Wilkes, B.J., McDowell, J.C. et al. 1994a, ApJS, 95, 1
14. Elvis, M., Fiore, F., Mathus, S. & Wilkes, B.J. 1994, ApJ, 425, 103
15. Elvis, M., Risaliti, G. & Zamorani, G. 2002, ApJ, 565, L75
16. Fadda, D., Marleau, F., Storrie-Lombardi, L.J. et al. 2006, AJ, 131, 2859
17. Fanidakis, N., Baugh, C.M., Benson, A.J. et al. 2012, MNRAS, 419, 2797
18. Freeman, P.E., Graziani, C., Lamb, D.Q. et al. 1999, ApJ, 524, 753
19. Gallagher, S. C., Richards, G. T., Lacy, M., et al. 2007, ApJ, 661, 30
20. Gilli, R., Comastri, A. & Hasinger, G. 2007, A&A, 463, 79
21. Graham, A.W. & Driver, S.P. 2007, MNRAS, 380, L15
22. Harko, T., Kovács, Z. & Lobo, F.S.N. 2009, Classical and Quantum Gravity, 26, 215006
23. Han, Y., et al. 2012, MNRAS, 423, 464.
24. Hasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441, 417
25. Hasinger, G. 2008, A&A, 490, 905
26. Hickox, R.C., Jones, C., Forman, W.R. et al. 2009, ApJ, 696, 891
27. Hopkins, P.F., Richards, G.T. & Hernquist, L. 2007, ApJ, 654, 731 (H07)
28. Hopkins, P.F., Kocevski, D.D. & Bundy, K. 2013, MNRAS, submitted, arXiv1309.6321
29. Jarvis, M.J., Rawlings, S., Willott, C.J., Blundell, K.M. & Lacy, M. 2001, MNRAS, 327, 907
30. Jarvis, M. J., Bonfield, D. G., Bruce, V. A., et al. 2013, MNRAS, 428, 1281
31. Juneau, S., Dickinson, M., Alexander, D.M. & Salim, S. 2011, ApJ, 736, 104
32. Juneau, S., Dickinson, M., Bournaud, F., Alexander, D.M., Daddi, E. et al. 2013, ApJ, 764, 176
33. Kalfountzou, E., Civano, F. M., Elvis, M. & Trichas, M. 2014, American Astronomical Society Meeting Abstracts, 223, #150.05
34. Lacy, M., Storrie-Lombardi, L.J., Sajina, A. et al. 2004, ApJS, 154, 166
35. Lacy, M., Wilson, G., Masci, F. et al. 2005, ApJS, 161, 41
36. Lacy, M., Petric, A.O., Sajina, A. et al. 2007, AJ, 133, 186
37. Lawrence, A. 1991, MNRAS, 248, 91
38. Lawrence, A. et al. 2007 MNRAS 379, 1599
39. Lusso, E., Hennawi, J.F., Comastri, A. et al. ApJ, 777, 86
40. Marshall, H.L., Avni, Y., Tannenbaum, H. & Zamorani, G. 1983, ApJ, 269, 35
41. Martínez-Sansigre, A., Rawlings, S., Lacy, M. et al. 2005, Natur, 436, 666
42. Martínez-Sansigre, A. & Taylor, A.M. 2009, ApJ, 692, 964
43. Matute, I., La Franca, F., Pozzi, F. et al. 2006, A&A, 451, 443 (M06)
44. Mauduit, J.-C., Lacy, M., Farrah, D. et al. 2012 PASP 124, 1135
45. Mendez, A., Coil, A.L., Aird, J. et al. 2013, ApJ, 770, 40
46. Messias, H., Afonso, J.M., Salvato, M., Mobasher, B. & Hopkins, A.M. 2014, A&A, 562, 144
47. Park, S.Q., Barmby, P., Fazio, G.G. et al. 2008, ApJ, 678, 744
48. Park, S.Q., Barmby, P., Willner, S.P. et al. 2010, ApJ, 717, 1181
49. Reyes, R., Zakamska, N.L., Strauss, M.A. et al. 2008, AJ, 126, 2373
50. Richards, G.T., Lacy, M., Storrie-Lombardi, L.J. et al. 2006, ApJS, 166, 470
51. Roseboom, I.G., Lawrence, A., Elvis, M., Petty, S., Shen, Y. & Hao, H. 2013, MNRAS, 429, 1494
52. Simpson, C. et al. 2006, MNRAS, 372, 741
53. Soltan, A. 1982, MNRAS, 200, 115
54. Stern, D. et al. 2005, ApJ, 631, 163
55. Thorne, K.S. 1974, ApJ 191, 507
56. Treister, E., Schawinski, K., Urry, C.M., Simmons, B.D. 2012, ApJ, 758, 39
57. Urrutia, T., Lacy, M. & Becker, R.H. 2008, ApJ, 674, 80
58. Willott, C.J., Rawlings, S., Blundell, K.M. & Lacy, M. 2000, MNRAS, 322, 356
59. Wright, E.L., Eisenhardt, P.R.M., Mainzer, A.K. et al. 2010, AJ, 140, 1868
60. Zakamska, N.L. et al. 2003, AJ, 126, 2125
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters