The AGN XLF at z=3-5

The X-ray luminosity function of Active Galactic Nuclei in the redshift interval z=3-5

Abstract

We combine deep X-ray survey data from the Chandra observatory and the wide-area/shallow XMM-XXL field to estimate the AGN X-ray luminosity function in the redshift range . The sample consists of nearly 340 sources with either photometric (212) or spectroscopic (128) redshift in the above range. The combination of deep and shallow survey fields also provides a luminosity baseline of three orders of magnitude, at . We follow a Bayesian approach to determine the binned AGN space density and explore their evolution in a model-independent way. Our methodology properly accounts for Poisson errors in the determination of X-ray fluxes and uncertainties in photometric redshift estimates. We demonstrate that the latter is essential for unbiased measurement of space densities. We find that the AGN X-ray luminosity function evolves strongly between the redshift intervals and . There is also suggestive evidence that the amplitude of this evolution is luminosity dependent. The space density of AGN with drops by a factor of 5 between the redshift intervals above, while the evolution of brighter AGN appears to be milder. Comparison of our X-ray luminosity function with that of UV/optical selected QSOs at similar redshifts shows broad agreement at bright luminosities, . At fainter luminosities X-ray surveys measure higher AGN space densities. The faint-end slope of UV/optical luminosity functions however, is steeper than for X-ray selected AGN. This implies that the type-I AGN fraction increases with decreasing luminosity at , opposite to trends established at lower redshift. We also assess the significance of AGN in keeping the hydrogen ionised at high redshift. Our X-ray luminosity function yields ionising photon rate densities that are insufficient to keep the Universe ionised at redshift . A source of uncertainty in this calculation is the escape fraction of UV photons for X-ray selected AGN.

keywords:
galaxies: active – galaxies: Seyferts – X-rays: diffuse background – quasars: general

1 Introduction

In recent years observations have established that supermassive black holes (SMBHs) are nearly ubiquitous in local spheroids (Magorrian et al., 1998; Kormendy & Ho, 2013). These relic black holes are believed to have grown their masses at earlier times mostly via accretion of material from larger scales (e.g. Soltan, 1982; Marconi et al., 2004). Questions that remain open are when during the lifetime of the Universe these events occurred and under what physical conditions black holes grow their masses. Moreover, observations show that in the local Universe correlations exist between the mass of SMBHs and the properties of the stellar component of the bulges in which they reside, such as velocity dispersion (e.g Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Gültekin et al., 2009; Graham et al., 2011), luminosity (e.g McLure & Dunlop, 2002; Marconi & Hunt, 2003; Gültekin et al., 2009), dynamical mass (e.g. Magorrian et al., 1998; Marconi & Hunt, 2003; Häring & Rix, 2004; Graham, 2012) and central light concentration (Graham et al., 2001; Savorgnan et al., 2013). Such correlations suggest a link between the growth of black holes and the formation of their host galaxies, although the exact nature of such an interplay is still not well understood. Processes that can establish such correlations include a common gas reservoir that both feeds the central black hole and forms stars on larger scales, outflows related to the energy output from the AGN itself that affect the Inter-Stellar Medium and regulate the formation of stars (Silk & Rees, 1998; Fabian, 1999; King, 2003, 2005; Di Matteo et al., 2005; Croton et al., 2006), or the merging history of galaxies and their SMBHs (e.g. Jahnke & Macciò, 2011).

One approach for improving our understanding of the formation of SMBHs as a function of cosmic time and their relation to their host galaxies is population studies of Active Galactic Nuclei (AGN), which signpost accretion events onto SMBHs. This requires a census of the AGN population across redshift to constrain for example, the accretion history of the Universe or study the incidence of active black holes among galaxies. In that respect, the AGN luminosity function, i.e. their comoving space density as a function of redshift and accretion luminosity, is one of the fundamental quantities that characterise the demographics of active black holes. The cosmic evolution of AGN leaves imprints on the shape and overall normalisation of the luminosity function. The total mass density locked into black holes at different epochs can be inferred by direct integration of the AGN luminosity function, under assumptions about the radiative efficiency of the accretion process and after applying appropriate bolometric luminosity corrections (e.g Marconi et al., 2004; Aird et al., 2010; Ueda et al., 2014). The space density of AGN split by host galaxy properties, such as stellar mass, morphology or level of star-formation, provides clues on the interplay between black hole accretion and galaxy evolution (Georgakakis et al., 2009, 2011; Aird et al., 2012; Bongiorno et al., 2012; Georgakakis et al., 2014).

Selection at UV/optical wavelengths (e.g. Richards et al., 2009; Ross et al., 2012) currently provides the largest spectroscopic AGN samples for luminosity function calculations (e.g. Ross et al., 2013). The downside is that the UV/optical continuum of AGN is sensitive to dust extinction along the line–of–sight and dilution by the host galaxy at faint accretion luminosities. Observations at X-ray wavelengths can mitigate these issues (e.g Brandt & Alexander, 2015). X-ray photons, particularly at rest-frame energies , can penetrate nearly unaffected large columns of intervening gas and dust clouds (), thereby providing samples least affected by obscuration biases. Moreover, the X-ray emission associated with stellar processes is typically 2 orders of magnitude fainter than the AGN radiative output and therefore contamination or dilution effects by the host galaxy are negligible at X-rays over a wide baseline of accretion luminosities. X-ray surveys also benefit from a well defined selection function that is relatively easy to quantify and account for in the analysis. The disadvantage of X-ray selection is that the detected AGN are often optically faint and therefore spectroscopic follow-up observations are often expensive. Nevertheless, intensive multiwavelength campaigns in recent years substantially increased the number of X-ray survey fields with sufficient quality ancillary data for reliable X-ray source identification and redshift measurements using either spectroscopy or photometric methods. Early results by Cowie et al. (2003) on the redshift evolution of the X-ray AGN space density and Ueda et al. (2003) on the hard band (2-10 keV) X-ray luminosity function and obscuration distribution of AGN have been expanded recently both in terms of data and analysis methodology (Yencho et al., 2009; Ebrero et al., 2009; Aird et al., 2010; Burlon et al., 2011; Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015; Miyaji et al., 2015). Although important details on the shape of the X-ray luminosity function and the obscuration distribution of AGN are still debated (Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015; Miyaji et al., 2015), the overall evolution of the X-ray luminosity function is reasonably well constrained at least to . The initial increase of the AGN X-ray luminosity density from to is followed by a broad plateau up to and a decline at higher redshift. However, the amplitude of the X-ray AGN evolution at is still not well constrained. Early studies suggested a moderate decline of the AGN space density at (Yencho et al., 2009; Ebrero et al., 2009; Aird et al., 2010), contrary to claims for a rapid drop (Brusa et al., 2009; Civano et al., 2011; Vito et al., 2013; Kalfountzou et al., 2014) that can be parametrised by an exponential law in redshift (Gilli et al., 2007) similar to the optical QSO space density evolution (Schmidt et al., 1995; Richards et al., 2006). Central to this debate is the typically small X-ray AGN sample sizes at . For example, there are 209 and 141 X-ray AGN with spectroscopic or photometric redshifts above in the most recent compilations of Kalfountzou et al. (2014) and Vito et al. (2014) respectively. These numbers should be compared with sample sizes of few thousands AGN at for the most recent X-ray luminosity function studies (e.g Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015; Miyaji et al., 2015)

Better constraints on the form and amplitude of the evolution of X-ray AGN at will have implications for the contribution of this population to the UV photon field density that is needed to keep the Universe ionised at high redshift. Haardt & Madau (2012) used the Ueda et al. (2003) X-ray luminosity function to predict a moderate contribution of AGN to the hydrogen ionising radiation field at , in broad agreement with constraints derived from UV/optical selected QSO luminosity functions (e.g. Fontanot et al., 2007; Masters et al., 2012; McGreer et al., 2013) and other X-ray AGN studies (Barger et al., 2003a; Grissom et al., 2014). There are also claims however, that AGN provide an important contribution to the photoionisation rate at high redshift (Fiore et al., 2012; Glikman et al., 2011; Giallongo et al., 2015). This discrepancy emphasises the need for further work to improve measurements of the AGN space density at high redshift and to better understand their role in the re-ionisation of the Universe.

In this paper we combine deep Chandra and wide-area/shallow XMM-Newton survey fields to compile one of the largest samples of X-ray selected AGN at to date. A Bayesian methodology is developed to correctly account for photometric redshift uncertainties and to determine in a non-parametric way the AGN comoving space density in the redshift intervals and and over 3 decades in X-ray luminosity [ in ]. Although parametric models are also fit to the data, we emphasise the importance of non-parametric estimates to determine in a model-independent way the shape and overall evolution of the X-ray luminosity function. Throughout this paper we adopt , and .

2 Data

For the determination of the X-ray luminosity function in the redshift interval we combine Chandra and XMM-Newton X-ray surveys with different characteristics in terms of area coverage and X-ray depth. These are the 4 Ms Chandra Deep Field South (CDFS; Xue et al., 2011; Rangel et al., 2013), the 2 Ms Chandra Deep Field North (CDFN; Alexander et al., 2003; Rangel et al., 2013), the Extended Groth Strip International Survey field (AEGIS, Davis et al., 2007; Laird et al., 2009; Nandra et al., 2015), the Extended Chandra Deep Field South (ECDFS; Lehmer et al., 2005), the Chandra COSMOS field (C-COSMOS, Elvis et al., 2009) and the equatorial field of the XMM-XXL survey.

2.1 Chandra survey fields

