Observational constraints on the specific accretion-rate distribution of X-ray selected AGN

TL;DR

This study analyzes the distribution and evolution of specific accretion rates in X-ray selected AGN, revealing trends in AGN activity, duty cycle evolution, and their relation to galaxy properties across cosmic time.

Contribution

It provides the first comprehensive observational constraints on the specific accretion-rate distribution of AGN across different redshifts and stellar masses.

Findings

Probability of hosting an AGN increases with decreasing specific accretion rate.

Evidence of a reversal in trend at very low accretion rates.

Steep decline in accretion probability above the Eddington limit.

Abstract

This paper estimates the specific accretion-rate distribution of AGN using a sample of 4821 X-ray sources from both deep and shallow surveys. The specific accretion-rate distribution is defined as the probability of a galaxy with a given stellar mass and redshift hosting an active nucleus with a certain specific accretion rate. We find that the probability of a galaxy hosting an AGN increases with decreasing specific accretion rate. There is evidence that this trend reverses at low specific accretion rates, $\lambda<10^{-4}-10^{-3}$ (in Eddington units). There is also a break close to the Eddington limit, above which the probability of an accretion event decreases steeply. The specific accretion-rate distribution evolves such that the fraction of AGN among galaxies drops toward lower redshifts. This decrease in the AGN duty cycle is responsible for the strong evolution of the accretion density of the Universe from redshift $z\approx1-1.5$ to the present day. Our analysis also suggests that this evolution is accompanied by a decoupling of accretion events onto black holes from the formation of stars in galaxies. There is also evidence that at earlier times the relative probability of high vs low specific accretion-rate events among galaxies increases. We argue that this differential redshift evolution of the AGN duty cycle with respect to $\lambda$ produces the AGN downsizing trend, whereby luminous sources peak at earlier epochs compared to less luminous ones. Finally, we also find a stellar-mass dependence of the specific accretion-rate distribution, with more massive galaxies avoiding high specific accretion-rate events.

Figures

Tables

Table 1. Table 1: Number of X-ray sources in the full-band selected sample

field	solid angle	Number of	Number of	Number of	photometric	spectroscopic
	(${\mathrm{deg}}^2$)	X-ray sources	optical/infrared IDs	sources used	redshifts	redshift
(1)	(2)	(3)	(4)	(5)	(6)	(7)
4Ms CDFS	$7.53\times {10}^{-2}$	417	413(1)	416	141	271
AEGIS-XD	$25.96\times {10}^{-2}$	818	793 (10)	805	475	308
COSMOS Legacy	1.53	2676	2641 (18)	2658	1087 (875)	1536
XMM-XXL	3.19	682	605 (0)	605	199	406 (206)
XMMSL	1165	479	440 (82)	337	4	333 (253)

Table 2. Table 2: cigale parameter settings for the determination of stellar masses of X-ray selected AGN.

Parameter description	values/range
Stellar population parameters
Initial Mass Function	Chabrier (2003)
Metallicity	0.02 (Solar)
Single Stellar Population Library	Bruzual & Charlot (2003)
Delayed star-formation history law parameters
e-folding time, $\tau$, in Gyrs	0.1, 0.3, 0.6, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 7.0, 9.0, 13.0
stellar population age, $t$, in Gyrs	0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0
Dust extinction of stellar light
Dust attenuation law	Calzetti et al. (2000)
Reddening $E\left(B-V\right)$	0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.2
$E\left(B-V\right)$ reduction factor between old and young stellar populations	0.44
Fritz et al. (2006) library parameters
Ratio between outer and inner dust torus radii	30, 100
$9.7\mu m$ equatorial optical depth	0.3, 3.0, 6.0, 10.0
Fritz et al. (2006) density profile parameter ${\beta }^{\star }$	-0.5
Fritz et al. (2006) density profile parameter ${\gamma }^{\star }$	0.0
Viewing angle relative to equatorial axis (deg)	0.001, 50.100, 89.990
AGN emission relative to infrared luminosity ($1-1000\hspace{0.17em}\mu m$)	0.0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 0.9

Table 3. Table 3: Galaxy X-ray luminosity function parameters

parameter	value
(1)	(2)
$\log K$	$-3.59$
$\log L_0$	$41.12$
$\alpha$	$0.81$
$\beta$	$2.66$
$z_c$	$0.82$

Equations

$$P(CR,T)=\frac{E^T{e}^{-E}}{T!}.$$

$$E=CR\times t+B.$$