The Chandra observations of the CDFS, CDFN, AEGIS, ECDFS and C-COSMOS were analysed in a homogeneous way by applying the reduction and source detection methodology described by Laird et al. (2009). Specific details on the analysis of the 4 Ms CDFS and the 2 Ms CDFN fields are presented by Rangel et al. (2013). The Chandra survey of the AEGIS field has two tiles. The wide and shallow one (AEGIS-W) consists of 8 Chandra pointings of 200 ks each. These data are described by Laird et al. (2009). The deep survey of the AEGIS field (AEGIS-XD) increased to a total of 800 ks the exposure time of the central regions of the AEGIS-W. The additional data overlap with the central 3 of the original 8 Chandra pointings observed as part of the AEGIS-W. The AEGIS-XD survey data reduction and source catalogue generation are described by Nandra et al. (2015).

The X-ray sources used in this paper are detected in the 0.5-7 keV (full) spectral band with Poisson false detection threshold . The count rates in the 0.5-7 keV band are converted to fluxes in the 0.5-10 keV band assuming a power-law spectral index with . This is steeper than the adopted for the X-ray flux estimation in the published catalogues of the CDFS, CDFN and AEGIS fields (Laird et al., 2009; Rangel et al., 2013; Nandra et al., 2015). The choice of the is motivated by the fact that at high redshift, , the observer-frame 0.5-7 keV band corresponds to harder rest-frame energies, which are least affected by obscuration. A diagnostic of the spectral shape of X-ray sources is their hardness ratio defined as HR=(H-S)/(H+S), where S, H are the observed count rates in the 0.5-2 and 2-7 keV spectral bands, respectively. For sources in the range we find a median hardness ratio . This is consistent with a power-law X-ray spectrum with , i.e. similar to the mean spectral index of the intrinsic AGN spectra (, Nandra & Pounds, 1994). We therefore choose to fix for the determination of fluxes. We note however, that the choice of (1.4 vs 1.9) has a small impact on the results presented in this paper.

Sensitivity curves, which measure the total survey area that is sensitive to sources of a particular flux are calculated following methods described in Georgakakis et al. (2008). The overlap between the ECDFS and the 4 Ms CDFS or the AEGIS-W and the AEGIS-XD is accounted for by defining independent spatial regions for each survey. Spatial masks that describe both the boundaries of the optical/infrared imaging of each field and regions of poor photometry because of bright stars (Aird et al., 2015) are also taken into account in the X-ray sensitivity calculations. The sensitivity curves in the 0.5-10 keV band are presented in Figure 1. Table 1 presents the number of X-ray sources in each field. The same spatial masks used for the construction of sensitivity maps are used to filter the X-ray source catalogue.

The optical identification of the X-ray sources in the CDFN, AEGIS-XD. AEGIS-W, ECDFS and C-COSMOS fields are based on the Likelihood Ratio method (Sutherland & Saunders, 1992) as implemented in Aird et al. (2015). The multiwavelength associations of the 4 Ms CDFS X-ray sources are presented by Hsu et al. (2014). They apply a Bayesian methodology, based on the work of Budavári & Szalay (2008), to different catalogues available in that field including the CANDELS/-band selected photometry presented by Guo et al. (2013), the Taiwan ECDFS Near-Infrared Survey (TENIS; Hsieh et al., 2012) and the MUSYC/-selected catalog of Cardamone et al. (2010). The number of X-ray sources with secure optical or infrared counterparts in each field are presented in Table 1.

Extensive spectroscopic campaigns have been carried out in the fields of choice. For the CDFN, ECDFS and AEGIS-W we use the compilation of spectroscopic redshifts presented by Aird et al. (2015). In the 4 Ms CDFS we use the spectroscopic redshifts compiled by Hsu et al. (2014). In the case of AEGIS-XD we use the spectroscopic redshift catalogue presented by (Nandra et al., 2015). Redshifts in the C-COSMOS are from the public releases of the VIMOS/zCOSMOS bright project (Lilly et al., 2009) and the Magellan/IMACS observation campaigns (Trump et al., 2009), as well as the compilation of redshifts for X-ray sources presented by Civano et al. (2012). The spectroscopic redshifts used in this paper have quality flags in the published catalogues from which they were retrieved that indicate a probability better than % of being correct.

For X-ray sources without spectroscopy, photometric redshifts are estimated using the multiwavelength photometric catalogues available for each survey field. The photometric redshifts of the X-ray sources in the 4 Ms CDFS, the AEGIS-XD and the COSMOS fields are determined following the methodology described by Salvato et al. (2009, 2011). Specific details can be found in Hsu et al. (2014; 4 Ms CDFS), Nandra et al. (2015; AEGIS-XD) and Salvato et al. (2011; COSMOS field). The estimated rms scatter of the X-ray AGN photometric redshifts is , 0.014 and 0.04 for the C-COSMOS, 4 Ms CDFS and AEGIS-XD samples, respectively. The corresponding outlier fraction, defined as , is about 5-6% in all three fields. In the case of ECDFS, AEGIS-W and CDFN we use the photometric redshifts estimated by Aird et al. (2015). For these fields and the outlier fraction is about 15%. The latter value is larger than in the C-COSMOS, 4 Ms CDFS and AEGIS-XD fields. This is related to differences in the methodology of estimating photometric redshifts and ultimately to the choice of template SEDs used in the calculation. Nevertheless, the photometric redshifts estimated by Aird et al. (2015) are assigned appropriately larger uncertainties, approximated by the corresponding Probability Distribution Functions (PDZ), that reflect the higher outlier fraction. Figure 2 plots spectroscopic vs photometric redshifts for the sample used in this paper and illustrates the overall quality of the photometric redshifts estimates. In the X-ray luminosity function calculations we use the full photometric redshift PDZ.These are typically unimodal but at increasing redshift they broaden and secondary peaks may also appear. Also the Aird et al. (2015) PDZs are typically broader than those in the C-COSMOS, 4 Ms CDFS and AEGIS-XD fields. Examples of PDZs used in this paper are shown in Figure 3. Table 1 presents the number of photometric and spectroscopic redshifts in each field, both total and in the intervals and . Sources without optical identifications and hence, without redshift estimates, are a minority in the Chandra surveys sample, 107 in total. These sources can be either moderate redshift () AGN with red SEDs because of e.g. obscuration and/or old stellar populations (Koekemoer et al., 2004; Schaerer et al., 2007; Rodighiero et al., 2007; Del Moro et al., 2009), or high redshift systems (). Since the redshift distribution of these sources is not known they are assigned a flat PDZ in the redshift interval and zero at other redshifts. Fixing the lower limit of the redshift range above to a value between and does not change the results and conclusions.

Additionally we have tested that the differences in the accuracy and outlier fraction of the photometric redshifts in the different fields used in this work do not affect the final results. For the COSMOS, AEGIS-XD and 4 Ms CDFS we can substitute the photometric redshift PDZs adopted in this paper (based on the methods of Salvato et al. 2009, 2011) with those estimated following the methodoloy of Aird et al. (2015). For these fields we can therefore estimate the X-ray luminosity function at (see next sections) using two different sets of photometric redshifts, those of Aird et al. (2015) and those determined following Salvato et al. (2009, 2011). We find no systematic differences between the X-ray luminosity functions at estimated using the two independent photometric redshift catalogues.

2.2 XMM-XXL survey data

The Chandra survey fields are complemented by the wide-area and shallow XMM-XXL survey (PI: Pierre). The XMM-Newton observations of this programme cover a total of about of sky split into two nearly equal area fields. In this paper we use the equatorial region of the XMM-XXL, which overlaps with the Canada-France-Hawaii Legacy Survey (CFHTLS) W1 field and covers an area of about . The approximate centre of this field lies at right ascension and declination .

The reduction of the XMM data, the construction of the X-ray catalogue and the association of the X-ray sources with optical counterparts are described by Liu et al. (2015) based on the methods presented in Georgakakis & Nandra (2011). In brief, the X-ray data reduction is carried out using the XMM Science Analysis System (SAS) version 12. We analyse XMM-Newton observations related to the XMM-XXL programme that were made public prior to 23 January 2012. XMM-XXL data observed after that date are not included in the analysis. As a result our final catalogue of the equatorial XMM-XXL field misses about worth of X-ray coverage. The epchain and emchain tasks of sas are employed to produce event files for the EPIC (European Photon Imaging Camera; Strüder et al., 2001; Turner et al., 2001) PN and MOS detectors respectively. Flaring periods resulting in elevated EPIC background are identified and excluded using a methodology similar to that described by Nandra et al. (2007). We use X-ray sources detected in the  keV spectral band with Poisson false detection probability of . The final sample consists of 7493 unique sources detected in 0.5-8 keV spectral band. The fluxes listed in the final source catalogue are in the 0.5-10 keV band assuming a power-law spectral energy distribution with . These X-ray sources are matched to the SDSS-DR8 photometric catalogue (Aihara et al., 2011) using the Maximum-Likelihood method (Sutherland & Saunders 1992). We assign counterparts to 3798 sources with Likelihood Ratio . At that cut the spurious identification rate is about 6% and the total number of 0.5-8 keV detected sources with optical counterparts is 3798 (see Table 1).

Redshifts for the XMM-XXL X-ray sources are from several follow-up spectroscopic campaigns. The XMM-XXL field overlaps with the SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al., 2013) programme, which provides spectroscopy for UV/optically selected broad-line QSOs and luminous red galaxies. Stalin et al. (2010) presented spectroscopy for X-ray sources selected in the original XMM-LSS survey (Clerc et al., 2014), which is part of the equatorial XMM-XXL survey field. Most of the redshifts however, are from a total of five special SDSS plates dedicated to follow-up spectroscopy of X-ray sources as part of the Ancillary Programs of SDSS-III. The overlap between those plates and the XMM-XXL survey region is . Targets were selected to have and , where is either the SDSS PSF magnitude in the case of optically unresolved sources (SDSS type=6) or the SDSS model magnitude for resolved sources. Specific details on these spectroscopic observations, including spectral classification, visual inspection and redshift quality assessment are presented by Menzel et al. (2015). The total number of XMM-XXL sources with secure spectroscopic measurements are presented in Table 1. Also shown in this table are the number of sources with spectroscopic redshifts in the interval .

Although when constructing the X-ray source catalogue of the XMM-XXL field fluxes are estimated for a power-law X-ray spectrum with spectral index , in the rest of the analysis we adopt for the calculation of fluxes, luminosities and sensitivity maps. This is because at the depth of the XMM-XXL survey AGN at are powerful QSOs with . The fraction of obscured AGN among such luminous sources is a decreasing function of luminosity (Ueda et al., 2003; Akylas et al., 2006; Merloni et al., 2014; Ueda et al., 2014; Buchner et al., 2015) and therefore a spectral index of , which represents the intrinsic unobscured power-law X-ray spectrum of both local Seyferts (e.g Nandra & Pounds, 1994) and luminous high-redshift QSOs (e.g. Vignali et al., 2005; Shemmer et al., 2005; Just et al., 2007) is appropriate for this population. The 0.5-10 keV fluxes estimated for are about 35% fainter than those for . The XMM-XXL sensitivity curve in the 0.5-10 keV band is shown in Figure 1 and is estimated following methods described in Georgakakis & Nandra (2011). In the calculation of the sensitivity curve we only consider the overlapping area between the SDSS-III Ancillary Programs spectroscopic plates used to target X-ray sources and the XMM-XXL survey region. We also take into account the flux limit for follow-up spectroscopy . This limit appears as a smooth drop in area rather than a sharp cut in Figure 1 because the Poisson probability of measuring a flux above this limit is used to determine the sensitivity curve.

3 Methodology

3.1 X-ray Luminosity function estimation: Chandra survey fields

A Bayesian approach is adopted for the determination of the X-ray luminosity function. The X-ray sources detected in a survey are essentially Poisson realisations of a parent sample and therefore the likelihood can be written as the product of the Poisson probabilities of individual sources. Following the works of Marshall et al. (1983), Loredo (2004), Aird et al. (2010) and Buchner et al. (2015) the likelihood can be written as

 L(di|θ)=e−λ×N∏i=1∫dlogLXdVdzdzp(di|LX,z)ϕ(LX,z|θ), (1)

where is the comoving volume per solid angle at redshift , signifies the dataset and represents the parameters of the luminosity function model, , that are to be estimated. The multiplication is over all sources, , and the integration is over redshift and X-ray luminosity. The quantity is the probability of a particular source having redshift and X-ray luminosity . This captures uncertainties in the determination of both redshifts (e.g. photometric redshifts measurements) and X-ray fluxes because of Poisson statistics and the Eddington bias. In equation 1, is the expected number of detected sources in a survey for a particular set of model parameters

 λ=∫dlogLXdVdzdzA(LX,z)ϕ(LX,z|θ). (2)

where, is the sensitivity curve that quantifies the survey area over which a source with X-ray luminosity and redshift (and hence flux ) can be detected. Note that the selection function term, , is not included within the integral of equation 1. The reader is referred to an extensive discussion in Loredo (2004) on that point.

The goal of this paper is the estimation of the X-ray luminosity function in the redshift interval . However, defining a sample of X-ray AGN in a relatively narrow redshift range is not straightforward. This is because the redshifts of many sources are determined by photometric methods and therefore have uncertainties, which are not negligible compared to the size of the redshift interval. It may happen for example, that the errors of the photometric redshift of a particular source straddle one (or both) boundaries of the redshift range of interest (see Fig. 3).

We deal with this difficulty by simply using all sources in the X-ray sample and splitting the luminosity function model into two terms

 Unknown environment '%' (3)

The first term refers to the luminosity function within the redshift interval of interest, , and has its own set of parameters, . The second term, , corresponds to the X-ray AGN space density outside that redshift range and has a different set of parameters, , which are treated as nuisance parameters. The advantage of this approach is that it allows an estimate of the AGN space density at nearly independent of the shape and form of the evolution of the X-ray luminosity function at lower redshift. In any X-ray selected sample the bulk of the AGN population lies at low redshift . This may introduce systematics in the determination of the AGN space density at , if the same evolutionary law is fit to the data across all redshifts. In this case it is possible that the determination of the model parameters is dominated by the regions of the redshift/luminosity parameter space with the most data. Our approach minimises the impact of this potential source of bias. In our analysis the term is modelled as a step function, i.e. the sum of constants that correspond to the AGN space densities at different luminosity and redshift bins. The values of these constants are determined by the data via equation 1. For the calculation of equation 1 assumptions need to be made on the redshift errors. For sources with photometric redshift determinations we adopt the corresponding redshift Probability Distribution Function (PDZ) as a measure of the uncertainty. Spectroscopic redshifts in the sample have reliabilities and therefore their corresponding PDZs are assumed to be delta functions at the spectroscopic redshift of the source. Sources without optical counterparts are assigned a flat PDZ in the redshift range (see Section 2.1).

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 is adopted in this calculation. The flux distribution in the 0.5-10 keV band is then convolved with the 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 . The two-dimensional probability distribution in and is the term of equation 1. In practice we use importance sampling (Press et al., 1992) to evaluate the integral of equation 1. For each source we draw and samples based on the PDZ and Poisson X-ray counts distribution of that source. The luminosity function is then evaluated for each sample point, . The integral of equation 1 is simply the average luminosity function of the sample.

3.2 X-ray Luminosity function estimation: XMM-XXL field

In the case of the XMM-XXL field, only spectroscopically confirmed sources in the redshift interval are used. For that sample there is no need to apply equation 3 to determine the corresponding AGN space density. The limitation of that sample however, is that it is both X-ray and optical flux limited because of the magnitude limit of the spectroscopic follow-up observations. Both cuts need to be accounted for to infer the X-ray luminosity function via equation 1. We do that by exploiting the fact that the XMM-XXL sample consists of powerful [] broad-line QSOs. We use the well-established correlation between monochromatic X-ray [] and UV [] luminosities of broad-line QSOs (e.g. Steffen et al., 2006; Just et al., 2007; Lusso et al., 2010) to link the observed optical magnitudes and X-ray fluxes of the sample and account for the selection effects. Figure 4 demonstrates the correlation between vs using X-ray selected broad-line QSOs from the XMM-XXL field in the redshift interval . The low redshift cut is to avoid X-ray AGN with relatively low luminosities, for which the UV/optical continuum shows non-negligible contribution from the host galaxy. We adopt the Lusso et al. (2010) bisector best-fit . At a given monochromatic optical luminosity we assume that the data points scatter around the above relation following a Gaussian distribution with standard deviation . Figure 5 shows that this is a reasonable assumption. From that figure we estimate . Broad Absorption Line (BAL) QSOs represent up to about 26% of optically selected samples (Hewett & Foltz, 2003; Reichard et al., 2003; Gibson et al., 2009), and are known to be X-ray faint either because of absorption or intrinsic X-ray weakness (e.g. Gallagher et al., 2006; Luo et al., 2014). Such sources do not appear to skew the distributions plotted in Figures 4 and 5. This is likely related to the bright X-ray flux limit of the XMM-XXL survey, which selects against X-ray faint Under assumptions on the optical Spectral Energy Distribution of QSOs, Figure 5 can be used to estimate the SDSS -band optical magnitude distribution of X-ray sources of given and and then determine the fraction of this distribution that is brighter than the spectroscopic magnitude limit of the survey, . In this case the expected number of detected sources with the surveyed area, , in equation 1 can be rewritten

 λ=∫dlogLXdVdzdzA(LX,z)B(LX,z,| r)×ϕ(LX,z|θ)η(r), (4)

where is the SDSS -band magnitude distribution of a source with and . The efficiency factor is the success rate of measuring secure redshifts for X-ray sources as a function of the SDSS -band magnitude. It accounts for the fact that not all X-ray sources with secure counterparts have successful spectroscopic redshift measurements. Collisions between SDSS fibers or the finite number of science fibers on the SDSS spectroscopic plates mean that not all candidate sources for follow-up spectroscopy can be assigned a fiber. Moreover, for the sources that are observed the rate of secure redshift measurement depends on their optical brightness. At fainter magnitudes the signal-to-noise ratio of the optical spectra decreases and therefore the ability to estimate redshifts is affected. The probability of a source being assigned a fiber is random, while the redshift success rate depends, at least to the first approximation, on optical magnitude. The factor is the number of spectroscopically confirmed X-ray sources in the magnitude interval divided by the total number of X-ray sources that are potential targets for follow-up spectroscopy in the same magnitude range. The parent X-ray sample of potential targets is selected to have and  mag (see Section 2.2). The -band magnitude dependence of is shown in Figure 6. For magnitudes in the range  mag the efficiency factor is nearly constant and larger than 80%. This fraction drops to about 50% at , the limiting magnitude for follow-up spectroscopy with the Sloan telescope, and is zero for  mag.

In equation 4 for the calculation of the optical k-corrections we use the simulated QSO SEDs of McGreer et al. (2013). They are generated assuming a double power-law continuum in the rest-frame UV/optical part of the SED with a break-point at . The short- and long-wavelength slopes are drawn from normal distributions with means of and , respectively. For both slopes the scatter is fixed to 0.3. Emission lines with luminosity dependent equivalent widths as well as Lya forest absorption are also added to the simulated SEDs as described in McGreer et al. (2013). A total of 60 000 model SEDs are generated in the redshift interval and for AGN luminosities . These are used to calculate the expected distribution of the observed -band optical magnitudes for a given redshift and AGN luminosity, and determine the term in equation equation 4. For the calculation of X-ray k-corrections we assume a power-law X-ray spectrum with .

3.3 X-ray luminosity function models