$$\begin{array}{rl}\mathcal{L}(d_i|\theta )=& e^{-\mu }\times \prod _{i=1}^N\int \mathrm{d}\log L_{\mathrm{X}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\\ & p(L_{\mathrm{X}},z|d_i){\varphi }_{\mathrm{T}}(L_{\mathrm{X}},z|{\theta }_{\mathrm{T}}),\end{array}$$

$${\varphi }_{\mathrm{T}}(L_{\mathrm{X}},z|{\theta }_{\mathrm{T}})=\varphi (L_{\mathrm{X}},z|\theta )+{\varphi }_{\mathrm{gal}}(L_{\mathrm{X}},z|{\theta }_{\mathrm{gal}}),$$

$$\mu =\int \mathrm{d}\log L_{\mathrm{X}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}zA(L_{\mathrm{X}},z){\varphi }_{\mathrm{T}}(L_{\mathrm{X}},z|{\theta }_{\mathrm{T}}),$$

$$\varphi (L_{\mathrm{X}},z)=\int \psi (M,z)P(\lambda (L_{\mathrm{X}},M),z)\mathrm{d}\log M,$$

$$\int P(\lambda (L_{\mathrm{X}},M),z)\mathrm{d}\log M=1.$$

$$\begin{split}\mathcal{L}(d_{i}\;|\;\omega,\theta_{\rm gal})=&e^{-\mu}\times\prod_{i=1}^{N}\Biggl{[}\int\mathrm{d}\log L_{\mathrm{X}}\,\,\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\,\,\mathrm{d}\log M\\ &p(L_{\mathrm{X}},M,z\,|\;d_{i})\;\psi(M,z)\,P(\lambda(L_{\mathrm{X}},M),z\;|\;\mathbf{\omega})\\ &+\int\mathrm{d}\log L_{\mathrm{X}}\,\,\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\\ &p(L_{\mathrm{X}},z|\;d_{i})\;\phi_{\rm gal}(L_{\mathrm{X}},z\;|\;\mathbf{\theta_{\rm gal}})\Biggr{]},\end{split}$$

$$\begin{array}{rl}\mu =& \int \mathrm{d}\log L_{\mathrm{X}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}z\mathrm{d}\log M\times \\ & A(L_{\mathrm{X}},z)\psi (M,z)P(\lambda (L_{\mathrm{X}},M),z|\omega )\\ & +\int \mathrm{d}\log L_{\mathrm{X}}\frac{\mathrm{d}V}{\mathrm{d}z}\mathrm{d}zA(L_{\mathrm{X}},z){\varphi }_{\mathrm{gal}}(L_{\mathrm{X}},z|{\theta }_{\mathrm{gal}}).\end{array}$$

$$\lambda =\frac{25L_X(2-10\mathrm{keV})}{1.26\times 10^{38}0.002M},$$

$$\phi_{\rm gal}(L_{X},z)=K\;\Bigg{(}\frac{L_{X}}{L_{X}^{\star}}\Bigg{)}^{-\alpha}\;e^{-\frac{L_{X}}{L_{X}^{\star}}}.$$

$$\log L_{X}^{\star}(z)=\left\{\begin{array}[]{ll}\log L_{0}+\beta\,\log(1+z),&z<z_{c},\\ &\\ \log L_{0}+\beta\,\log(1+z_{c}),&z\geq z_{c}.\\ \end{array}\right.$$

$$\overline{\lambda (z)}=\frac{{\int }_{{\lambda }_l}^{\inf }P(\lambda ,z)\lambda \mathrm{d}\log \lambda }{{\int }_{{\lambda }_l}^{\inf }P(\lambda ,z)\mathrm{d}\log \lambda },$$

$$\varphi (M,z)={\int }_{L_X>L_{cut}}\psi (M,z)P(\lambda (L_{\mathrm{X}},M),z)\mathrm{d}\log \lambda ,$$

Full text

††pubyear: 2016††pagerange: Observational constraints on the specific accretion-rate distribution of X-ray selected AGN–Observational constraints on the specific accretion-rate distribution of X-ray selected AGN

1 Introduction

X-ray surveys in the last 15 years have provided an excellent census of the Active Galactic Nuclei (AGN) population and describe with increasing accuracy how the space density of these sources depends on redshift, accretion luminosity and amount of line-of-sight obscuration (e.g., Brandt & Alexander, 2015). The diverse studies of the accretion history of the Universe have established a number of observational facts and trends. We now know for example, that about 80% of the black-hole growth in the Universe is taking place behind obscuring gas and dust clouds and that the obscured AGN fraction is higher at earlier cosmic times (e.g., Ueda et al., 2014; Buchner et al., 2015). The space density of AGN increases rapidly from the local Universe to redshift $z\approx 1.5-2$ (Ueda et al., 2003; Hasinger et al., 2005; Aird et al., 2010, 2015), and then shows a strong decline, at least for the luminous part of the population, beyond redshift $z\approx 3$ (e.g., Aird et al., 2008; Vito et al., 2014; Georgakakis et al., 2015; Marchesi et al., 2016b). These evolutionary patterns are not unique to X-ray AGN but are also found in studies of optically selected QSOs (Osmer, 1982; Schmidt et al., 1995; Richards et al., 2006; Croom et al., 2009; Ikeda et al., 2011; Ross et al., 2013; McGreer et al., 2013). X-ray selected AGN also appear to evolve in a anti-hierarchical manner, which is often termed downsizing, whereby more luminous sources have space densities that peak at earlier epochs compared to less luminous ones (Ueda et al., 2003; Hasinger et al., 2005; Barger et al., 2005; Ebrero et al., 2009). Striking similarities are also observed between the global evolution of the accretion-rate density and the star-formation rate density of the Universe (Ueda et al., 2003; Aird et al., 2010; Zheng et al., 2009; Aird et al., 2015).

Work is also ongoing to interpret the evolutionary patterns described above in a physical context. One approach is to study the accretion events onto supermassive black holes in relation to their host galaxies to explore which galaxy types and environments are conducive to black-hole growth. The diversity of AGN host galaxy morphologies for example, provides constraints on the role of major-galaxy mergers or secular processes as triggering mechanisms of the observed nuclear activity (e.g., Georgakakis et al., 2009; Gabor et al., 2009; Cisternas et al., 2011; Kocevski et al., 2012; Cheung et al., 2015; Cisternas et al., 2015). The level of star-formation in AGN hosts is used as a proxy of gas availability to explore the necessary conditions for black-hole accretion events (e.g., Santini et al., 2012; Rovilos et al., 2012; Rosario et al., 2013; Mullaney et al., 2015). The position of AGN on the cosmic web and their clustering properties relative to galaxies may hold important clues on the diversity of black-hole fuelling modes and the role of Mpc-scale environment on accretion events (Allevato et al., 2011; Fanidakis et al., 2013; Allevato et al., 2014).

The studies above explore AGN host galaxy properties in relation to accretion luminosity measured by either X-ray or multiwavelength observations. An alternative important parameter of the accretion process that provides more physical insights into the fuelling mode of AGN is the Eddington ratio, which measures the rate a black hole grows relative to its maximum capacity. Different AGN fuelling modes or triggering mechanisms may produce events with similar accretion luminosities but distinct Eddington-ratio distributions (e.g., Fanidakis et al., 2012; Hirschmann et al., 2014). One approach to place constraints on the Eddington-ratio distribution of AGN as a function of redshift is via the continuity equation of the black-hole mass function using the AGN luminosity function as boundary condition (Merloni & Heinz, 2008; Shankar et al., 2013; Aversa et al., 2015). A direct measurement of the Eddington ratio distribution of AGN populations is limited to broad optical emission-line sources only (e.g., Kelly et al., 2010; Schulze et al., 2015), which are needed to estimate the black-hole mass of individual systems. A relatively recent development that has opened new opportunities in studies of black-hole growth has been the determination of the specific accretion rate of large samples of AGN. This quantity measures the rate a black-hole accretes material relative to the stellar mass of its host galaxy and under certain assumptions it can be viewed as a proxy of the Eddington ratio. The advantage is that the specific accretion rate can be measured for AGN over a wide range of accretion luminosities and levels of line-of-sight obscuration, thereby providing a more representative view of the AGN population. Early observational studies suggested a nearly universal shape of the specific accretion rate distribution of X-ray selected AGN, approximated by a power-law with a slope that does not depend strongly on either redshift or host galaxy stellar mass (Aird et al., 2012), and a break close to the Eddington limit (Bongiorno et al., 2012). However, more recent investigations suggest a more complex picture, where the shape of the specific accretion rate distribution may depend on the redshift, the stellar mass (Bongiorno et al., 2016), and possibly the star-formation rate of AGN hosts (Kauffmann & Heckman, 2009; Georgakakis et al., 2014; Azadi et al., 2015). Therefore a detailed study of how the specific accretion rate distribution of AGN changes with redshift, and/or diagnostics of the physical conditions on larger scales is important to get insights into the processes at play. A limitation of current observational studies outside the local Universe is the small sample size that affects the statistical reliability of the results and the ability to explore trends with redshift.

In this paper we combine the largest sample to date of X-ray sources with stellar mass estimates for their host galaxies to constrain the specific accretion-rate distribution of AGN as a function of redshift. The sample is compiled from X-ray survey fields with different sizes and depths, from deep and pencil-beam, such as the 4 Ms Chandra Deep Field South with an area of about $\rm 0.08\,deg^{2}$ (Xue et al., 2011), to wide and shallow, such as a subset of the XMM Slew Survey (Saxton et al., 2008) with approximately 6 s exposure over about a thousand square degrees. A Bayesian methodology is developed that uses the stellar mass function of the galaxy population as a boundary condition to estimate the specific accretion-rate distribution of AGN. Our analysis is non-parametric, i.e., no shape is imposed on the specific accretion-rate distribution of AGN. This approach also has the advantage that uncertainties in the properties of individual AGN that are used in the analysis, such as stellar masses, redshifts and X-ray luminosities are correctly propagated into the analysis. The combination of samples with diverse X-ray depths and sky areas provide a satisfactory coverage of the luminosity-redshift plane and yields robust constraints on the specific accretion rate distribution of AGN over a range of redshifts and specific accretion rates.

Parallel to the work presented in this paper, Aird et al. (in prep.) have developed an independent approach to infer the specific accretion-rate distribution of black holes within galaxy samples binned in redshift, stellar mass and star-formation rate. Although we share some of the X-ray data with the Aird et al. (in prep.) paper, the two works differ in the supporting multiwavelength photometry and the overall methodology adopted to infer the specific accretion-rate distribution, i.e., the probability of a galaxy hosting an AGN with a certain specific accretion rate. In our work we use ony X-ray selected AGN and require that the convolution of the specific accretion-rate distribution and the galaxy stellar mass function yields the X-ray luminosity function of AGN. Aird et al. (in prep.) start with large galaxy samples, binned by stellar mass and redshift, and use the X-ray data to infer the probability of a galaxy hosting a black hole with a certain specific accretion rate (see also Aird et al., 2017). The results from these two independent approaches are compared in the Appendix of our paper. In the calculations that follow we adopt cosmological parameters $\rm H_{0}=70\,km\,s^{-1}\,Mpc^{-1}$, $\Omega_{M}=0.3$, $\Omega_{\Lambda}=0.7$. This choice is to facilitate comparisons with the results from previous studies.

2 Data Products

2.1 Chandra and XMM-Newton survey fields

The X-ray survey fields used in this paper include the 4 Ms Chandra Deep Field South (CDFS; Xue et al., 2011; Rangel et al., 2013), the deep ($\rm\approx 800\,ks$) Chandra survey of the All-wavelength Extended Groth Strip (AEGIS-XD, Nandra et al., 2015), the recently completed Chandra COSMOS-Legacy survey (Civano et al., 2016) and the subregion of the XMM-XXL (Pierre et al., 2015) equatorial field that overlaps with the VISTA (Visible and Infrared Survey Telescope for Astronomy) Deep Extragalactic Observations survey (VIDEO; Jarvis et al., 2013).

The Chandra observations of the survey fields above are analysed in a homogeneous way by applying the data processing steps and the source detection described in Laird et al. (2009) and Nandra et al. (2015) as well as the sensitivity map generation methods described in Georgakakis et al. (2008). In this paper we use X-ray sources selected in the 0.5-7 keV (full) band with a Poisson false probability threshold of $<4\times 10^{-6}$. The source catalogue of the XMM-XXL equatorial field is presented by Liu et al. (2016) based on the data reduction and analysis pipeline described in Georgakakis & Nandra (2011). We use sources in that field selected in the 2-8 keV (hard) band to the same Poisson false probability threshold as above.

The identification of X-ray sources with optical and mid-infrared counterparts follows the likelihood ratio method (Sutherland & Saunders, 1992; Ciliegi et al., 2003; Brusa et al., 2007). Specific details on the identification of sources in the Chandra survey fields used in this work can be found in Aird et al. (2015). The association of XMM-XXL sources with multiwavelength photometric catalogues is described by Georgakakis et al. (2017) and Menzel et al. (2016). In this paper we use a subregion of the XMM-XXL equatorial field with total area of $\rm 3.6\,deg^{2}$, within which there is available near-infrared photometry from the 4th Data Release (DR4) of the VIDEO survey produced by the Vista Science Archive (VSA; Irwin et al., 2004; Hambly et al., 2008; Cross et al., 2012). We adopt this approach because near-infrared data are needed for the estimation of stellar masses and specific accretion rates for AGN host galaxies (see next section).

The compilation of spectroscopic redshifts for X-ray sources in the Chandra fields is presented in Georgakakis et al. (2015) and in the XMM-XXL in Georgakakis et al. (2017) and Menzel et al. (2016). These papers also describe the estimation of photometric redshifts for X-ray sources based on methods developed by Salvato et al. (2009, 2011). For the new Chandra COSMOS-Legacy survey, spectroscopic redshifts are from the public releases of the VIMOS/zCOSMOS bright project (Lilly et al., 2009), the compilation of spectroscopy for X-ray sources presented by Civano et al. (2012) and Marchesi et al. (2016a), the MOSFIRE Deep Evolution Field (MOSDEF, Kriek et al., 2015), the PRIsm MUlti-object Survey (PRIMUS, Coil et al., 2011), the VIMOS Ultra-Deep Survey (VUDS, Le Fèvre et al., 2015), the 3D-HST programme (Skelton et al., 2014; Momcheva et al., 2015) and the Sloan Digital Sky Survey (SDSS; Alam et al., 2015). Photometric redshifts for the Chandra COSMOS-Legacy X-ray sources are from Marchesi et al. (2016a). For the small number of X-ray sources in our reduction of the Chandra COSMOS-Legacy data that do not appear in the Marchesi et al. (2016a) catalogue photometric redshifts are estimated following the methods described in Aird et al. (2015).

Figure 1 presents the distribution on the luminosity-redshift plane of the X-ray sources detected in the Chandra and XMM-XXL survey fields described above. Table 1 provides further information on the number of X-ray sources in each field as well as the corresponding spectroscopic or photometric redshift completeness. Figure 2 presents the sensitivity curves of the Chandra and XMM-XXL surveys.

2.2 XMM Slew Survey

The data from the Chandra and XMM dedicated survey programmes are supplemented by the wide area and shallow XMM Slew Survey (XMMSL, Saxton et al., 2008). Spectroscopic redshifts for the XMMSL sources are primarily from the SPIDERS (SPectroscopic IDentification of eROSITA Sources, Dwelly et al., 2017) programme of the SDSS-IV (Smee et al., 2013; Dawson et al., 2016) using the SDSS telescope (Gunn et al., 2006) and spectrographs (Smee et al., 2013). Our analysis uses X-ray sources in the “clean” version 1.6 of the XXMSL111http://www.cosmos.esa.int/web/xmm-newton/xmmsl1d-ug catalogue that overlap with a total of 3159 SDSS spectroscopic plates observed up to 11 May 2016 (inlcuding 392 observed as part of the SDSS-IV). Spectroscopic data obtained up to that date will be included in the SDSS Data Release 14 (DR14), which will be made public in summer 2017. The total overlap between the XMMSL v1.6 and the DR14 plates is $\rm 1165\,deg^{2}$. Within that area there are 942 detections which correspond to 816 unique XMMSL sources. Figure 3 shows the sky distribution of the DR14 spectroscopic plates within the SDSS-IV footprint and the spatial distribution of the XMMSL exposures within these plates.

The XMMSL catalogue lists sources detected independently in three energy bands, $0.2-2$, $0.2-12$ and $2-12$ keV. In this paper we use $0.2-2$ keV detected sources from the ‘clean’ version of the XMMSL catalogue. The choice of energy interval is motivated by the higher sensitivity and lower expected spurious detection fraction (0.7% Saxton et al., 2008) compared to e.g., the 2-12 keV selected sample. The soft-band selection introduces a bias against obscured sources. The XMMSL sample however, is dominated by luminous AGN ($L_{X}(\rm 2-10\,keV)\ga 10^{44}\,erg\,s^{-1}$l; see Fig. 1), which are expected to include a large, if not dominant, fraction of unobscured sources (e.g., Akylas et al., 2006; Ueda et al., 2014).

Because of the small exposure time per pixel (typically few seconds) the XMMSL sources typically consist of a small number of photons detected against an almost zero background. The estimated count rates and fluxes are therefore severely affected by the Eddington bias. We account for this effect using the source and background counts as well as the exposure times listed in the XMMSL source catalogue. Under the assumption of Poisson statistics the probability of a source having a observed count rate $CR$ is

[TABLE]

where $T$ is the observed number of counts (source and background) and $E$ is the expected number of counts assuming that the source count rate on the EPIC-PN detector is $CR$222The EPIC-MOS detectors are not used in slew mode., the exposure time is $t$, and the level of background at the source position is $B$

[TABLE]

The conversion of 0.2-2 keV count rates to fluxes in the same energy interval assumes a power-law spectral energy distribution with index $\Gamma=1.9$ modified by the Galactic neutral hydrogen column density, $N_{H}$, in the direction of each source333This is different from the XMMSL v1.6 catalogue, where the listed fluxes are determined assuming a fixed Galactic hydrogen column density $N_{H}=\rm 3.0\times 10^{20}cm^{-2}$.. For the 0.2-2 keV XMMSL sample that overlaps with the SDSS DR14 plates the mean Galactic neutral hydrogen column density is $\log[N_{H}/\rm cm^{-2}]=20.3$ and the standard deviation of the distribution is 0.26.

Although sensitivity maps are not among the data products of the XMMSL, it is possible to generate the selection function post-processing (e.g., Warwick et al., 2012). We construct the 0.2-2 keV sensitivity curve under the assumption that at the flux limit of the XMMSL ($f_{X}(\rm 0.2-2\,keV)\ga 5.7\times 10^{-13}\,erg\,s^{-1}\,cm^{-2}$; Saxton et al., 2008) the differential number count distribution of the extragalactic component of the XMMSL catalogue is well represented by a power-law with the Euclidean slope, i.e., $dN/dCR\propto CR^{-2.5}$. This is a reasonable assumption given that the X-ray number counts in the 0.5-2 keV band only deviate from the Euclidean slope below the flux level of $f_{X}(\rm 0.5-2\,keV)\approx 10^{-14}\,erg\,s^{-1}\,cm^{-2}$ (e.g., Georgakakis et al., 2008). The sensitivity curve of the sample can then be determined as the ratio of the observed extragalactic differential number counts, $(dN/dCR)_{obs}$, i.e., before accounting for the source detection limit selection effects, and the expected true differential number count distribution, i.e., $dN/dCR\propto CR^{-2.5}$. The $(dN/dCR)_{obs}$ is constructed by summing up the $P(CR,T)$ probability density functions of individual sources estimated from Equation 1. Figure 2 presents the resulting sensitivity curve.

The XMMSL sources are associated with mid-infrared counterparts in the AllWISE catalogue (Cutri et al., 2013) using a search radius of 60 arcsec. This angle should be compared with the typical $1\,\sigma$ rms positional uncertainty of the XMMSL sources of 8 arcsec (Saxton et al., 2008). Full details of the targeting of the XMMSL sources are presented in Dwelly et al. (2017). In brief, a novel Bayesian cross-matching algorithm is used (Salvato et al. in prep.), which is based on the techniques introduced by Budavári & Szalay (2008) and Rots & Budavári (2011). In addition to astrometric uncertainties, the identification is guided by prior information on the expected distribution of X-ray sources in the mid-infrared colour-magnitude space defined by the WISE $[W2]$ magnitude and the $[W1-W2]$ colour. The prior is constructed using bright X-ray sources [$f_{X}(\rm 0.2-2\,keV)\ga 10^{-13}\,erg\,s^{-1}\,cm^{-2}$] from the 3XMM survey (Watson et al., 2009; Rosen et al., 2016). The output of the Bayesian cross-matching algorithm is the posterior probability, $P_{post}$, of an XMMSL source being associated with an AllWISE counterpart. In the case of multiple potential associations the one with the highest $P_{post}$ is identified as the counterpart. XMMSL associations with $P_{post}>0.01$ are among the candidates for follow-up spectroscopy as part of the SPIDERS programme of the SDSS-IV (Dwelly et al., 2017). At this $P_{post}$ threshold the expected spurious identification rate is about 12%. This value is derived empirically using a sample of bright X-ray sources from the 3XMM Watson et al. (2009), the Swift XRT Point Source Catalogue (1SXPS; Evans et al., 2014) and the Chandra Source Catalog v1.1 (CSC; Evans et al., 2010). Optical counterparts are searched for in the SDSS photometric catalogue using the AllWISE positions and a search radius of 1.5 arcsec. The spectroscopic redshifts for the XMMSL sample are from the SDSS-IV SPIDERS (Dwelly et al., 2017) and eBOSS (extended Baryon Oscillation Spectroscopic Survey, Myers et al., 2015) programmes, as well as previous SDSS spectroscopic surveys, i.e., SDSS-III BOSS (Baryon Oscillation Spectroscopic Survey, Dawson et al., 2013) and the SDSS-I (Abazajian et al., 2009).

The extragalactic sample used in this work is selected from a total of 595 XMMSL sources that are detected in the 0.2-2 keV spectral band (B5 in XMMSL conventions) above the maximum likelihood threshold det_ml=10 and overlap with the SDSS-DR14 plates. As a first step we exclude sources that lie within the bright-star masks defined as part of the tiling process of the SDSS-III BOSS (Dawson et al., 2013) and SDSS-IV eBOSS (Dawson et al., 2016). This filtering reduces the number of 0.2-2 keV selected XMMSL sources to 479.

Next we apply a cut in the Bayesian identification posterior probability, $P_{post}>0.01$, i.e., the same threshold used for defining the targets for follow-up spectroscopy with SPIDERS. This criterion reduces the XMMSL sample to 440 sources with AllWISE counterparts. The 39 sources with $P_{post}<0.01$ include AllWISE associations without optical photometric counterparts in the SDSS (total of 16), stars (total of 3), galaxy cluster candidates (total of 2). These 39 sources are excluded from the analysis.

Among the 440 sources with $P_{post}>0.01$ there is a large fraction of non-AGN. Stars are identified either as bright and saturated point sources on the SDSS images (33), as objects with significant non-zero proper motions (total of 13; Munn et al., 2004), or as sources with stellar optical spectra observed as part of the SDSS spectroscopic surveys (9). X-ray sources that lie within 1 arcmin of a catalogued galaxy cluster or group are also marked as non-AGN (18). We also flag a total of 9 X-ray sources that are associated with BL Lacertae type objects (BL Lac) based on spectroscopic information from either the SDSS or the literature. The spectral energy distribution of these radio-loud sources is dominated by emission from a relativistic jet, believed to be aligned to the line-of-sight of the observer. The X-ray properties of these sources cannot be used as a proxy of the intrinsic accretion luminosity of the active nucleus. The above classes of sources (candidate stars, candidate clusters/groups and BL Lacs; total of 82) are excluded from the analysis. This selection process reduced the 0.2-2 keV selected XMMSL sample to 358 sources.

A total of 337 of this sample have redshift information from either the SDSS spectroscopic surveys up to DR14 (301) or the literature (36; NASA Extragalactic Database). The majority of these redshifts are spectroscopic (333) although a small number (4) of photometric redshift estimates (Csabai et al., 2007) are also included. XMMSL sources without spectroscopic redshift measurements (total of 21) include failed redshifts (total of 3), optically bright X-ray sources ($r\la 17$ mag; total of 7), which are not targeted to reduce cross-talk between neighbouring fibres, sources with no SDSS photometric counterparts (total of 6), which may include spurious associations, and optically fainter X-ray sources that have not been targeted because of fibre-collision problems, the faint magnitude limit of SPIDERS target selection ($r=22.5$) or because they are scheduled to be observed after 11 May 2016 and are therefore not part of DR14.

In summary, there are potentially 397 0.2-2 keV selected XMMSL X-ray sources likely to be associated with AGN (358 + 39 sources; the latter being those with $P_{post}<0.01$). We use in the analysis a total 337, for which redshift information is available. These objects represent 85% of the parent sample of candidate AGN. The 15% incompleteness also includes residual Galactic star contamination, X-ray galaxy clusters, as well as spurious X-ray detections and associations. Figure 1 presents the distribution on the X-ray luminosity vs redshift space of all X-ray sources used in this paper, i.e., from the Chandra, XMM-XXL and XMMSL surveys. This figure demonstrates the complementarity of the different samples.

2.3 Stellar mass estimates for AGN host galaxies

Stellar masses are estimated using the cigale code version 0.8 (Code Investigating GALaxy Emission, Noll et al., 2009) to fit templates to the observed multiwavelength photometry (UV to mid-infrared) of X-ray selected AGN. cigale convolves input star-formation histories with Single Stellar Population (SSP) models to synthesise the Spectral Energy Distribution (SED) of galaxies. Attenuation of stellar radiation by dust and the re-emission of the absorbed energy in the infrared are also accounted for in a self consistent manner. The emission associated with AGN is modelled using the Fritz et al. (2006) library of template SEDs as described in Ciesla et al. (2015). These templates correspond to the transmitted radiation of a central source through obscuring material with toroidal geometry. Free parameters of the model include the ratio between the outer and inner radius of the torus, the opening angle of torus, the optical depth and density profile of the obscuring material and the viewing angle of the observer relative to the equatorial axis. The choice of model parameters for the AGN templates follow Ciesla et al. (2015). Table 2 presents the assumptions and range of cigale parameters adopted in this work to generate galaxy and AGN/galaxy composite templates. We use the delayed star-formation history model, in which the star-formation rate (SFR) as a function of time, $t$, is parametrised as $SFR(t)\propto t\times\exp^{-t/\tau_{1}}$, where $\tau_{1}$ is the e-folding time of the young stellar population. This parametric form provides a more realistic description of the time-averaged star-formation rate of galaxies in semi-analytic simulations (Ciesla et al., 2015). We also adopt the Bruzual & Charlot (2003) SSPs with solar metallicity, the Chabrier (2003) Initial Mass Function (IMF) and the Calzetti et al. (2000) dust extinction law. The UV/optical stellar emission absorbed by dust is remitted in the IR assuming the Dale et al. (2014) templates.

The input photometry to cigale is different for each field. In the case of the CDFS we use the photometric catalogue of Hsu et al. (2014) that combines photometry from the the far-ultraviolet ($\rm\approx 1500$Å) to the mid-infrared ($\rm 8.0\mu m$). Only the broad-band filters listed in Table 1 of Hsu et al. (2014) are used to determine stellar masses. For the AEGIS-XD the multiwavelength ($\rm\approx 1500$Å  to $\rm 8.0\mu m$) photometric catalogue presented by Nandra et al. (2015, see their Table 7) was used in the calculation of stellar masses. In the COSMOS-Legacy field the fluxes in the CFHT $u^{\star}$, SUBARU $Vg^{+}r^{+}i^{+}z^{+}$ (Capak et al., 2007), UKIRT WFCAM $J$ (McCracken et al., 2010), CFHT WIRCAM $Ks$ (Capak et al., 2007) and the four IRAC (3.6, 4.5, 5.8 and $8.0\,\mu m$) filters were used. For X-ray sources detected in the XMM-XXL survey we used optical photometry ($ugriz$) from the CFHTLenS (Canada-France Hawaii Telescope Lensing Survey; Heymans et al., 2012; Erben et al., 2013), near-infrared fluxes ($YJHK$) from the VISTA (Visible and Infrared Survey Telescope for Astronomy) Deep Extragalactic Observations survey (VIDEO; Jarvis et al., 2013) based on the 4th Data Release (DR4) of the Vista Science Archive (VSA), and the mid-infrared photometry in the WISE (Wide-field Infrared Survey Explorer mission; Wright et al., 2010) W1, W2, W3, and W4 bands with central wavelengths of 3.4, 4.6, 12, and $\rm 22\,\mu m$, respectively from the AllWISE data release. The photometry for the XMMSL Survey sample is from the GALEX (Galaxy Evolution Explorer) All-Sky Imaging Survey Data Release 5 (GALEX-GR5; Bianchi et al., 2011), the Sloan Digitial Sky Survey DR13 (SDSS Collaboration et al., 2016), the Two Micron All Sky Survey point source catalogue (Skrutskie et al., 2006) and the AllWISE catalogue.

cigale provides probability distribution functions for the output parameters determined from the SED fits. For each model template the $\chi^{2}$ is determined based on the input photometry of each galaxy. The probability distribution function of the derived parameters of interest (e.g., stellar mass, star-formation rate) for each source is then estimated based on the $\chi^{2}$ values of the full set of model SEDs. Uncertainties associated with e.g., photometric errors or degeneracies among different templates are therefore accounted for, at least to the first approximation, and can be propagated into the analysis. Figure 4 shows the distribution of the AGN sample on the stellar mass vs redshift parameter space.

A fraction of the X-ray source population has photometric redshift estimates (40%, see Table 1) with uncertainties described by the corresponding probability distribution functions (PDZs). When estimating stellar masses we account for these errors. The PDZ of each source is first sampled in steps of $\Delta z=0.1$ between redshifts $z=0$ and $z=6$. For each bin $i$ at redshift $z_{i}$ the cigale code is applied to the observed photometry with redshift fixed to $z_{i}$ to yield a stellar mass probability distribution function. The normalisation of the stellar mass probability distribution function is not unity but the value of the PDZ at the corresponding redshift bin $i$. This approach produces 2-dimensional probability distribution maps in redshift and stellar mass, which provide a measure of the uncertainty in the determination of both quantities. We use these maps in the subsequent analysis to determine the specific accretion-rate distribution of X-ray selected AGN. Figure 5 shows an example of such a 2-dimensional probability distribution map.

2.4 Black hole mass estimates for broad line QSOs

The XMM-XXL and XMMSL survey fields, because of their shallower depth compared to the Chandra surveys used in this work, include a large fraction of type-1 AGN with broad optical-emission lines and blue continua. The determination of stellar masses for the host galaxies of these sources is challenging because the AGN continuum and emission lines contaminate or even dominate the observed broad-band fluxes. Stellar emission of the underlying host galaxy may however, be significant at longer wavelengths (i.e., rest-frame near-infrared, Bongiorno et al., 2012) and therefore the cigale decomposition of the observed Spectral Energy Distribution into galaxy and AGN contributions can at least partially mitigate AGN contamination issues. Nevertheless, some level of systematic uncertainties should be expected and it is therefore desirable to investigate their impact on the results and conclusions. For broad-line AGN it is also possible to have an independent estimate of the mass scale of the system via the determination of single-epoch virial black-hole masses (e.g., Shen & Liu, 2012). This calculation comes with its own set of systematics and random errors, which however are likely to be different from those associated with template fits to the observed broad-band Spectral Energy Distribution (e.g., Ciesla et al., 2015). Virial black-hole mass estimates therefore allow us to explore the sensitivity of our results and conclusions to the difficulties in determining masses for type-1 sources. Black-hole masses for the XMM-XXL and the XMMSL samples are estimated for broad-line AGN with SDSS spectroscopy based on methods similar to those described in Shen et al. (2011) and Shen & Liu (2012). For sources in the XMM-XXL the black hole masses are presented in Liu et al. (2016). For the XMMSL AGN they are estimated by one of us (AS).

Figure 6 displays stellar mass estimated via SED fits against black-hole mass determined from single-epoch optical spectroscopy for broad-line sources in the XMM-XXL and XMMSL surveys. Also shown is the local relation between (bulge) stellar-mass and black-hole mass $M_{BH}=0.002M$ (e.g., Marconi & Hunt, 2003; Häring & Rix, 2004). Broad line AGN are systematically offset by a factor of about three toward larger black-hole masses compared to the local relation. The direction and amplitude of this offset has been reported in previous studies (e.g Lauer et al., 2007; Matsuoka et al., 2014). It may indicate redshift evolution of the $M_{BH}-M$ relation (e.g Treu et al., 2004; Peng et al., 2006; Canalizo et al., 2012), differences in the calibration and normalisation of this relation compared to local quiescent black-holes (e.g Matsuoka et al., 2014) or the impact of selection biases (Lauer et al., 2007; Salviander et al., 2007; Schulze & Wisotzki, 2011). Although the redshift evolution of the $M_{BH}-M$ relation is debated (e.g., Shields et al., 2003; Shen et al., 2008; Merloni et al., 2010), it is generally accepted that selection effects play an important role. For example, In Figure 6 we plot the relation between the mean stellar mass of galaxies for a given black-hole mass estimated by convolving the local $M_{BH}-M$ relation with the stellar mass function of galaxies (Schulze & Wisotzki, 2011). The deviation from the local relation at the high-mass end is a consequence of the exponentially decreasing space density of massive galaxies and the scatter of the $M_{BH}-M$ relation (Lauer et al., 2007; Schulze & Wisotzki, 2011). At fixed black-hole mass there are more galaxies that host overmassive black holes for their stellar masses. Any additional sample selection effects, e.g. flux limits, luminosity cuts, will further impact the distribution of AGN in Figure 6. It also interesting that in the recent work of Reines & Volonteri (2015) the oposite trend from Figure 6 is observed, i.e. their local sample ($z<0.055$) of low-luminosity broad-line AGN selected from the SDSS has lower black hole masses for their stellar mass compared to local relations. This difference is also likely related to sample selection effects.

It is not the goal of this study to explore the $M_{BH}-M$ relation of active galaxies, but rather to investigate systematics in the inferred specific accretion-rate distribution of AGN by testing independent methods of measuring the mass scale of the system. In the analysis that follows, the stellar masses estimated by cigale are replaced by the virial black-hole mass estimates to explore systematic differences in the results. Because our analysis requires a measure of the stellar mass of the AGN host (see section 3), virial black-hole masses are converted to stellar masses using the relation $M_{BH}=0.002M$ (Marconi & Hunt, 2003). We mitigate issues related to the intrinsic scatter of this scaling relation or possible evolution with redshift by assigning a Gaussian uncertainty of $\sigma=0.5$ dex (e.g. Shen, 2013) to the logarithmic stellar mass of each source determined by scaling the virial black-hole masses.

In the analysis that follows, broad-line AGN in the XMM-XXL (total of 206 sources) and XMMSL (total of 253 sources) surveys are assigned stellar mass probability distribution functions that are determined indirectly by scaling the single-epoch virial black hole mass estimates. We have tested that our main results and the conclusions do not change if instead we use stellar masses determined from template fits to the broad-band spectral energy distribution of individual sources. The comparison of the specific accretion-rate distributions estimated using the two different approaches to approximate the stellar masses of broad-line AGN is presented in the Appendix. Small differences are observed at the high specific-accretion rate end of the distribution in the redshift interval $z=0.0-0.5$ (see Fig. 21). This level of difference indicates that current uncertainties in the determination of stellar masses for QSOs (e.g. factor of about 3 on the average in Fig. 6) is adequate for our purposes. We also emphasise the relatively small number of broad-line quasars within our X-ray selected sample, which also plays a role in the level of difference between the specific accretion-rate distributions compared in Figure 21.

3 Methodology

3.1 The specific accretion-rate distribution

The determination of the specific accretion-rate distribution of AGN is based on a modified version of the standard Poisson likelihood used for the calculation of the X-ray luminosity function (Aird et al., 2010, 2015; Buchner et al., 2015; Georgakakis et al., 2015)

[TABLE]

where ${\rm d}V/{\rm d}z$ is the comoving volume per unit solid angle and unit redshift at redshift $z$, $d_{i}$ signifies the data available for source $i$ (e.g. X-ray counts, redshift and stellar mass estimates) and the summation is over all sources, $N$, in the sample. The quantity $\phi_{\rm T}(L_{\mathrm{X}},z\;|\;\mathbf{\theta_{\rm T}})$ is the luminosity function model with parameters $\theta_{\rm T}$ that describes the space density as a function of redshift and X-ray luminosity of all extragalactic X-ray populations. The latter includes mostly AGN but also a small but non-negligible contribution from normal galaxies with X-ray emission dominated by stellar processes and/or hot interstellar gas (Hornschemeier et al., 2002; Georgantopoulos et al., 2005; Georgakakis et al., 2007; Lehmer et al., 2016). The total luminosity function of all extragalactic populations can therefore be expressed as the sum of the AGN and normal galaxy luminosity functions

[TABLE]

where $\phi(L_{\mathrm{X}},z\;|\;\mathbf{\theta})$ is the AGN luminosity function with parameters $\theta$ and $\phi_{\rm gal}(L_{\mathrm{X}},z\;|\;\mathbf{\theta_{\rm gal}})$ is the normal galaxy luminosity function described by a different set of parameters $\theta_{\rm gal}$. In equation 3 the multiplication is over all sources, $N$, in the sample, and the integration is over redshift and X-ray luminosity. The quantity $p(L_{\mathrm{X}},z|\;d_{i})$ is the probability of a particular source having redshift $z$ and X-ray luminosity $L_{\rm X}$ given the observational data, $d_{i}$. This approach captures uncertainties in the determination of both redshifts (e.g., photometric redshifts) and X-ray fluxes because of Poisson statistics and the Eddington bias. In equation 3 the parameter $\mu$ is the expected number of detected sources in the survey area for a particular set of luminosity function model parameters $\theta_{T}$

[TABLE]

where, $A(L_{\mathrm{X}},z)$ is the X-ray selection function that measures the probability of detecting a source with X-ray luminosity $L_{X}$ and redshift $z$ within the survey area.

The equations above can be modified by introducing the Eddington ratio distribution of AGN, $P(\lambda,z)$, which is defined as the probability of a galaxy at redshift $z$ and with stellar mass $M$ hosting an AGN accreting at an Eddington ratio $\lambda$. Given this definition the AGN luminosity function can be written

[TABLE]

where $\psi(M,z)$ is the stellar mass function of galaxies at redshift $z$. The Eddington ratio $\lambda$ is a function of the X-ray luminosity and stellar mass, under the assumption that the black hole mass is related to the stellar mass of the AGN host galaxy and that the X-ray luminosity can be used as proxy of the bolometric AGN luminosity. The quantity $P(\lambda,z)$ as defined above is a probability density function, i.e., integrates to unity at any guiven redshift

[TABLE]

Equation 14 is used to express the AGN luminosity function $\phi(L_{X},z)$ in Equations 3 and 5 in terms of the galaxy mass function and the Eddington ratio distribution of AGN

[TABLE]

where the likelihood is explicitly split into AGN and normal galaxy components. It is also assumed that each source in the survey has an equal probability of being a galaxy or an AGN, i.e., no prior knowledge such as morphology or X-ray–to–optical flux ratio is used to inform the nature of individual sources (see also Aird et al., 2015). In the equation above $\omega$ represents the set of parameter related to the Eddington ratio probability distribution function model that are to be determined from the observations. The quantity $p(L_{\mathrm{X}},z,M|\;d_{i})$ is the probability of a particular source having redshift $z$, X-ray luminosity $L_{\rm X}$ and hosted by a galaxy with stellar mass $M$. The expected number of detected sources in a survey for a particular set of model parameters $\omega$

[TABLE]

In the equations above it is assumed that the stellar mass function of galaxies is known to a reasonable level of accuracy. We adopt the parametrisation presented by Ilbert et al. (2013) for the total galaxy mass function. They used two Schechter functions (Press & Schechter, 1974) with parameters evolving with redshift to represent the mass function of galaxies in the redshift interval $z=0-4$. Using other galaxy mass functions available in the literature (e.g., Muzzin et al., 2013) does not change our main results.

Assumptions also must be made to estimate the Eddington ratio of AGN at redshift $z$, with X-ray luminosity $L_{X}$, hosted by galaxies with stellar mass $M$. Ideally direct estimates of black hole masses are needed for the determination of Eddington ratios. This is a challenging task however, particularly in the absence of broad emission lines in the optical spectra of large fraction of X-ray selected AGN because of e.g., obscuration or AGN light dilution by stellar emission. In this paper a Magorrian-type relation is adopted Magorrian et al. (1998) to scale the stellar mass of the galaxy to black hole mass. There are lengthy discussions in the literature on the relation between galaxy bulge mass proxies and black hole mass, the level of scatter of such correlations (e.g., Kormendy & Ho, 2013) or possible evolution with redshift (e.g., Treu et al., 2004; Peng et al., 2006; Canalizo et al., 2012). For the shake of simplicity, the baseline results presented in this paper adopt a single and redshift-independent scaling relation between black hole mass and bulge mass, $M_{BH}=0.002M_{bulge}$ (Marconi & Hunt, 2003). We further assume no scatter in that relation and that the AGN hosts are bulge dominated, i.e., that the bulge mass can be approximated by the stellar mass of the AGN host galaxy, $M=M_{bulge}$. A single X-ray bolometric conversion factor is further adopted, $L_{bol}=25\,L_{X}(\rm 2-10\,keV)$ (Elvis et al., 1994). Under these assumptions the Eddington ratio of AGN is estimated as

[TABLE]

where the X-ray luminosity is in units of $\rm erg\,s^{-1}$ and the stellar mass is in solar units. The $\lambda$ parameter is essentially a scaled version of the specific accretion rate defined as $L_{X}/M$ (e.g., Bongiorno et al., 2012; Aird et al., 2012, 2013; Bongiorno et al., 2016) expressed in Eddington-ratio units. From here on we refer to $\lambda$ as the specific accretion rate and $P(\lambda,z)$ as the specific accretion-rate distribution of AGN, to emphasize that these quantities are derived using the stellar mass as a proxy of black-hole mass for the bulk of the X-ray selected sample. We also underline that the only purpose of Equation 10 is to express the specific accretion rate in Eddington ratio units, which can be interpreted in a physical context. Our results and conclusions are insensitive to the assumptions adopted to derive Equation 10, e.g. bolometric correction. Nevertheless we have also explored the impact on the results and conclusions of different definitions of the $\lambda$ quantity. We adopt for example, the luminosity-dependent bolometric corrections of Marconi et al. (2004), instead of the fixed conversion factor of Equation 10. We also introduce scatter of $\sigma=0.38$ in the relation between the black hole mass and bulge mass of Marconi & Hunt (2003). Although these different definitions of $\lambda$ yield $P(\lambda,z)$ distributions that differ in detail from our baseline model assumptions, the main characteristics of these distributions as well as our results and conclusions on the redshift evolution of the AGN duty-cycle (see Results and Discussion sections) remain unchanged.

The equations above ignore the impact of obscuration on the estimation of the luminosity of individual AGN and the determination on the specific accretion-rate distribution. A statistically robust treatment of the obscuration requires the estimation of neutral hydrogen column densities, $N_{H}$, from the X-ray spectra of individual AGN and the inclusion of model terms that describe the dependence of the specific accretion-rate distribution of AGN on $N_{H}$, i.e., $P(\lambda(L_{X},M),z,N_{H})$. We defer this multi-dimensional analysis to future work. Accounting for obscuration effects will increase the space density of AGN we recover as by-product from our analysis, particularly at intermediate X-ray luminosities. The higher normalisation of the recovered X-ray luminosity function translates to an increase of the AGN duty cycle, i.e. the probability of a galaxy to host an AGN, and hence a higher amplitude for the the inferred $P(\lambda,z)$ at a given redshift. The impact of obscuration effects on our results will also decrease with increasing redshift because of the nature of the X-ray k-corrections.

Poisson statistics are used to determine the flux distribution that is consistent with the extracted source and background counts. A power-law X-ray spectrum with $\Gamma=1.9$ is adopted in this calculation. The resulting flux distribution is combined with the redshift information, spectroscopic or photometric redshift PDZ, to estimate the luminosity distribution of each source at rest-frame energies 2-10 keV. The relevant k-corrections also assume a power-law X-ray spectrum with $\Gamma=1.9$. Equation 10 is used to determine the $\lambda$ distribution of individual sources given their probability density distributions in X-ray luminosity stellar mass, and redshift (see also section 2.3).

A non-parametric approach for the determination of $P(\lambda(L_{X},M),z)$ is adopted. Our baseline model assumes that the specific accretion-rate distribution of AGN is independent of the stellar mass of the host galaxy (e.g., Aird et al., 2012, 2013). A two-dimensional grid in redshift and specific accretion rate is defined and $P(\lambda(L_{X},M),z)$ is assumed to be a constant within each grid pixel with dimensions ($\log\lambda\pm\mathrm{d}\log\lambda$, $z\pm\mathrm{d}z$). The value of the AGN specific accretion-rate distribution in each grid pixel is determined via equation 8. The edges of the grid pixels in each of the two dimensions are $\log\lambda=$ ($-10$, $-6$, $-5.5$, $-5.0$, $-4.5$, $-4.0$, $-3.5$, $-3.0$, $-2.5$, $-2.0$, $-1.5$, $-1.0$, $-0.5$, $0.0$, $0.5$, $1.0$, $1.5$, $2.0$) and $z=$(0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0). In Section 4 we will also discuss modifications to the baseline model to explore the stellar-mass dependence of the $P(\lambda(L_{X},M),z)$.

Importance sampling (Press et al., 1992) is used to evaluate the integrals in equation 8. For each source we draw random ($\log\lambda$, $z$, $M$) points from the full three-dimensional probability distribution. The specific accretion-rate distribution model is then evaluated for each sample point. The integral of equation 8 is simply the average $P(\lambda(L_{X},M),z)$ of the sample. The Hamiltonian Markov Chain Monte Carlo code Stan444http://mc-stan.org/shop/ (Carpenter et al., 2016) is used for Bayesian statistical inference.

The Appendix compares our results on the specific accretion-rate distribution of AGN with those of Aird et al. (in prep.). The later work developed parallel to our study and uses an independent methodology to infer the $P(\lambda,M,z)$ and a different set of multiwavelength photometric catalogues to estimate galaxy properties, such as stellar mass. The agreement between the two distinct approaches compared in Figure 22 is remarkable.

3.2 The X-ray luminosity function of normal galaxies

The normal galaxy (non-AGN) X-ray luminosity function parameters can in principle be determined simultaneously with the specific accretion-rate distribution parameters. For the sake of speed however, we follow the alternative approach of adopting a parametric model for the X-ray luminosity function of galaxies with parameters fixed to values from the literature. This approach is motivated by the fact that the X-ray luminosity function and evolution with redshift of this class of sources have been studied to some detail in the past (e.g., Georgantopoulos et al., 2005; Georgakakis et al., 2006; Ptak et al., 2007; Tzanavaris & Georgantopoulos, 2008; Aird et al., 2015). These previous studies can therefore inform our expectations for the contribution of normal galaxies to the current sample. Moreover, normal galaxies represent only a small component of current X-ray survey source populations (e.g., Lehmer et al., 2016). We adopt the parametrisation of Aird et al. (2015) and use a Schecter form for the X-ray luminosity function of galaxies

[TABLE]

where the characteristic luminosity, $L_{X}^{\star}$, evolves as

[TABLE]

The parameters are fixed to the maximum-likelihood values determined by Aird et al. (2015, their Table 10). For completeness they are also listed in Table 3. Equations 8 and 9 are then used to determine the specific accretion-rate distribution of AGN by keeping the normal galaxy luminosity function parameters fixed. This approach ignores uncertainties in the determination of the galaxy X-ray luminosity function. The advantage is that the required CPU time for the determination of the $P(\lambda,z)$ is significantly reduced.

3.3 The non-parametric luminosity function of X-ray selected AGN

A by-product of our analysis is the reconstructed AGN X-ray luminosity function derived from Equation 14 by convolving the $P(\lambda,z)$ distribution with the stellar mass function of galaxies. It is desirable to assess the quality of this reconstruction by comparing with independent estimates of the AGN space density as a function of redshift and X-ray luminosity.

For this exercise we derive non-parametric estimates of the X-ray luminosity function of AGN. A 2-dimensional grid in X-ray luminosity and redshift is defined with pixel edges $\log L_{X}(\rm 2-10\,keV)=$ ($38$, $39$, $40$, $41$, $42$, $43$, $43.5$, $44$, $44.5$, $45$, $46$), $z=$(0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0). Equation 8 is then used to determine the AGN space density, $\phi(L_{X},z)$, assuming that this quantity is constant within each grid pixel. The Hamiltonian Markov Chain Monte Carlo code Stan is used to sample the likelihood function described by Equation 3 and estimate the model parameters. We account for the contribution of normal galaxies to the X-ray source population by adopting the Schechter model described above for the luminosity function of these sources and fixing the parameters to the values listed in Table 3. Appendix C compares the reconstructed and non-parametric determinations of the AGN X-ray luminosity function.

4 Results

4.1 The specific accretion-rate distribution of AGN

As a first approximation we ignore a possible dependence of the specific accretion-rate distribution of AGN on the stellar mass of their host galaxies; such a dependence will be explored later in this section. Equation 8 is applied to the full sample of X-ray selected AGN without a division into stellar mass bins. The resulting constraints on the $P(\lambda,z)$ are displayed in Figure 7, where the different panels correspond to the redshift intervals defined in section 3. The width of the hatched regions in Figure 7 corresponds to the 90% confidence interval and shows the region of the parameter space where the data provide meaningful constraints to the specific accretion-rate distribution of AGN. As expected, the sample becomes insensitive to low specific accretion-rate AGN with increasing redshift. This behaviour is manifested by a threshold in the specific accretion-rate below which the width of the hatched regions is comparable to the full range of the $P(\lambda,z)$ prior (i.e., many orders of magnitude). This specific accretion-rate threshold moves to higher values with increasing redshift.

In each redshift panel of Figure 7 the $P(\lambda,z)$ distributions, i.e. the probability of a galaxy hosting an accretion event, increase with decreasing specific accretion rate, at least to the limit $\lambda\approx 10^{-4}-10^{-3}$. Also, the slopes of the distributions steepen close to the Eddington limit, $\lambda\approx 0$. These results are consistent with studies that approximate the specific accretion-rate distribution of AGN with a single power-law and find evidence for a break at approximately the Eddington limit (Aird et al., 2012; Bongiorno et al., 2012). This break can be parametrised by a second steeper power-law distribution (Aird et al., 2013). At very low specific accretion rates, $\lambda\la 10^{-4}-10^{-3}$ there is evidence for a flattening of the $P(\lambda,z)$ distribution in all panels of Figure 7. This trend is driven by the data and is only partially imposed by the requirement that the $P(\lambda,z)$ distribution integrates to unity.

Figure 7 also compares our results with the parametric specific accretion-rate distribution of Aird et al. (2012). They model the $P(\lambda,z)$ distribution with a single power-law that has a normalisation that evolves with redshift as $\propto(1+z)^{\gamma}$. These results are formally valid for AGN with specific accretion-rate in the range $\lambda=10^{-4}-1$ and for redshifts $z\la 1$. We choose not to extrapolate the Aird et al. (2012) results beyond $z>1.5$ and outside the specific accretion-rate interval $\lambda=10^{-4}-1$. Our analysis suggests a $P(\lambda,z)$ shape that is more complex than a single power-law. Nevertheless the functional form derived by Aird et al. (2012) roughly approximates the slope and overall redshift evolution of the low Eddington-ratio end of the distributions plotted in Figure 7. Also shown in this figure are the power-law fits to the specific accretion-rate distributions estimated by Bongiorno et al. (2012) for their samples at redshift intervals $z=0.3-0.8$, $0.8-1.5$ and $1.5-2.5$. The specific accretion rates in the Bongiorno et al. work are converted to the Eddington-ratio units adopted in our paper using equation 10. At the lowest redshift bin there is fair agreement between the slope found by Bongiorno et al. (2012) and our results. At higher redshift their slope appears to be in better agreement with the high specific accretion-rate end of our results. The specific accretion-rate distributions estimated by Bongiorno et al. (2016) are also plotted in Figure 7, after scaling their specific accretion rate to the units adopted in our work using Equation 10. We plot the specific accretion-rate term of Bongiorno et al. (2016) (their equation 11) estimated at the average mass $\log M/M_{\odot}=10.5$ of our sample. Because of the different methodology adopted by Bongiorno et al. (2016) their specific accretion rate distribution has no normalisation. Their curves in Figure 7 have been arbitrarily normalised to roughly match the amplitude of the $P(\lambda,z)$ distributions derived in our work. Nevertheless, the overall shape of the distributions, particularly at $z\ga 1$, is similar. We also caution that the Bongiorno et al. (2016) results are an extrapolation at redshifts below $z\approx 0.5$.

The redshift dependence of the $P(\lambda,z)$ distribution can be explored by comparing the different panels of Figure 7. It can be argued that the fundamental characteristics of the $P(\lambda,z)$ do not change drastically with redshift. Despite small differences in the slopes, which will be discussed later in this section, the general shape of the $P(\lambda,z)$ remains roughly the same across redshift. The most striking change is that of the normalisation of the $P(\lambda,z)$, which increases toward higher redshift. This trend appears to saturate at $z>2$, i.e., in the redshift bins $z=2.0-2.5$, $z=2.5-3.0$ and $z=3.0-4.0$. This behaviour is further demonstrated in Figure 8, which compares the $P(\lambda,z)$ distribution derived in the redshift interval $z=0.0-0.5$ to the distributions determined at higher redshifts. In each panel of this figure a vertical shift has been applied to the $P(\lambda,z=0.0-0.5)$ distribution to facilitate the comparison. The amplitude of this shift is indicated in each panel of Figure 8. We emphasise that these shifts are empirically visually estimated to demonstrate that the observed redshift evolution of the specific accretion-rate distribution of AGN can be represented in a rough manner by a change of the normalisation of the $P(\lambda,z=0.0-0.5)$.

An implication of the evolutionary scheme in which the fundamental shape of the $P(\lambda,z)$ is roughly redshift independent is that the typical specific-accretion rate of the AGN population does not evolve strongly with redshift. We show this point using the the quantity $\lambda\,P(\lambda,z)$, in which the linear power-law slope of the specific accretion-rate distribution is removed to highlight the value of $\lambda$ that dominates the AGN population at a given redshift. This is plotted in Figure 9. For comparison in that figure the distribution derived in the redshift interval $z=0.0-0.5$ is overplotted on the distributions determined at higher redshifts. In each panel of this figure the same vertical shift as in Figure 8 has been applied to the $\lambda\,P(\lambda,z=0.0-0.5)$ distribution to facilitate the comparison. In Figure 9 the $\lambda\,P(\lambda,z)$ shows a plateau or a peak at $\log\lambda\approx-1$ that does not change significantly with increasing redshift. Another way to demonstrate this point is to estimate the mean specific accretion-rate of the distributions in Figure 7, i.e., the integral of the quantity $\lambda\,P(\lambda,z)$

[TABLE]

where the $\lambda_{l}$ is the lower limit of the integration. For this exercise the integration is performed in the $\lambda$ interval over which the $P(\lambda,z)$ distributions are well constrained. We choose as an integration lower limit $\lambda_{l}=10^{-3}$ and restrict this calculation to redshift $z<2.5$. At higher redshift our $P(\lambda,z)$ constraint suffers large uncertainties for $\lambda<10^{-3}$. Different choices of the integration lower limit will yield different mean specific accretion rates but will not alter the redshift dependence of this quantity, which is the focus of this exercise. This is presented in Figure 10 to redshift $z=2.5$. For comparison also shown in this figure is the expected amplitude of the evolution of the specific star-formation rate of galaxies in the same redshift interval (Elbaz et al., 2011). This comparison is interesting because the star-formation rate density and the accretion density of the Universe follow similar patterns (e.g., Aird et al., 2010, 2015), and links are often proposed between the formation of stars in galaxies and the growth of black holes at their centres (e.g., Santini et al., 2012; Rosario et al., 2012). Figure 10 demonstrates than the mean specific accretion rate of the X-ray AGN population evolves mildly with redshift when compared with the level of evolution of the specific star-formation rate of galaxies in the same redshift range. At relatively low redshift, $z\la 1$ the mean specific accretion rate is nearly constant, while at higher redshifts it appears to track the overall evolution of the specific star-formation rate of galaxies. This can be interpreted as evidence that accretion events onto supermassive black holes decouple from the formation of new stars in galaxies below redshift $z\approx 1$.

Next, the zero-order approximation of a universal shape for the specific accretion-rate distribution is scrutinised by investigating changes with redshift of the basic shape characteristics of the $P(\lambda,z)$. This is easier done using the quantity $\lambda\,P(\lambda,z)$ plotted in Figure 9, which has the linear power-law slope of the specific accretion-rate distributions removed to highlight subtle differences in slopes among redshift bins. In Figure 9 there is evidence that at high redshift, $z\ga 1$, the slope of the specific-accretion rate distributions below $\log\lambda\approx-2$ is flatter than that at low redshifts, $z=0.0-0.5$. This trends is manifested by the offset between the red (lowest redshift bin; $\lambda\,P(\lambda,z=0.0-0.5)$) and the blue (high redshift bin) distributions at $\log\lambda\la-2$. This difference becomes more pronounced at $z>1.5$ and suggests that the value of the specific accretion rate below which the distributions of Figure 7 turn over depends on redshift, i.e., shifts to higher $\lambda$ with increasing redshift. Similar trends are claimed by Bongiorno et al. (2016) and Aird et al. (in prep.).

Figure 11 explores subtle differences in the redshift evolution of the probability of a galaxy experiencing accretion events of a given specific accretion rate. The left panel of that figure plots the $P(\lambda,z)$ in different $\lambda$-bins as a function of redshift. The right panel faciliates the comparison between different $\lambda$-bins by normalising them to a fixed value at low redshift. Low and moderate specific accretion-rate bins ($\log\lambda\la-0.5$) evolve at a similarly high rate between redshifts $z=0$ and $z\approx 1$. At higher redshift however, $z\approx 1-3$, the evolution rate of low specific accretion-rate bins ($\log\lambda\la-1.5$) slows compared to higher specific accretion-rate intervals ($-1.5\la\log\lambda\la 0$). This effect is manifested for example, in the right panel of Figure 11 by the fact that the (normalised) probability of an accretion event with $\log\lambda=-0.75\pm 0.25$ lies above that for $\log\lambda=-1.25\pm 0.25$ or $\log\lambda=-1.75\pm 0.25$ in the redshift interval $z\approx 1-3$. It is also interesting that close to the Eddington limit, i.e., for $\lambda=-0.25\pm 0.25$ and $\lambda=+0.25\pm 0.25$, the evolution rate of the $P(\lambda,z)$ starts off slow at low redshift, contrary to other $\lambda$-bins, but then increases and broadly follows the evolution of the high specific accretion-rate intervals. This type of evolution implies that the probability (or duty-cycle) of moderate and high specific accretion-rate events increases relative to that of low specific accretion-rate episodes toward higher redshift. Therefore one expects an increasing fraction of moderate and high specific accretion-rate AGN with increasing redshift. It therefore appears that the evolution, in which only the overall normalisation of the $P(\lambda,z)$ distributions changes with redshift, needs to be revised to include a $\lambda$-dependent differential evolution. We emphasise neverthless, that this a subtle effect, with the main evolution pattern being that of an overall increase in normalisation.

We further explore the stellar-mass dependence of the $P(\lambda,z)$ distribution by splitting the sample into stellar mass bins and determining the specific accretion-rate distribution for each separately. The results are shown in Figure 12, which compares the $P(\lambda,z)$ distributions in three mass bins: $\log M/M_{\odot}=10.0-10.5$, $\log M/M_{\odot}=10.5-11.0$ and $\log M/M_{\odot}=11.0-11.5$. There is a systematic trend whereby more massive galaxies are less likely to host AGN with high specific accretion rates approaching the Eddington limit, $\log\lambda\approx 0$. This trend is more evident at redshifts $z\la 2$, where the uncertainties of the $P(\lambda,z)$ distributions in individual mass bins are smaller. At the low specific accretion-rate end of the $P(\lambda,z)$ distributions ($\log\lambda\la-1$) there is fair agreement between the different mass samples, at least for the redshift panels of Figure 12, where the uncertainties of each mass sample are reasonably small.

In the analysis that follows we present results using the specific accretion-rate distributions that do not depend on the stellar mass. This because of the relatively large uncertainties of the $P(\lambda,z)$ constraints for individual mass bins. Nevertheless, where necessary we will also discuss the impact on the results and conclusions of the mass-dependent specific accretion-rate distributions plotted in Figure 12.

4.2 X-ray AGN mass function and duty cycle

The redshift evolution of the $P(\lambda,z)$ distributions of Figure 7, which can be approximated by an increase of the normalisation with increasing redshift, indicates that galaxies are more likely to host an AGN at earlier cosmic times.

This result can be demonstrated by estimating the stellar mass function of AGN hosts above a given X-ray luminosity cut following relation

[TABLE]

where the integration is over the range of specific accretion rates that yield X-ray luminosities (for a given stellar mass) larger than the cut $L_{cut}$. Figure 13 presents the stellar mass function of AGN at different X-ray luminosity cuts in comparison with the mass function of the galaxy population determined by Ilbert et al. (2013). Also displayed in Figure 13 are the stellar mass functions of X-ray selected AGN with $L_{X}(\rm 2-10\,keV)>10^{41}\,erg\,s^{-1}$ derived by Georgakakis et al. (2011) at mean redshifts $z\approx 0.1$ and $z\approx 0.8$. Although that study shares some of the data used in the present analysis, the mass functions have been estimated using a very different methodology compared to the one presented here. The reasonably good agreement between the mass functions of Georgakakis et al. (2011) and those derived in this paper further underlines the robustness of the results. A feature of Figure 13 is that with decreasing X-ray luminosity cut the mass function of AGN becomes typically broader and extends to lower stellar masses. Another important point to be taken from Figure 13 is that at fixed stellar mass and X-ray luminosity cut the fraction of AGN among galaxies increases toward higher redshift and even approaches unity for certain regions of the stellar mass, redshift and X-ray luminosity threshold parameter space.

Another way of presenting this latter result is via the AGN duty cycle, i.e., the ratio of the mass functions of AGN and galaxies above a particular X-ray luminosity cut. This relation is presented as a function of stellar mass and for different X-ray luminosity cuts in Figure 14. For comparison, independent duty-cycle observational constraints from the literature (Shi et al., 2008; Haggard et al., 2010; Georgakakis et al., 2011) are also shown at the appropriate redshift-interval panels of Figure 14. At fixed stellar mass and X-ray luminosity threshold the duty cycle of AGN increases with increasing redshift.

4.3 X-ray AGN downsizing

AGN downsizing is defined as the trend whereby the space density of luminous AGN peaks at earlier cosmic times compared to systems with lower accretion luminosities (Barger et al., 2005; Hasinger et al., 2005). This differential evolution with X-ray luminosity is shown in Figure 15, which plots as a function of redshift our non-parametric constraints on the space density of AGN in different X-ray luminosity bins. In this section we discuss this observational trend in the context of the specific accretion rate distributions derived in our work.

Figure 16 displays in different redshift intervals the space density of AGN as a function of stellar mass (mass function) and specific accretion rate ($\lambda$ function) in the same luminosity intervals as in Figure 15. These quantities are estimated by convolving the mass function of galaxies with the $P(\lambda,z)$ distribution and integrating over the appropriate X-ray luminosity range. The mass functions in the top set of panels of Figure 16 become broader with decreasing X-ray luminosity and the peaks of the distributions shift to somewhat lower stellar mass with increasing redshift. The latter point is further demonstrated in the left panel of Figure 17, which presents the median stellar mass of the mass functions of Figure 16. The decreasing trend of the median stellar mass with increasing redshift is more pronounced for lower luminosity AGN. Nevertheless, the variation in stellar mass is about 0.5 dex or less over a wide range of redshifts, even for the lowest luminosity interval $\log L_{X}(\rm 2-10\,keV)=43.0-43.5$. A similar amplitude of redshift variation is therefore expected for the specific accretion rate of the X-ray luminosity subsamples. The specific accretion rate function is shown in the bottom set of panels of Figure 16. There is a clear trend whereby lower luminosity AGN are systematically shifted to lower specific accretion rates, albeit with a tail extending to high $\lambda$ values. The median specific accretion rate of these distributions is displayed in the middle panel of Figure 17. There is a clear offset in the median specific accretion rate of the subsamples selected by X-ray luminosity. AGN with luminosities $\log L_{X}(\rm 2-10\,keV)=44.5-45.0$ are typically associated with black holes accreting with $\log\lambda\approx-0.7$, while those with luminosities $\log L_{X}(\rm 2-10\,keV)=43.0-43.5$ typically have accretion rates $\log\lambda\approx-1.7$, with some level of variation across redshift. The trend whereby X-ray selected AGN split by luminosity have distinct median specific accretion rates is also evident when using the mass-dependent $P(\lambda,z)$ distributions to convolve with the mass function of galaxies and produce the curves in Figure 17. This result is therefore robust to the stellar mass dependence of the specific accretion-rate distribution of AGN.

The AGN downsizing of Figure 15 can therefore be understood as the result of two effects. The first is the stratification in median specific accretion rate of AGN samples selected by X-ray luminosity, as shown in Figure 17. The second effect is that high specific accretion-rate events are more likely relative to low specific-accretion ones with increasing redshift. In other words the probability of a galaxy hosting an AGN with a certain value of the specific accretion rate depends on both $z$ and $\lambda$. These trends are indicated in Figure 11. They are further demonstrated in the right panel of Figure 17, which plots the duty-cycle of AGN in the luminosity intervals defined in Figure 15. The duty-cycle is the probability of galaxy with a given stellar mass hosting an AGN with X-ray luminosity in a given range. This is estimated by dividing the the AGN stellar mass functions plotted in the top set of panels of Figure 16 with the total galaxy stellar mass function and then integrating along the stellar mass axis. AGN with $\log L_{X}(\rm 2-10\,keV)=43.0-43.5$ include a large fraction of sources with $\lambda\la-1.5$ and therefore their duty-cycle evolution slows at redshifts $z\ga 0.5-1$ relative to AGN with $\log L_{X}(\rm 2-10\,keV)=44.5-45.0$. The latter population includes a large fraction of accretion events with $\lambda\ga-1.0$, which are characterised by a higher rate of duty-cycle evolution to redshifts $z\approx 2-2.5$. The duty-cycle evolution of all accretion events eventually slows at redshift beyond $z\ga 2.5$ resulting in the observed decrease of the space density of AGN at early cosmic times.

5 Discussion

In this paper we combine wide-area/shallow and pencil-beam/deep X-ray surveys to compile one of the largest samples of X-ray selected AGN to date, with a total of 4821 sources. The X-ray data are supplemented by multiwavelength observations to derive the specific accretion-rate distribution of AGN, $P(\lambda,z)$. This quantity is defined as the probability of a galaxy at a given redshift, $z$, hosting an AGN with specific accretion rate $\lambda$. Our approach uses the mass function of galaxies as a boundary condition and has the advantage that is non-parametric in nature, i.e., no assumptions are made on the underlying shape or model that describes the $P(\lambda,z)$ distribution. Additionally, uncertainties in the determination of the properties of individual AGN in our sample, such as host galaxy stellar masses, X-ray luminosities and redshifts are propagated into the analysis.

5.1 AGN evolution and

relation to star-formation

Our non-parametric method recovers a specific accretion-rate distribution of AGN (see Fig. 7) that increases toward low specific accretion rates to $\log\lambda\approx-3$ and shows a break close to the Eddington limit, above which the probability of a galaxy hosting an accretion event drops steeply. Below $\log\lambda\approx-3$ there is evidence that the increasing trend of the $P(\lambda,z)$ is inverting and the distribution drops. There is also evidence that the turnover point at low specific accretion rates, $\log\lambda\la-3$, is a function of redshift, i.e., it shifts to higher $\lambda$ values with increasing redshift. These characteristics are broadly consistent with previous observational studies that either derive the specific accretion-rate distribution of AGN via non-parametric methods (Aird et al., 2012; Bongiorno et al., 2012), or adopt an analytic model for the specific accretion rate distribution and recover the relevant parameters by requiring that the convolution of this model with the AGN or galaxy mass functions yields the AGN luminosity function (Aird et al., 2013; Bongiorno et al., 2016). Recent hydrodynamical cosmological simulations of galaxy formation and black-hole growth also produce Eddington ratio distributions that are qualitatively similar, in terms of shape and redshift evolution, to the observational results shown in Fig. 7 (Sijacki et al., 2015).

We find that the specific accretion-rate distribution of AGN evolves strongly with cosmic time so that the probability of a galaxy hosting an AGN increases toward higher redshifts. In other words, galaxies were on average more nuclear-active in the past. This trend translates to a higher AGN duty cycle at earlier epochs and for more massive galaxies (see Fig. 14). At $z>1.5$ and $\rm M>10^{11}\,M_{\odot}$, for example, nearly 100% of the galaxies are predicted to host an AGN with $L_{X}(\rm 2-10\,keV)>10^{41}\,erg\,s^{-1}$. It is interesting to compare the high duty-cycles inferred by our analysis with independent observational estimates of the incidence of AGN among galaxies at high redshift. Genzel et al. (2014) studied the kinematics of massive galaxies $\rm\log M/M_{\odot}>10.9$ at $z=1-3$ and found that 2/3 of their sample show broad nuclear emission-line components $\approx\rm 450-5300\,km\,s^{-1}$. These objects may be associated with black-hole accretion events and suggest a high duty-cycle of AGN outflows at the redshifts and stellar-mass intervals above. Such a high fraction is not inconsistent with the duty cycles implied by our analysis.

The rapid increase of the AGN duty-cycle from redshift $z=0$ to $z\approx 1-1.5$ is also driving the strong evolution of the AGN space and luminosity densities between these redshifts (see also Aird et al., 2013). In Figure 14 the change in amplitude of the AGN duty-cycle at fixed stellar mass and luminosity cut from low redshift to $z\approx 1$ is about 1 dex. This is similar to the increase of the corresponding integrated AGN space or luminosity densities (e.g. Ueda et al., 2014; Aird et al., 2015).

The observed increase of the AGN duty-cycle toward higher redshift may be related to the evolution of the star-formation rate in galaxies. As a first step we explore such a relation in a qualitative manner. Star-forming galaxies are known to occupy a relatively tight sequence on the star-formation vs stellar mass plane (Noeske et al., 2007; Speagle et al., 2014; Whitaker et al., 2012; Tomczak et al., 2016). The dominant evolution pattern is a shift of this sequence toward higher star-formation rates with increasing redshift, i.e., a change of the sequence’s normalisation (e.g., Speagle et al., 2014). Put differently, at fixed stellar mass galaxies form stars at a higher rate with increasing redshift (e.g., Speagle et al., 2014; Ilbert et al., 2015). Figure 18 shows our observational results on the mean X-ray luminosity of galaxies as a function of stellar mass and for different redshift intervals. Our constraints are compared with X-ray stacking analysis results from the literature (Mullaney et al., 2012). At fixed stellar mass galaxies have systematically higher X-ray luminosities with increasing redshift, i.e., qualitatively similar behaviour as the star-formation rate. Similar trends are evident when using the mass-dependent specific accretion-rate distributions of Figure 12 to determine the mean X-ray luminosity as a function of stellar mass and redshift. In this case however, the slopes of the $L_{X}-M$ relations at $z\la 2$ are shallower than unity, i.e. the slope derived using the mass-independent $P(\lambda,z)$ distribution. The difference is more striking at the lowest redshift-bin, $z=0.0-0.5$ and becomes less pronounced with increasing redshift.

Next we explore quantitatively the relation between the growth of black holes and the star-formation rate in galaxies. Figure 19 compares the accretion rate onto the black hole, $dm/dt$, with the star-formation rate (SFR) of Main Sequence galaxies at the corresponding redshift. For the conversion of the mean X-ray luminosities of Fig. 12 to black-hole accretion rates we assume radiative efficiency of 10% and a fixed bolometric correction of 25. The relation between SFR, stellar mass and redshift is from Aird et al. (2017). Figure 19 shows that the overall normalisations of the Main Sequence of star-formation and mean black-hole accretion rate evolve similarly with redshift. However, the mass dependence of the Main Sequence SFR and $dm/dt$ have different slopes. Adopting a stellar-mass independent $P(\lambda,z)$ distribution to derive the mean X-ray luminosity of galaxies yields a slope of unity, compared to slopes in the range $0.6-1$ for the star-formation rate (e.g. Speagle et al., 2014). Taking into account the mass-dependence of the $P(\lambda,z)$ introduces a tilt into the $dm/dt$ vs stellar mass relation, the amplitude of which is a function of redshift. At $z=0.0-0.5$ the mass-dependence of the black-hole accretion rate is flatter than that of the star-formation rate, while at $z=1.5-2$ it is steeper. This suggests that the relation between star-formation and black-hole accretion events is redshift dependent.

A different manifestation of this effect is presented in Figure 10, which plots the redshift dependence of the mean specific accretion rate, $\overline{\lambda}$, of the $P(\lambda,z)$ distribution. We find that this quantity remains nearly constant out to redshift $z\approx 1$, and then at higher redshifts follows the evolution of the specific star-formation rate of Main Sequence galaxies. There is therefore evidence for two distinct evolutionary regimes of the mean specific accretion rate of AGN. At high redshift the quantity $\overline{\lambda}$ couples to the Main Sequence of star-formation, while at lower redshift it appears to separate from the evolution of the specific star-formation rate of the galaxy population. This differential behaviour with redshift can be interpreted as evidence for multiple black-hole accretion mechanisms that dominate at different cosmic epochs. This scenario is consistent with some recent cosmological Semi-Analytic Models (Fanidakis et al., 2012; Hirschmann et al., 2012), which invoke different black-hole fuelling modes to reproduce the evolution of AGN and galaxies. In these models the different accretion modes take place in distinct environments with respect to the level of star-formation and become important at different redshifts.

Finally, we note that the different stellar-mass dependence of black-hole accretion rate and star-formation rate contradicts recent results by Mullaney et al. (2012). These authors estimated the mean X-ray and infrared luminosities of galaxies and find that these quantities have remarkably similar dependence on galaxy stellar mass and redshift. This was interpreted in the context of a universal ratio between mean black hole growth rate and SFR. The discrepancy between Figure 19 and the Mullaney et al. (2012) results may be related to sample selection effects (e.g. X-ray AGN vs infrared galaxies) or biases associated with the stacking analysis (e.g., Eddington bias, see Mullaney et al. 2012). It is also noted that the trends observed in Figure 19 are not sensitive to the adopted conversion from X-ray to bolometric luminosity. The use of the Marconi et al. (2004) bolometric corrections yields similar results.

5.2 Constraints on black-hole fuelling modes

Our results suggest that the probability of galaxy hosting an accretion event onto its supermassive black hole depends on both the redshift and the specific accretion-rate of the event. Figure 11 shows that high specific accretion-rate events ($\lambda\ga-1.5$) become more likely relative to low specific accretion-rate ones with increasing redshift. Bongiorno et al. (2016) parametrise the specific accretion-rate distribution of AGN by a double power-law with a slope at the low-$\lambda$ end that depends on redshift. They then infer a behaviour similar to that found by our non-parametric approach, i.e., a low-$\lambda$ slope that flattens with increasing redshift resulting in a differential evolution of the specific accretion-rate distribution. Merloni & Heinz (2008) solved the continuity equation of the black-hole mass function using the locally determined one as boundary condition and assuming that the evolution of the X-ray luminosity function of AGN traces the growth of black-holes as a function of cosmic time. They recover an Eddington-ratio distribution for AGN that flattens at higher redshift, i.e., indicating an increasing relative importance of highly accreting objects at earlier cosmic epochs. We emphasise that this differential evolution is a second order effect in our analysis. The bulk effect is an overall increase of the normalisation of the $P(\lambda,z)$ toward higher redshift.

The differential evolution of the $P(\lambda,z)$ with respect to $\lambda$ may be related to the fact that galaxies at higher redshift are, on average, more gas-rich compared to lower redshift counterparts having the same stellar mass (e.g., Genzel et al., 2010; Tacconi et al., 2010, 2013; Magdis et al., 2012; Santini et al., 2014). This description however, cannot be the only effect. If a larger reservoir of gas is available to the supermassive black hole, and the same black-hole fuelling mode or AGN triggering mechanism is in operation across all redshifts, then one might expect a higher frequency of accretion events with increasing redshift over the full range of specific accretion rates. Our analysis suggests the existence of additional physical processes that suppress lower specific accretion-rate events relative to high specific accretion-rate ones at higher redshift. One possibility is that there are multiple AGN fuelling or triggering modes that dominate at different cosmic times and specific accretion-rate regimes (see also Section 5.1). For example, one may imagine a simple scenario, in which the growth of black holes in gas-rich galaxies proceeds at relatively high Eddington rates, while in gas-poor early-type galaxies black holes grow at typically slower rates relative to the Eddington limit (e.g., Kauffmann & Heckman, 2009; Fanidakis et al., 2012). The latter fuelling mode may become less important with increasing redshift, as galaxies become, on average, more gas-rich at fixed stellar mass, thereby leading to the observed differential evolution of the $P(\lambda,z)$ distribution with respect to the specific accretion-rate. This scenario is consistent with the distribution of the AGN space density on the 2-dimensional parameter space of stellar mass and X-ray luminosity plotted in Figure 20. At low redshift the AGN population is dominated by low specific accretion-rate systems, $\log\lambda\la-2$. At higher redshift moderate specific accretion-rate events, $-2\la\log\lambda\la-1$, become increasingly prominent toward $z\approx 2$. At higher redshift the AGN population is dominated by accretion events with $-2\la\log\lambda\la-1$. This also partly because of the lack of observational constraints below $log\lambda\la-2$ at $z>2$.

Variations of the specific accretion-rate distribution of AGN as a function of the star-formation rate of their host galaxies, a proxy to gas availability, is an important diagnostic of multiple black-hole fuelling modes (Kauffmann & Heckman, 2009; Georgakakis et al., 2014; Azadi et al., 2015). There is evidence that high specific accretion-rate events are less common among quiscent galaxies compared to star-forming ones (Aird et al. in prep). Moreover, recent work suggests that the signatures of multiple black-hole fuelling modes are imprinted on the large scale environment of AGN and are manifested as systematic variations of the clustering amplitude with redshift and accretion luminosity (Fanidakis et al., 2013; Koutoulidis et al., 2013; Mountrichas et al., 2016).

Semi Analytic Models (SAM) that simulate the formation of galaxies and the growth of black holes at their centres must also adopt different AGN triggering processes to match the observational results on the cosmological evolution and demographics of AGN (e.g., Fanidakis et al., 2012; Hirschmann et al., 2012, 2014; Gutcke et al., 2015). AGN triggering mechanisms in these simulations include galaxy mergers, disk instabilities and hot-gas accretion, each of which becomes important at different cosmic times, dark matter halo environments and accretion luminosity regimes. For example, there are striking similarities between our Figure 20, which displays the space density of AGN in the 2-dimensional stellar mass vs X-ray luminosity space, and the Hirschmann et al. (2014) SAM predictions shown in their Figure 12, which presents the simulated AGN space density distribution in the black-hole mass and bolometric luminosity plane. In that figure of Hirschmann et al., there is a clear shift of the AGN population, in terms of space density, from low Eddington-ratio and high black-hole mass systems at low redshift (radio-mode accretion) to high-Eddington ratio and lower black-hole mass sources (quasar-mode accretion) with increasing redshift, i.e. similar to the trends in Figure 20. This behaviour of the simulated AGN in the Hirschmann et al. (2014) work is because of the interplay of the two black hole fuelling modes adopted in that SAM. The quasar-mode corresponds to high Eddington-ratio events, $\lambda_{Edd}>0.01$, while the radio-mode corresponds to $\lambda_{Edd}<0.01$. The difference between the two modes is the level of feedback, which is assumed to be more efficient in the case of the radio-mode black-hole fuelling. This assumption represents the impact on the interstellar medium of radio jets from massive and slow-accreting black holes hosted by early-type and massive galaxies. In that model the radio-mode component becomes significant below redshift $z\approx 1$, while at higher redshift the quasar-mode increasingly dominates the growth of black holes. Therefore in this simulation there is differential evolution of the Eddington ratio distribution, qualitatively similar to the evolution of $P(\lambda,z)$ we find in our work (e.g. see Fig. 11). High Eddington-ratio AGN evolve faster at moderate and high redshift compared to low Eddington-ratio ones.

Our analysis also suggests differences in the specific accretion-rate distribution of AGN split by stellar mass. There is evidence that more massive galaxies are less likely to host AGN with an accretion rate approaching the Eddington limit (Fig. 12). This trend is stronger at $z<2$ and disappears at higher redshift, although the uncertainties at $z>2$ are also larger. The observed variations in the shape of the $P(\lambda,z)$ with stellar mass may be related to gas availability in massive galaxies that also live, on average, in denser environments. Bongiorno et al. (2016) model the specific accretion rate distribution of AGN in the COSMOS survey using a double power-law and introduce a stellar mass dependence for the break point between the two power-law slopes. They find that the break moves to lower specific accretion rates with increasing stellar mass, similar to our non-parametric approach. Schulze et al. (2015) suggest that the type-I quasar duty-cycle decreases with increasing black-hole mass at least for redshift $z<1$, indicating a suppression of type-I AGN activity at higher mass systems, in qualitative agreement with our results.

5.3 Implications on the downsizing of X-ray selected AGN

Our analysis also provides observational clues on the origin of AGN downsizing. We find that this observational trend is primarily related to the redshift evolution of the AGN duty-cycle, i.e. the probability of a galaxy hosting an accretion event within a given X-ray luminosity interval.

Firstly, it is shown that the median stellar mass of AGN split by accretion luminosity differs by less than a factor of six for X-ray luminosities spanning two orders of magnitude between $L_{X}=10^{43}$ and $\rm 10^{45}\,erg\,s^{-1}$ (Fig. 17). There is also no strong redshift dependence of the median stellar mass between $z=0$ and $z=3$ for the luminosity selected sub-samples (maximum factor of three for AGN with $\log L_{X}=43-43.5$). These findings argue against black hole mass being primarily responsible for the AGN downsizing, under the assumption that stellar mass tracks black-hole mass. Moreover, our analysis suggests that AGN samples selected by luminosity differ in their median specific accretion rate, i.e., at least one order of magnitude between $L_{X}=10^{43}$ and $\rm 10^{45}\,erg\,s^{-1}$ (Fig. 17). More luminous sources are associated typically with higher specific accretion rates compared to their lower luminosity counterparts. There is also little redshift dependence of the median specific accretion rate for the luminosity selected AGN sub-samples. Accretion events within a given luminosity interval have the same median specific accretion rate across redshift. The parameter that is found to evolve differentially with redshift is the duty-cycle of AGN split by X-ray luminosity. The right panel of Figure 17 shows that the probability of a galaxy hosting a moderate X-ray luminosity event relative to a high X-ray luminosity one decreases with increasing redshift in the interval $z\approx 1-2.5$ (see also Fig. 11). We therefore find that the AGN downsizing trend is primarily a duty-cycle effect, i.e. it is associated with the differential evolution of the probability of a galaxy to host an accretion event. This conclusion does not change if we take into account the stellar-mass dependence of the specific accretion-rate distribution shown in Figure 12, whereby massive galaxies are less likely to host AGN with high accretion rates.

The conclusions above on AGN downsizing are compared to the recent work of Bongiorno et al. (2016). These authors argue that the main effect in AGN downisizing is the differential redshift evolution of the specific accretion-rate distribution of AGN, with the space density of high specific accretion-rate events peaking at earlier epochs. This is similar to our results, although in detail there are differences. We do not find evidence that the probability of galaxies hosting accretion events with different $\lambda$ values peaks at different redshifts (see Fig. 11). Instead, the evolution of the $P(\lambda,z)$ for different specific accretion-rate bins shows a broad plateau in the redshift interval $z\approx 1-3$. Also our conclusions on the origin of AGN downsizing differ from those of Schulze et al. (2015). They studied the Eddington-rate (not specific accretion rate) distribution of type-I quasars and found evidence that the black hole mass is the main driver of downsizing. This apparent discrepancy may indicate observational selection effects, i.e., broad-line QSOs vs X-ray detections, or the biases associated with the use of stellar mass to approximate the Eddington ratio of AGN.

Our results on the origin of the AGN downsizing are at least qualitatively consistent with recent simulations by Hirschmann et al. (2014). In this latter work it is the interplay between the quasar and radio mode accretion as well as the heavy black-hole seeding at early times that drives the downsizing trend of the AGN luminosity function. Therefore at least in that SAM the AGN downsizing is directly related to the presence of multiple modes for fuelling the black holes of galaxies.

Finally, our results can be used to explore the origin of the break $L_{X}^{*}$ of the X-ray luminosity function of AGN. Our non-parametric analysis suggests a break luminosity in the range $\log L_{X}^{*}/\rm erg\,s^{-1}\approx 43.5-44.5$ depending on redshift (see Fig. 23). Such AGN are associated with median specific-accretion rates $\lambda\approx 10^{-1.5}-10^{-1}$ (see Fig. 17). These range of values brackets the position where the slope of the specific accretion-rate distributions of Figure 7 steepens toward high $\lambda$ values. Additionally AGN with $\log L_{X}/\rm erg\,s^{-1}\approx 43.5-44.5$ are associated with median galaxy stellar masses in the range $\log M/M_{\odot}=10.5-11.0$. This interval brackets the break of the stellar mass function of the galaxy population $\log M/M_{\odot}\approx 10.7$ (Ilbert et al., 2013). The $L_{X}^{*}$ of the X-ray luminosity function is therefore related to the breaks of the stellar mass function and the specific accretion-rate distributions (Aird et al., 2013). It is also possible to show analytically that these quantities are related under the assumptions of a Schechter form for the stellar mass function and a double power-law for the specific accretion rate distribution (Caplar et al., 2015).

Conclusions

The main conclusions from our work are summarised below

• •

The non-parametric approach for the determination of the specific accretion-rate distribution of AGN recovers a shape that is characterised by an increase toward low specific accretion rates, a break close to the Eddington limit, above which the probability of a galaxy hosting an AGN drops steeply, and a turnover at very low specific accretion rates, $\log\lambda\la-3$.

• •

The redshift evolution of the specific accretion-rate distribution is such that the fraction of AGN among galaxies at a given accretion-luminosity limit (duty-cycle) increases with both increasing stellar mass and redshift.

• •

The rapid increase of the AGN duty-cycle with increasing redshift is primarily responsible for the strong evolution of the AGN space and luminosity densities between the local Universe and $z=1-1.5$.

• •

There is evidence for differential evolution of the AGN duty-cycle with X-ray luminosity or specific accretion rate. We argue that this trend is primarily responsible for the AGN downsizing.

• •

Our analysis suggests that the relation between accretion events onto supermassive black holes and star-formation episodes is redshift dependent. At redshift $z<1$ there is evidence that the evolution of the specific accretion rate of AGN is disjoint from the evolution of the specific star-formation rate of galaxies. At higher redshift the two quantities evolve similarly. This is consistent with suggestions that there are multiple black-hole fuelling modes that dominate at different redshifts.

• •

The specific accretion-rate distribution of AGN likely also depends on stellar-mass. More massive galaxies are systematically shifted to lower specific accretion rates compared to less massive systems.

Acknowledgements

The authors thank the anonymous referee for their careful reading of the paper and their constructive comments. This work benefited from the thales project 383549 that is jointly funded by the European Union and the Greek Government in the framework of the programme “Education and lifelong learning”. JA acknowledges support from ERC Advanced Grant FEEDBACK 340442. A.S. acknowledges support by JSPS KAKENHI Grant Number 26800098. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatory of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autonoma de Mexico, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Appendix A The impact of different approaches for assessing the stellar mass of AGN hosts

Two methods are adopted in our work to estimate the host galaxy stellar masses of AGN with broad optical-emission lines. The first fits hybrid AGN/galaxy templates to the observed broad-band spectral energy distribution of X-ray sources to decompose the emission of the underlying host from that of the nuclear source (Section 2.3). The second approach approximates the stellar mass of the galaxy by simply scaling the virial black-hole mass estimated from single-epoch spectroscopy (Section 2.4). The two methods are affected by different systematics and random errors and provide a handle to assess the impact of uncertainties in broad-line quasar stellar masses to the results and conclusions on the specific accretion-rate distribution.

Figure 21 compares the $P(\lambda,z)$ for the two different approaches of approximating the stellar mass of broad-line quasars. There is a steepening of the high specific accretion-rate end of the distributions, when using virial black-hole masses to approximate the stellar mass of the underlying galaxy. Most affected are the specific accretion-rate distribution in the two lowest redshift bins. This behaviour arises because these X-ray AGN populations include a large number of luminous broad-line quasars from the XMMSL sample. Nevertheless, the overall effect is moderate and has no impact on the any of our final results and conclusions.

Appendix B Comparison with Aird et al.

In a study that developed parallel to our work, Aird et al. (in prep.) investigated the specific accretion-rate distribution of X-ray selected AGN, using as starting point galaxy samples from the CANDELS-3DHST (GOODS-S, GOODS-N, AEGIS, COSMOS) and the COSMOS/UltraVISTA field. In the common fields between the Aird et al. (in prep.) and our studies the same X-ray observations are used. Aird et al. (in prep.) used the Bayesian mixture model described in Aird et al. (2017) to determine the specific accretion-rate distribution of the galaxy population in different stellar mass and redshift bins. The comparison between the specific accretion-rate distribution derived in our work and in Aird et al. (in prep.) is shown in Figure 22. Overall the agreement is good, considering the different multiwavelength data used in the two studies and the independent methodologies for inferring the specific accretion-rate distribution of AGN. At the high-$\lambda$ end of the distributions the Aird et al. (in prep.) results are systematically higher that ours, although we emphasise that the tension is not significant given the large uncertainties. It is likely that this difference arises because of the different methodologies, the smoothness priors adopted by Aird et al. (in prep.), and possibly the propagation of stellar mass and redshift errors in our work.

Appendix C The X-ray Luminosity Function of AGN

A by-product of our analysis deriving the $P(\lambda,z)$ distribution are new observational constraints on the X-ray luminosity function of AGN by combining Chandra and XMM surveys with different characteristics in terms of area and depth (see Table 1). This section presents these constraints using both a parametric and a non-parametric approach. It also compares direct estimates of the X-ray luminosity function with that reconstructed by convolving the $P(\lambda,z)$ distributions derived in Section 4 with the stellar mass function of galaxies using Equation 14.

The parametric X-ray luminosity function estimation is using the LADE model as described in Section 3. We also account for the contribution of non-AGN among X-ray sources using the Schechter parameterisation of the X-ray luminosity function of normal galaxies presented in Section 3. For the non-parametric estimate, a dimensional grid in X-ray luminosity and redshift is defined and the AGN space density is assumed to be constant within each grid pixel (see Section 3). We account for the contribution of normal galaxies as described in Section 3, i.e., by adopting a Schechter model for the luminosity function of these sources and fixing the parameters to those listed in Table 3.

Figure 23 compares the parametric and non-parametric constraints to the X-ray luminosity function of both AGN and galaxies. Also plotted in this Figure for comparison is the LADE model of Aird et al. (2015) for their soft-band selected sample. This figure shows the level of agreement between parametric and non-parametric luminosity function estimates, the contribution of the non-AGN population to the estimated space densities and the comparison with previous constraints from the literature.

Figure 24 compares the reconstructed AGN luminosity function based on the $P(\lambda,z)$ distributions derived in Section 4 with the non-parametric X-ray luminosity function derived after accounting for the contribution of normal galaxies. The agreement between the two estimates is reasonable.