The Bayesian framework outlined above explicitly requires a model with a set of free parameters that are constrained by the observations. We consider both non-parametric and parametric models for the X-ray luminosity function in the redshift range . The non-parametric model simply assumes that the space density of AGN is constant within a given luminosity and redshift interval. In this case the free model parameters to be estimated are the AGN space densities in each luminosity and redshift bin. This is equivalent to the widely used 1approach for the determination of binned luminosity functions. The advantage of the non-parametric approach is that it allows investigation of the form and amplitude of the X-ray luminosity function evolution in a model-independent way.

We also consider four parametric models for the X-ray luminosity function and its redshift evolution, which have been extensively used in the literature. These are the Pure Luminosity Evolution (PLE), Pure Density Evolution (PDE), Luminosity Dependent Density Evolution (LDDE) and Luminosity And Density Evolution (LADE) models. The parametrisation of each model follows Vito et al. (2014). The X-ray luminosity function in the redshift range is defined as the space density of AGN per logarithmic luminosity bin and is described by a double power-law of the form

 Missing or unrecognized delimiter for \Big (5)

where is the normalization, and are the faint and bright-end power-law slopes, respectively, and is the break luminosity. The PLE model assumes that is a function of redshift and evolves according to the relation

 L⋆(z)=L⋆(z0)×(1+z1+z0)p, (6)

where we fix and the exponent parametrises how fast the break luminosity evolves with redshift. In the case of PDE it is assumed that only the normalisation of the luminosity function evolves with the redshift according to the relation

 K(z)=K(z0)×(1+z1+z0)q, (7)

where and the exponent parametrises the speed of the normalisation factor evolution. The LADE model adopted here is similar to the Independent Luminosity and Density Evolution (ILDE) model described by Yencho et al. (2009). We use equation 6 to parametrise the redshift evolution of and also add a normalisation evolution term of the form

 K(z)=K(z0)×10q(z−z0). (8)

Finally we consider the LDDE parametrisation of Hasinger et al. (2005), where the normalisation factor of equation 5 changes with redshift as

 K(z)=K(z0)×(1+z1+z0)q+β(logLX−44), (9)

where as previously and the rate of the density evolution also depends on the X-ray luminosity via the parameter .

Recall that in addition to the model parameters that describe the X-ray luminosity function in the redshift range the likelihood function in equation 1 also includes terms that correspond to the total AGN space density outside the redshift range of interest, or (see equation 3). A non-parametric model is adopted to describe that term. This to minimise the impact of parametric model assumptions to the results and conclusions. We use 3 redshift bins in the range with size . Each redshift bin is split into 5 logarithmic luminosity bins in the range (units of ) with width  dex. One additional term is used to model the luminosity function in the range . A constant AGN space density is assumed within the above luminosity and redshift intervals. We therefore use a total of 16 nuisance parameters to model the AGN X-ray luminosity function at or .

The MultiNest multimodal nested sampling algorithm (Feroz & Hobson, 2008; Feroz et al., 2009) is used for both Bayesian parameter estimation and the calculation of the Bayesian evidence, , of each model, i.e. the integral of the model likelihood over the parameter space allowed by the priors. The Bayesian evidence is used for model comparison, i.e. to select among different models the one that better describes the data.

3.4 Potential obscuration effects

The analysis presented in this paper uses the full-band selected sample to determine the X-ray luminosity function of AGN at . This is because of the higher sensitivity of the 0.5-10 keV band compared to e.g. the 2-10 keV band, which translates to a larger number of sources. A potential issue however, is that the analysis described above ignores the impact of obscuration on the observed AGN flux. Figure 7 shows that for AGN at moderate line-of-sight columns, (), suppress the observed 0.5-10 keV flux of AGN by less than about 25%. Higher column densities however, have a larger impact on the 0.5-10 keV flux and therefore result in incompleteness in the sample that is not accounted for in our analysis. Nevertheless, we are primarily interested in the differential evolution of AGN in the redshift intervals and . Figure 7 shows that for fixed obscuration and intrinsic AGN luminosity there is little difference in the level of flux suppression between redshifts and . Incompleteness related to obscuration is therefore expected to be similar across the redshift range . This allows direct comparison of the inferred AGN space densities in the redshift bins and , under the assumption that the distribution of AGN in obscuration does not change dramatically between these redshift intervals. Recent studies on the obscuration distribution of AGN (Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015) show that the obscured AGN fraction increases with redshift, at least to . In Buchner et al. (2015) for example, the obscured fraction in Compton thin AGN is about 50% at and increases to 70% at . At higher redshifts however, current constraints on the evolution of the obscured AGN fraction are still limited by small number statistics. Nevertheless, there are suggestions that the obscured AGN fraction remains roughly constant with redshift at (Buchner et al., 2015).

Obscuration related effects potentially have a larger impact on the analysis of the XMM-XXL data. For these sources it is explicitly assumed that their optical/UV continua are not obscured or reddened. This assumption allows us to place them on the correlation for type-I AGN and correct for the spectroscopic magnitude cut as explained above. This poses a problem when comparing the Chandra with the XMM-XXL constraints on the X-ray luminosity function. Nevertheless, the shallow XMM-XXL survey is only sensitive to powerful QSOs [] at . At lower redshift at least, it is well established that the obscured AGN fraction is rapidly decreasing with increasing accretion luminosity (e.g. Ueda et al., 2003; Akylas et al., 2006; Ueda et al., 2014; Buchner et al., 2014; Merloni et al., 2014). Obscuration incompleteness corrections are therefore likely to be small for luminosities . We explore this assumption further in the next section by directly comparing the Chandra and XMM-XXL constraints in overlapping luminosity bins.

4 Results

4.1 The non-parametric X-ray luminosity function determination

Figure 8 plots the non-parametric estimate of the X-ray luminosity function in two redshift intervals, and . For the low redshift sub-sample, , we also plot separately the AGN space density determined independently from the Chandra deep fields and the XMM-XXL survey. As expected these datasets probe different luminosity ranges. The Chandra fields provide constraints to (, but do not probe sufficient volume to detect luminous and rare sources with luminosities . The XMM-XXL survey can constrain the bright-end of the luminosity function but is not sensitive enough to detect AGN below . Nevertheless, in the interval () both the Chandra and XMM data provide meaningful constraints to the AGN space density. These independent non-parametric estimates are in good agreement within the 68% errors plotted in Figure 8. This suggests that obscuration effects do not have a strong impact on the AGN space densities estimated from the XMM-XXL data. In Figure 8 the panel does not plot separately the space density estimates from the Chandra and then XMM-XXL survey fields. This is because at this redshift interval there is virtually no overlap in the luminosities probed by the two datasets. The XMM-XXL constrains only the brightest luminosity bin, (units ). The Chandra fields include only sources that are fainter than .

A striking result illustrated in Figure 8 is the strong evolution of the AGN population between redshifts and . This is further demonstrated in Figure 9, which plots as a function of luminosity the ratio of the non-parametric AGN space densities in the redshift intervals and . At luminosities there is a factor of 5 decrease in the AGN space density from to . There is also evidence for luminosity dependent evolution. AGN with in Figure 9 appear to evolve faster than intrinsically brighter sources. We quantify the statistical significance of the evidence for luminosity dependent evolution using the quantity

 R=logϕ(LX=1045−1046,z=4−5)ϕ(LX=1045−1046,z=3−4)−logϕ(LX=1043−1045,z=4−5)ϕ(LX=1043−1045,z=3−4), (10)

where is the binned (non-parametric) luminosity function. is the logarithmic difference between the high luminosity () data point of Figure 9 and the sum of the two moderate luminosity data points () of the same plot. We bin together moderate luminosity AGN to simplify the problem and also improve the statistical reliability of the results. We use the full probability density distribution of and find that at the 90% probability , i.e. there is differential evolution between moderate () and powerful () X-ray AGN from to . We further investigate this using a simple parametric approach. A double power-law function (equation 5) is fit independently to the data in the redshift bins and . In this exercise evolutionary effects within each of the two redshift intervals are ignored. The best-fit parameters are presented in Table 2. The ratio between the two double power-laws in the redshift bins and is plotted with the shaded region in Figure  9. The apparent increase of this ratio toward faint luminosities is because the faint-end slope of the sample is poorly constrained and on the average steeper than that of the sample. Nevertheless, in the interval , where constraints on the AGN space density are available for both the and the sub-samples, the ratio between the two power-laws is consistent with luminosity-dependent evolution. Larger samples are needed however, particularly in the redshift interval , to reduce the uncertainties in Figure 9 and further explore the evidence for luminosity dependent evolution.

The amplitude of the AGN X-ray luminosity function evolution between the redshift intervals and in Figure 9 is further compared with the predictions of luminosity function parametrisations that include an exponential decline at (e.g. Gilli et al., 2007). Such models are motivated by the rapid evolution of optical QSO space density at high redshift (e.g. Schmidt et al., 1995; Richards et al., 2006). We use the Gilli et al. (2007)2 LDDE parametrisation of the X-ray luminosity function with and without an additional exponential cutoff to predict the space density of AGN at and . The ratio between these predictions is plotted in Figure 9. The amplitude of the AGN evolution inferred in this paper for luminosities is consistent with the Gilli et al. (2007) LDDE parametrisation that includes an exponential cutoff. AGN with luminosities in the range lie in between the Gilli et al. (2007) model predictions with and without an exponential cutoff. This is suggestive of milder evolution, albeit at the confidence level.

4.2 Parametric X-ray luminosity function determination

Next we use Bayesian model comparison to assess which of the evolutionary models outlined in Section 3.3 provides a better description of the observations in the redshift interval . The PLE, PDE, LADE and LDDE parametric models of Section 3.3 are fit to the combined Chandra and XMM-XXL dataset. Table 3 presents the best-fit parameters for each parametric model. The Bayes factor (ratio between evidences) of each model relative to the model with the highest evidence (PDE) is also shown in that table.

The model with the highest evidence in Table 3 is the PDE, with the second best being the LDDE. The Bayes factor of the two models is . Based on the Jeffreys interpretation of the Bayes factor (Jeffreys, 1961) this difference suggests that both models describe equally well the evolution in the redshift interval of the X-ray selected AGN sample presented in this paper. The bulk of the AGN in the present sample have X-ray luminosities . Figure 9 shows that such sources experience similar reduction in their space density between and , i.e. consistent with pure density evolution. More luminous AGN [], which appear to experience milder evolution in Figure 9, represent only a small fraction of the present sample. This combined with the fact that the PDE is a simpler model compared to the LDDE (5 vs 6 free parameters) results in similar Bayesian evidences for these two parametric models.

Table 3 also shows that the LADE parametrisation adopted in this work performs worse than the LDDE. The Bayes factor of the two models is . This difference is strong evidence in favor of the LDDE Jeffreys (1961).

The one model that performs significantly worse than the rest is the PLE. The Bayes factor of that model relative to the one with the highest evidence (PDE) is . Based on the Jeffreys interpretation of the Bayes factor (Jeffreys, 1961) this is decisive evidence against the PLE model.

Although the PDE model is favoured by our analysis for the evolution of AGN in the redshift interval in the next sections we use the LDDE model to compare with previous studies and to determine the contribution of X-ray AGN to the ionisation of the Universe. This is because of the small difference in the evidences of the PDE vs the LDDE and the long literature on the LDDE evolutionary model. Nevertheless, the results and conclusions are not sensitive to the particular choice of AGN evolutionary model.

4.3 Comparison with previous studies

A number of studies on the X-ray luminosity function of AGN have appeared in the literature recently. These include works focusing on the space density of X-ray AGN at high redshifts (Civano et al., 2010; Kalfountzou et al., 2014; Vito et al., 2014) and results on the global X-ray luminosity function evolution across redshift (Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015; Miyaji et al., 2015). In this section we compare our X-ray luminosity function with previous studies using X-ray samples selected in the 0.5-2 or 0.5-10 keV bands, i.e. similar to the one presented in this paper. Vito et al. (2014) compiled one of the largest samples to date in the redshift interval . They select sources in the 0.5-2 keV band and also present non-parametric () estimates of the X-ray luminosity function. These aspects of the Vito et al. analysis methodology are similar to ours. We also use recent results of Aird et al. (2015) as a representation of parametric approaches to fit the X-ray luminosity function of AGN to the full redshift interval accessible to current X-ray selected samples, .

Vito et al. (2014) use a sample similar to ours in size to determine the X-ray luminosity function in the range . They use both photometric and spectroscopic redshifts (total of 141) and apply corrections to account for sources without photometric redshift determinations (total of 65). We compare our best-fit LDDE model with their results in Figure 10. The redshift intervals of each panel are the same as in Figure 7 of Vito et al. (2014). We plot their non-parametric estimates that include their redshift incompleteness corrections. Our LDDE parametric models in Figure 10 is estimated at the middle of each redshift interval. We find that our X-ray luminosity function determination is systematically lower than the Vito et al. (2014) data points. This discrepancy ( significance in e.g. the redshift interval of Figure 10) may be related to the fact that Vito et al. (2014) account for X-ray obscuration in the determination of X-ray luminosities and the calculation of the of individual sources. A total of 36 of their 141 sources have column densities , which translates to suppression of their fluxes compared to the unobscured () case by % (see Figure 7). Such corrections are ignored in the analysis presented here. Alternatively the discrepancy may be related to how sources with photometric redshifts or sources without any redshift information are treated. Vito et al. (2014) use only the best-fit photometric redshift solution without taking into account the corresponding uncertainties. It is further assumed that the 65 sources in their sample without photometric redshifts all lie in the redshift interval and that they follow the redshift distribution of the X-ray sources with redshift determinations (photometric or spectroscopic) in the range . The amplitude of this correction is and  dex increase at luminosities and respectively. In that respect Vito et al. (2014) determine maximal X-ray luminosity functions. It is likely that some of the X-ray sources without redshift information lie outside the redshift range . Heavily obscured AGN or AGN hosted by early type hosts at moderate redshifts, , may also have red SEDs, similar to those of high redshift sources (Koekemoer et al., 2004; Schaerer et al., 2007; Rodighiero et al., 2007; Del Moro et al., 2009).

We further investigate this issue by adopting a methodology similar to that of Vito et al. (2014) to determine the X-ray luminosity function. We use the best-fit photometric redshifts only, i.e. ignoring the photometric redshift probability distribution functions. Sources without optical identifications in the sample are assigned random redshifts in the range based on the redshift distribution of sources with spectroscopic or photometric redshift measurements. The results are presented in Figure 11. It shows that an approach that ignores photometric redshift errors results in an overestimation of the AGN space density. This is not only because of sources without optical identifications and therefore without photometric redshift (total of 107 in our sample) estimates being forced to lie in the range . The overestimation is mainly the result of ignoring photometric redshift uncertainties. In the likelihood equation 1 the probability of a source having redshift and luminosity is weighted by the luminosity function. Sources with broad (see Fig. 3) redshift probability distribution functions are more likely to lie at low redshift and moderate luminosities simply because the space density of AGN is higher there. Ignoring this weighting in equation 1 by e.g. fixing photometric redshifts to a single (best-fit) value overestimates the luminosity function at high redshifts.

Also plotted in Figure 10 are the flexible double power-law parametric model of Aird et al. (2015) for their 0.5-2 selected sample before applying corrections for obscuration (i.e. directly comparable to our analysis). The binned X-ray luminosity function estimates of Aird et al. (2015) are also shown. These data points are estimated using the method developed by Miyaji et al. (2001) and are therefore not independent of the underlying parametric model plotted in Figure 10. Nevertheless, these binned estimates provide a measure of the associated uncertainties and can be used to identify redshift and luminosity intervals where the parametric model provides a poorer fit to the data. The Aird et al. (2015) X-ray luminosity function estimated from their 0.5-2 keV selected sample is in fair agreement with our results. This may be not surprising given the overlap between the two datasets and the similar approaches. Figure 12 plots the X-ray luminosity density evolution of AGN using our LDDE model and the Aird et al. (2015) total X-ray luminosity function that includes both obscured (Compton thin and thick) and unobscured sources. The X-ray luminosity function estimated in this work accounts for about 40% of the total X-ray luminosity density determined by Aird et al. (2015).

Finally, we also compare in Figure 10 the X-ray luminosity function with high redshift determinations of the optical QSO luminosity function. The conversion from UV/optical to X-rays ultimately depends on the scatter in the relations between bolometric and X-ray or UV luminosities. There are suggestions that the bolometric–to–X-ray luminosity ratio has a larger scatter than the bolometric–to–UV luminosity ratio (Hopkins et al., 2007). We therefore convert from UV/optical to X-rays by convolving the best-fit parametric model of UV/optical QSOs from the relevant publications with the relation of Lusso et al. (2010) assuming a scatter of 0.4 dex (e.g. Figure 5). This calculation also requires assumptions on the shape of the UV SED of QSOs at wavelengths  Å. We assume a power-law of the form (Vanden Berk et al., 2001; Telfer et al., 2002). A steeper slope, (Lusso et al., 2015), translates to an X-ray luminosity function that is offset by  dex compared to . At the and panels of Figure 10 we overplot the Masters et al. (2012) best-fit double power-law models at and . At the panel of Figure 10 we show the Ross et al. (2013) LEDE (Luminosity Evolution and Density Evolution) model fit to the BOSS Stripe82 QSO data (Palanque-Delabrouille et al., 2011). We choose not to extrapolate this model past redshift , the limiting redshift of the BOSS Stripe82 QSO sample. At the redshift interval (mean redshift ) we transform to X-rays the optical QSO luminosity function determined by McGreer et al. (2013). The optical luminosity functions are plotted by shaded regions in Figure 10. The shape and normalisation of the UV/optical QSO luminosity functions in Figure 10 at luminosities are sensitive to the scatter of the relation used in the convolution.

A striking result from Figure 10 is the steep faint-end slope of the Masters et al. (2012) best-fit double power-law models at and compared to the faint-end slope of the X-ray luminosity function. The ratio between the UV/optical and the X-ray luminosity functions is a proxy of the type-I fraction among AGN, , and can be used to explore the luminosity dependence of this quantity at high redshift. In this section we define type-I AGN as sources with blue UV/optical continua, typical of broad-line QSOs. We use the Masters et al. (2012) luminosity function and the LDDE parametrisation for the X-ray luminosity function to determine and plot the results in the case of in Figure 13. The luminosity dependence of the at is very similar and is not shown. We find evidence that the type-I AGN fraction at is a non-monotonic function of luminosity. There is a minimum at , followed by a steep increase toward fainter luminosities.

The behaviour of in Figure 13 is opposite to studies that find a drop in the type-I or X-ray unobscured () AGN fraction with decreasing luminosity or equivalently that type-II or obscured () AGN dominate at faint luminosities (Ueda et al., 2003; Akylas et al., 2006; Ueda et al., 2014; Merloni et al., 2014; Aird et al., 2015). Studies on the obscuration distribution of AGN in the local Universe support a picture where the obscured AGN fraction increases with decreasing luminosity. Nevertheless, they also find evidence for a turnover (drop) of the obscured AGN fraction at very faint X-ray luminosities, below about (e.g Burlon et al., 2011; Brightman & Nandra, 2011b). Recently, Buchner et al. (2015) extended these results to higher redshift and found evidence that the obscured AGN fraction peaks at a redshift-dependent luminosity and then drops at both brighter and fainter luminosities. Figure 13 overplots the obscured AGN fraction derived by Buchner et al. (2015) in the redshift interval . For this comparison we assume that . This is a simplistic assumption because the definition of Type-I QSOs in the case of UV/optical samples and unobscured/obscured AGN in X-ray samples like in Buchner et al. (2015) is different. Nevertheless, to the first order the quantity should be at least loosely related to . In Figure 13 there is qualitative agreement between the Buchner et al. (2015) parameter and our definition of .

At bright luminosities the quantity in Figure 13 increases with increasing . This in agreement with previous studies based on X-ray (Ueda et al., 2014, e.g.), optical (e.g. Simpson, 2005) or infrared (e.g. Assef et al., 2013) data. At the brightest luminosities probed by our data, , the type-I fraction is about . This number is relevant to the population of powerful QSOs (bolometric ) with reddened UV/optical continua [extinction ] identified in recent wide-area infrared surveys (Stern et al., 2014). These reddened/obscured sources correspond to X-ray luminosities (Marconi et al., 2004) and are suggested to represent up to 50% of the QSO population at these luminosities (Assef et al., 2014). We caution however, that there are uncertainties on the inferred bolometric luminosities of these sources and hence, the obscured AGN fraction, depending on the assumed geometry, physical scale and covering fraction of the obscuring material (Assef et al., 2014). Additionally, the X-ray properties of these infrared selected AGN are poorly known. There are suggestions that they represent heavily obscured, even possibly Compton thick, systems (Stern et al., 2014), in which case they are expected to be underrepresented in our sample.

4.4 Contribution of AGN to ionisation of the Universe

The X-ray luminosity function can be used to set limits on the contribution of the AGN to the radiation field that keeps the Universe ionised at redshift . The advantage of X-ray selection is that it provides a better handle on the faint-end of the AGN luminosity function compared to UV/optically selected samples. The downside is that assumptions on the escape fraction of UV photons have to be made to convert AGN space densities to ionising photon densities. At the very least however, X-rays surveys can set strict upper limits on the contribution of AGN to the ionising photon field, under the assumption that all photons emitted by the central source escape to the Inter-Galactic Medium.

The rate of hydrogen ionising photons is estimated by integrating the AGN Spectral Energy Distribution in the energy range 1-4 ryd

 ˙n=∫4ryd1rydL(ν)hνdν, (11)

where is the AGN monochromatic luminosity at frequency . The Lusso et al. (2010) relation is used to convert the 2-10 keV X-ray luminosity to UV monochromatic luminosity at 2500 Å. We then extrapolate to the wavelength range 227Å  (4 ryd) and 910Å  (1 ryd) assuming a double power law for the AGN SED of the form

 L(ν)∝{ν−0.5(1100\,\AA<λ< 2500\,\AA)ν−1.7λ< 1100\,\AA (12)

The spectral slopes in the above relation are from Vanden Berk et al. (2001) and Telfer et al. (2002). At any given redshift the comoving density of the hydrogen ionising rate is then estimated by integrating the X-ray luminosity function

 N=∫ϕ(LX,z)fesc(LX,z)˙ndlogLX. (13)

The LDDE parametrisation of the X-ray AGN luminosity function is adopted. The integration limits are set to and . The photon escaping factor, , accounts for the fraction of obscured AGN, which likely depends on both redshift and accretion luminosity. In these sources the ionising photons are absorbed locally and therefore do not contribute to the cosmic ionisation radiation field. Our analysis does not constrain the distribution of AGN in obscuration. Additionally at redshift there are still discrepancies among different studies on the obscured AGN fraction and its dependence on luminosity (e.g. Ueda et al., 2014; Buchner et al., 2015; Aird et al., 2015). For example, in Figure 13 we find evidence that the type-I AGN fraction and by proxy the obscured AGN fraction, have a different dependence on luminosity compared to relations established at lower redshift. We choose to present our baseline results for the luminosity dependent type-I AGN fraction determined by Merloni et al. (2014). In that study type-I refers to AGN with either blue UV/optical continua and/or broad emission lines. This definition is more releavant to the determination of the hydrogen ionising photon rates of AGN, compared to e.g. the standrard X-ray unobscured AGN definition, . Adopting the Merloni et al. (2014) relation for also allows direct comparison with previous X-ray studies that also used monotonically increasing obscured AGN fraction with decreasing X-ray luminosity to approximate the escape fraction of UV photons in high redshift AGN. Figure 14 plots as a function of redshift for our baseline model for the photon escaping fraction. In that figure we also place an upper limit in by setting , i.e. the extreme case that all photons escape and contribute to the hydrogen ionising radiation. This translates to a net increase of by a factor of about 1.4 compared to the Merloni et al. (2014) type-I AGN relation for . We further explore how the results in Figure 14 change if is approximated by X-ray definitions of the unobscured () AGN fraction. We adopt the luminosity dependence of the unobscured X-ray AGN fraction derived by Buchner et al. (2015) in the redshift interval (see Fig. 13). This assumption yields photon rate densities that are a factor of about 1.7 smaller compared to the baseline results that use the Merloni et al. (2014) type-I AGN fraction as proxy of . For clarity these results are not plotted in Figure 14.

The constraints above should be compared to the minimum photon rate density required to keep the Universe ionised at a given redshift. This is estimated from equation (26) of Madau et al. (1999) after scaling it to our cosmology. We also assume that the ionised hydrogen clumping factor in that relation, which is a measure of the inhomogeneity of the medium, evolves with redshift as

 C=1+43×z−1.71. (14)

This relation is based on cosmological simulations (Pawlik et al., 2009) and is adopted by Haardt & Madau (2012) to synthesise the evolving spectrum of the diffuse UV/X-ray radiation field. Figure 14 plots the redshift dependence of the photon rate density required to keep the Universe ionised. This figure shows that for the baseline model (i.e. Merloni et al., 2014, as proxy of ) the AGN contribution to the photon rate density required to keep the Univese ionised decreases from 70% at to about 20% at . Assuming the fractions above translate to upper limits at and at . These numbers are in broad agreement with some previous studies on the role of X-ray AGN in the ionisation of the Universe (Barger et al., 2003a; Haardt & Madau, 2012; Grissom et al., 2014). The decreasing ionising photon rate density fractions with increasing redshift is a direct consequence of the strong evolution of the X-ray AGN luminosity function between redshifts and .

Finally in Figure 14 the model constraints on the ionising photon rate density at from our X-ray sample are compared to previous works based on either UV/optical QSO surveys (Glikman et al., 2011; Masters et al., 2012; McGreer et al., 2013) or UV/X-ray selected samples (Giallongo et al., 2015). For the calculation of the ionising photon rate densities the luminosity functions estimated in these works are integrated between absolute magnitudes and  mag by adopting the UV spectrum of Equation 12. Our baseline model assuming the Merloni et al. (2014) escaping fraction is consistent with the lower-normalisation UV/optical data points in Figure 14.

5 Discussion

In this paper we explore the evolution of the AGN X-ray luminosity function in the redshift interval by combining deep Chandra surveys with the wide-area and shallow XMM-XXL sample. This dataset provides sufficient depth and volume to determine the space density of AGN over 3 dex in luminosity, , in the redshift intervals and . The analysis methodology we develop takes into account Poisson errors in the determination of X-ray fluxes and luminosities as well as uncertainties in photometric redshift measurements. We demonstrate that the latter is critical for unbiased measurements of the AGN space density at high redshift. Ignoring photometric redshift errors overestimates the X-ray luminosity function. We also choose to follow a non-parametric approach and determine the space density of AGN in luminosity and redshift bins. This allows us to explore the evolution of the AGN population independent of model assumptions. Additionally when a model parametrisation is applied to AGN at all redshifts the fit may be driven by the redshift and luminosity intervals that contain most data. This can introduce biases at high redshift, , where AGN samples are typically small as a result of the rapid evolution of the AGN population as well as survey sensitivity and volume limitations. We account for this potential issue by splitting the luminosity function into independent redshift components, i.e. , .

We confirm previous studies for a strong evolution of the AGN population in the redshift interval (Brusa et al., 2009; Vito et al., 2013; Civano et al., 2012; Kalfountzou et al., 2014). We also find suggestive evidence for luminosity dependent evolution of the X-ray luminosity function. The space density of AGN in the luminosity interval decreases faster than more luminous sources between redshifts and , albeit at the 90% significance level. A similar evolution pattern is also observed in the optical luminosity function of QSOs between redshifts and (Masters et al., 2012). Our finding can be interpreted as evidence that the formation epoch of the most powerful QSOs [] precedes that of lower luminosity systems. This is similar to AGN downsizing trends established at lower redshifts (Ueda et al., 2003; Hasinger, 2008; Aird et al., 2010; Miyaji et al., 2015). A strong evolution of the faint-end of the AGN luminosity function with increasing redshift is consistent with the absence of X-ray selected AGN at in the CANDELS (Grogin et al., 2011; Koekemoer et al., 2011) subregion of the Chandra Deep Field South (Weigel et al., 2015). Extrapolating the LDDE model of Table 3 we predict AGN at redshift in that field.

Our analysis also places limits on the contribution of AGN to the UV photon field needed to keep the hydrogen ionised at high redshift. Using empirical relations for the type-1 AGN fraction as a function of luminosity (Merloni et al., 2014) we show that AGN dominate or at least contribute a sizable fraction of the required UV photons to redshift . At higher redshift the evolution of the X-ray luminosity function translates to a decreasing contribution of X-ray AGN to the UV photon field required to keep the hydrogen ionised. The extreme assumption of a photon escaping fraction of unity for all AGN sets an upper limit of 30% to the contribution of AGN to the UV photon rate density required to keep the hydrogen ionised at . Barger et al. (2003a) use multicolour optical data in the 2 Ms Chandra Deep Field North (Barger et al., 2003b) and conclude that the X-ray selected AGN candidates at are too few to ionise the intergalactic medium at those redshifts. Haardt & Madau (2012) estimated the contribution of AGN to hydrogen ionisation rate using the Ueda et al. (2003) X-ray luminosity function and AGN obscuration distribution. They found that AGN do not play an important role as a source of ionising photons above redshifts . Grissom et al. (2014) determine the contribution of AGN to the ionisation of the hydrogen in the Universe by taking into account secondary collisional ionisations from the X-ray radiation. They extrapolate to high redshift () the Hiroi et al. (2012) hard X-ray luminosity function and conclude that AGN only contribute a small fraction of the photon rate densities required to ionise the Universe at these redshifts. Our results are in agreement with the above studies and do not support claims for a dominant role of AGN to the ionisation of the hydrogen in the Universe at redshift (Glikman et al., 2011; Giallongo et al., 2015). This discrepancy is likely related to the way different groups select their samples and subsequently account for this selection in the analysis. It also highlights the need for further research to better constrain the impact of AGN radiation to the ionisation of the Universe. Glikman et al. (2011) estimate type-I QSO space densities at that are a factor of 3-4 higher than those determined by Masters et al. (2012) or Ikeda et al. (2011) at similar redshifts and luminosities. Giallongo et al. (2015) combined X-ray and HST optical/near-IR data in the CANDELS GOODS-S region to identify optical sources with photometric or spectroscopic redshifts and then study their X-ray properties following methods described by Fiore et al. (2012). Their approach allows them to identify faint AGN with X-ray luminosities as a low as in the redshift interval . Strictly speaking the Giallongo et al. (2015) photon-rate densities in Figure 14 are upper limits. The UV photon escape fraction is set to one, part of the observed UV radiation of individual sources may be associated with the AGN host galaxy, the photometric redshift uncertainties of at least some sources in the sample are large. The increase of the X-ray depth in the CANDELS GOODS-S from 4 to 7 Ms (PI Brandt) region will help better constrain the faint-end of the AGN luminosity function at high redshift and their role in the ionisation of the Universe.

Finally the parametric X-ray luminosity functions derived in this paper are used to make predictions on the number of AGN that eROSITA (Merloni et al., 2012) surveys will detect. Our LDDE parametrisation is convolved with the expected X-ray sensitivity of the 4-year eROSITA All Sky Survey. The number of AGN in logarithmic luminosity bins of size is plotted as a function of 2-10 keV luminosity in Figure 15. Predictions are presented for 3 redshift intervals, , and . This shows that surveys by eROSITA will provide tight constraints on the evolution of bright AGN and will allow us to explore with high statistical significance the evidence for luminosity dependence of the AGN population at high redshift. This however, would also require a dedicate follow-up program to identify high redshift AGN among the eROSITA population. High multiplex spectroscopic facilities that are able to simultaneously observe large number of targets over a wide field of view are essential for eROSITA X-ray source follow-ups. The SDSS/BOSS spectrographs (Smee et al., 2013) at the Apache Point SDSS telescope (Gunn et al., 2006) and in the future the ESO/4MOST facility (4-metre Multi-Object Spectroscopic Telescope, de Jong et al., 2014) are well suited for follow-up observations of the eROSITA sky.

6 Conclusions

X-ray data from Chandra deep surveys and the shallow/wide XMM-XXL sample are combined to explore the evolution of the X-ray luminosity function at high redshift, . Our analysis accounts for Poisson errors in the calculation of fluxes and luminosities as well as photometric redshift uncertainties. We also show that the latter point is crucial for unbiased AGN space density measurements. The sample used in the paper consists of nearly 340 sources with either photometric (212) or spectroscopic (128) redshift in the redshift range . The luminosity baseline of the sample is at . Our main findings are

• the AGN population evolves strongly between the redshift intervals and

• there is evidence, significant at the 90% level, that the amplitude of the AGN evolution depends on X-ray luminosity. Sources with luminosities appear to evolve faster than brighter ones.

• the faint-end slope of UV/optical QSO luminosity functions is steeper than that of the X-ray selected AGN samples. This implies an increasing fraction of type-I AGN with decreasing X-ray luminosity at .

• X-ray AGN may dominate or at least contribute substantially to the UV photon rate density required to keep the Universe ionised to . At higher redshift the contribution of AGN to UV hydrogen ionising field decreases.

7 Acknowledgements

This work has made use of the Rainbow Cosmological Surveys Database, which is operated by the Universidad Complutense de Madrid (UCM), partnered with the University of California Observatories at Santa Cruz (UCO/Lick,UCSC). 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”. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Footnotes

1. In the case of a step function (i.e. non-parametric) model for and vanishing uncertainties for the redshifts and luminosities of individual sources, it is straightforward to show that the maximum likelihood value of in equation 1 reduces to the Page & Carrera (2000) binned luminosity function estimator.
2. http://www.bo.astro.it/~gilli/counts.html

References

1. Aihara H., et al., 2011, ApJS, 193, 29
2. Aird J., Coil A. L., Georgakakis A., Nandra K., Barro G., Perez-Gonzalez P. G., 2015, ArXiv e-prints 1503.01120
3. Aird J., Coil A. L., Moustakas J., Blanton M. R., Burles S. M., Cool R. J., Eisenstein D. J., Smith M. S. M., Wong K. C., Zhu G., 2012, ApJ, 746, 90
4. Aird J., et al., 2010, MNRAS, 401, 2531
5. Akylas A., Georgantopoulos I., Georgakakis A., Kitsionas S., Hatziminaoglou E., 2006, A&A, 459, 693
6. Alexander D. M., et al., 2003, AJ, 126, 539
7. Assef R. J., Eisenhardt P. R. M., Stern D., Tsai C.-W., Wu J., Wylezalek D., Blain A. W., Bridge C. R., Donoso E., Gonzales A., Griffith R. L., Jarrett T. H., 2014, ArXiv e-prints 1408.1092
8. Assef R. J., Stern D., Kochanek C. S., Blain A. W., Brodwin M., Brown M. J. I., Donoso E., Eisenhardt P. R. M., Jannuzi B. T., Jarrett T. H., Stanford S. A., Tsai C.-W., Wu J., Yan L., 2013, ApJ, 772, 26
9. Barger A. J., Cowie L. L., Capak P., Alexander D. M., Bauer F. E., Brandt W. N., Garmire G. P., Hornschemeier A. E., 2003a, ApJ, 584, L61
10. Barger A. J., Cowie L. L., Capak P., Alexander D. M., Bauer F. E., Fernandez E., Brandt W. N., Garmire G. P., Hornschemeier A. E., 2003b, AJ, 126, 632
11. Bongiorno A., et al., 2012, MNRAS, 427, 3103
12. Brandt W. N., Alexander D. M., 2015, A&A Rev., 23, 1
13. Brightman M., Nandra K., 2011a, MNRAS, 413, 1206
14. —, 2011b, MNRAS, 414, 3084
15. Brusa M., et al., 2009, A&A, 507, 1277
16. Buchner J., Georgakakis A., Nandra K., Brightman M., Menzel M.-L., Liu Z., Hsu L.-T., Salvato M., Rangel C., Aird J., Merloni A., Ross N., 2015, ArXiv e-prints 1501.02805
17. Buchner J., Georgakakis A., Nandra K., Hsu L., Rangel C., Brightman M., Merloni A., Salvato M., Donley J., Kocevski D., 2014, A&A, 564, A125
18. Budavári T., Szalay A. S., 2008, ApJ, 679, 301
19. Burlon D., Ajello M., Greiner J., Comastri A., Merloni A., Gehrels N., 2011, ApJ, 728, 58
20. Cardamone C. N., Urry C. M., Schawinski K., Treister E., Brammer G., Gawiser E., 2010, ApJ, 721, L38
21. Civano F., et al., 2010, ApJ, 717, 209
22. —, 2011, ApJ, 741, 91
23. —, 2012, ApJS, 201, 30
24. Clerc N., Adami C., Lieu M., Maughan B., Pacaud F., Pierre M., Sadibekova T., Smith G. P., Valageas P., Altieri B., Benoist C., Maurogordato S., Willis J. P., 2014, MNRAS, 444, 2723
25. Cowie L. L., Barger A. J., Bautz M. W., Brandt W. N., Garmire G. P., 2003, ApJ, 584, L57
26. Croton D. J., et al., 2006, MNRAS, 365, 11
27. Davis M., et al., 2007, ApJ, 660, L1
28. Dawson K. S., et al., 2013, AJ, 145, 10
29. de Jong R. S., et al., 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 0
30. Del Moro A., Watson M. G., Mateos S., Akiyama M., Hashimoto Y., Tamura N., Ohta K., Carrera F. J., Stewart G., 2009, A&A, 493, 445
31. Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
32. Ebrero J., Carrera F. J., Page M. J., Silverman J. D., Barcons X., Ceballos M. T., Corral A., Della Ceca R., Watson M. G., 2009, A&A, 493, 55
33. Elvis M., et al., 2009, ApJS, 184, 158
34. Fabian A. C., 1999, MNRAS, 308, L39
35. Feroz F., Hobson M. P., 2008, MNRAS, 384, 449
36. Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
37. Ferrarese L., Merritt D., 2000, ApJ, 539, L9
38. Fiore F., et al., 2012, A&A, 537, A16
39. Fontanot F., Cristiani S., Monaco P., Nonino M., Vanzella E., Brandt W. N., Grazian A., Mao J., 2007, A&A, 461, 39
40. Gallagher S. C., Brandt W. N., Chartas G., Priddey R., Garmire G. P., Sambruna R. M., 2006, ApJ, 644, 709
41. Gebhardt K., Bender R., Bower G., Dressler A., Faber S. M., Filippenko A. V., Green R., Grillmair C., Ho L. C., Kormendy J., Lauer T. R., Magorrian J., Pinkney J., Richstone D., Tremaine S., 2000, ApJ, 539, L13
42. Georgakakis A., Nandra K., 2011, ArXiv: 1101.4943
43. Georgakakis A., Nandra K., Laird E. S., Aird J., Trichas M., 2008, MNRAS, 388, 1205
44. Georgakakis A., et al., 2009, MNRAS, 397, 623
45. —, 2011, MNRAS, 418, 2590
46. —, 2014, MNRAS, 440, 339
47. Giallongo E., et al., 2015, ArXiv e-prints 1502.02562
48. Gibson R. R., Jiang L., Brandt W. N., Hall P. B., Shen Y., Wu J., Anderson S. F., Schneider D. P., Vanden Berk D., Gallagher S. C., Fan X., York D. G., 2009, ApJ, 692, 758
49. Gilli R., Comastri A., Hasinger G., 2007, A&A, 463, 79
50. Glikman E., Djorgovski S. G., Stern D., Dey A., Jannuzi B. T., Lee K.-S., 2011, ApJ, 728, L26
51. Graham A. W., 2012, ApJ, 746, 113
52. Graham A. W., Erwin P., Caon N., Trujillo I., 2001, ApJ, 563, L11
53. Graham A. W., Onken C. A., Athanassoula E., Combes F., 2011, MNRAS, 412, 2211
54. Grissom R. L., Ballantyne D. R., Wise J. H., 2014, A&A, 561, A90
55. Grogin N. A., et al., 2011, ApJS, 197, 35
56. Gültekin K., et al., 2009, ApJ, 698, 198
57. Gunn J. E., et al., 2006, AJ, 131, 2332
58. Guo Y., et al., 2013, ApJS, 207, 24
59. Haardt F., Madau P., 2012, ApJ, 746, 125
60. Häring N., Rix H.-W., 2004, ApJ, 604, L89
61. Hasinger G., 2008, A&A, 490, 905
62. Hasinger G., Miyaji T., Schmidt M., 2005, A&A, 441, 417
63. Hewett P. C., Foltz C. B., 2003, AJ, 125, 1784
64. Hiroi K., Ueda Y., Akiyama M., Watson M. G., 2012, ApJ, 758, 49
65. Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 654, 731
66. Hsieh B.-C., Wang W.-H., Hsieh C.-C., Lin L., Yan H., Lim J., Ho P. T. P., 2012, ApJS, 203, 23
67. Hsu L.-T., et al., 2014, ApJ, 796, 60
68. Ikeda H., et al., 2011, ApJ, 728, L25
69. Jahnke K., Macciò A. V., 2011, ApJ, 734, 92
70. Jeffreys H., 1961, Theory of Probability (3rd Edition). Oxford University Press, New York
71. Just D. W., Brandt W. N., Shemmer O., Steffen A. T., Schneider D. P., Chartas G., Garmire G. P., 2007, ApJ, 665, 1004
72. Kalfountzou E., Civano F., Elvis M., Trichas M., Green P., 2014, MNRAS, 445, 1430
73. King A., 2003, ApJ, 596, L27
74. —, 2005, ApJ, 635, L121
75. Koekemoer A. M., Alexander D. M., Bauer F. E., Bergeron J., Brandt W. N., Chatzichristou E., Cristiani S., Fall S. M., Grogin N. A., Livio M., Mainieri V., Moustakas L. A., Padovani P., Rosati P., Schreier E. J., Urry C. M., 2004, ApJ, 600, L123
76. Koekemoer A. M., et al., 2011, ApJS, 197, 36
77. Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
78. Laird E. S., et al., 2009, ApJS, 180, 102
79. Lehmer B. D., et al., 2005, ApJS, 161, 21
80. Lilly S. J., et al., 2009, ApJS, 184, 218
81. Liu Z., et al., 2015, MNRAS submitted
82. Loredo T. J., 2004, in American Institute of Physics Conference Series, Vol. 735, American Institute of Physics Conference Series, Fischer R., Preuss R., Toussaint U. V., eds., pp. 195–206
83. Luo B., et al., 2014, ApJ, 794, 70
84. Lusso E., Worseck G., Hennawi J. F., Prochaska J. X., Vignali C., Stern J., O’Meara J. M., 2015, ArXiv e-prints 1503.02075
85. Lusso E., et al., 2010, A&A, 512, A34
86. Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648
87. Magorrian J., Tremaine S., Richstone D., Bender R., Bower G., Dressler A., Faber S. M., Gebhardt K., Green R., Grillmair C., Kormendy J., Lauer T., 1998, AJ, 115, 2285
88. Marconi A., Hunt L. K., 2003, ApJ, 589, L21
89. Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169
90. Marshall H. L., Tananbaum H., Avni Y., Zamorani G., 1983, ApJ, 269, 35
91. Masters D., et al., 2012, ApJ, 755, 169
92. McGreer I. D., Jiang L., Fan X., Richards G. T., Strauss M. A., Ross N. P., White M., Shen Y., Schneider D. P., Myers A. D., Brandt W. N., DeGraf C., Glikman E., Ge J., Streblyanska A., 2013, ApJ, 768, 105
93. McLure R. J., Dunlop J. S., 2002, MNRAS, 331, 795
94. Menzel M.-L., et al., 2015, MNRAS submitted
95. Merloni A., et al., 2012, ArXiv e-prints1209.3114
96. —, 2014, MNRAS, 437, 3550
97. Miyaji T., Hasinger G., Schmidt M., 2001, A&A, 369, 49
98. Miyaji T., et al., 2015, ArXiv e-prints 1503.00056
99. Nandra K., O’Neill P. M., George I. M., Reeves J. N., 2007, MNRAS, 382, 194
100. Nandra K., Pounds K. A., 1994, MNRAS, 268, 405
101. Nandra K., et al., 2015, arXiv 1503.09078
102. Page M. J., Carrera F. J., 2000, MNRAS, 311, 433
103. Palanque-Delabrouille N., Yeche C., Myers A. D., Petitjean P., Ross N. P., Sheldon E., Aubourg E., Delubac T., Le Goff J.-M., Pâris I., Rich J., Dawson K. S., Schneider D. P., Weaver B. A., 2011, A&A, 530, A122
104. Pawlik A. H., Schaye J., van Scherpenzeel E., 2009, MNRAS, 394, 1812
105. Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The art of scientific computing. Cambridge: University Press, —c1992, 2nd ed.
106. Rangel C., Nandra K., Laird E. S., Orange P., 2013, MNRAS, 428, 3089
107. Reichard T. A., Richards G. T., Schneider D. P., Hall P. B., Tolea A., Krolik J. H., Tsvetanov Z., Vanden Berk D. E., York D. G., Knapp G. R., Gunn J. E., Brinkmann J., 2003, AJ, 125, 1711
108. Richards G. T., Myers A. D., Gray A. G., Riegel R. N., Nichol R. C., Brunner R. J., Szalay A. S., Schneider D. P., Anderson S. F., 2009, ApJS, 180, 67
109. Richards G. T., et al., 2006, AJ, 131, 2766
110. Rodighiero G., Cimatti A., Franceschini A., Brusa M., Fritz J., Bolzonella M., 2007, A&A, 470, 21
111. Ross N. P., et al., 2012, ApJS, 199, 3
112. —, 2013, ApJ, 773, 14
113. Salvato M., et al., 2009, ApJ, 690, 1250
114. —, 2011, ApJ, 742, 61
115. Savorgnan G., Graham A. W., Marconi A., Sani E., Hunt L. K., Vika M., Driver S. P., 2013, MNRAS, 434, 387
116. Schaerer D., Hempel A., Egami E., Pelló R., Richard J., Le Borgne J.-F., Kneib J.-P., Wise M., Boone F., 2007, A&A, 476, 97
117. Schmidt M., Schneider D. P., Gunn J. E., 1995, AJ, 110, 68
118. Shemmer O., Brandt W. N., Vignali C., Schneider D. P., Fan X., Richards G. T., Strauss M. A., 2005, ApJ, 630, 729
119. Silk J., Rees M. J., 1998, A&A, 331, L1
120. Simpson C., 2005, MNRAS, 360, 565
121. Smee S. A., et al., 2013, AJ, 146, 32
122. Soltan A., 1982, MNRAS, 200, 115
123. Stalin C. S., Petitjean P., Srianand R., Fox A. J., Coppolani F., Schwope A., 2010, MNRAS, 401, 294
124. Steffen A. T., Strateva I., Brandt W. N., Alexander D. M., Koekemoer A. M., Lehmer B. D., Schneider D. P., Vignali C., 2006, AJ, 131, 2826
125. Stern D., et al., 2014, ApJ, 794, 102
126. Strüder L., et al., 2001, A&A, 365, L18
127. Sutherland W., Saunders W., 1992, MNRAS, 259, 413
128. Telfer R. C., Zheng W., Kriss G. A., Davidsen A. F., 2002, ApJ, 565, 773
129. Trump J. R., et al., 2009, ApJ, 696, 1195
130. Turner M. J. L., et al., 2001, A&A, 365, L27
131. Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104
132. Ueda Y., Akiyama M., Ohta K., Miyaji T., 2003, ApJ, 598, 886
133. Vanden Berk D. E., et al., 2001, AJ, 122, 549
134. Vignali C., Brandt W. N., Schneider D. P., Kaspi S., 2005, AJ, 129, 2519
135. Vito F., Gilli R., Vignali C., Comastri A., Brusa M., Cappelluti N., Iwasawa K., 2014, MNRAS, 445, 3557
136. Vito F., Vignali C., Gilli R., Comastri A., Iwasawa K., Brandt W. N., Alexander D. M., Brusa M., Lehmer B., Bauer F. E., Schneider D. P., Xue Y. Q., Luo B., 2013, MNRAS, 428, 354
137. Weigel A. K., Schawinski K., Treister E., Urry C. M., Koss M., Trakhtenbrot B., 2015, MNRAS, 448, 3167
138. Xue Y. Q., et al., 2011, ApJS, 195, 10
139. Yencho B., Barger A. J., Trouille L., Winter L. M., 2009, ApJ, 698, 380