High-redshift AGN in the Chandra Deep Fields: the obscured fraction and space density of the sub-L_* population Free

Abstract

We investigate the population of high-redshift (3 ≤ z < 6) active galactic nuclei (AGN) selected in the two deepest X-ray surveys, the 7 Ms Chandra Deep Field-South and 2 Ms Chandra Deep Field-North. Their outstanding sensitivity and spectral characterization of faint sources allow us to focus on the sub-L_* regime (logL_X ≲ 44), poorly sampled by previous works using shallower data, and the obscured population. Taking fully into account the individual photometric-redshift probability distribution functions, the final sample consists of ≈102 X-ray-selected AGN at 3 ≤ z < 6. The fraction of AGN obscured by column densities logN_H > 23 is ∼0.6–0.8, once incompleteness effects are taken into account, with no strong dependence on redshift or luminosity. We derived the high-redshift AGN number counts down to F_{0.5–2 keV} = 7 × 10⁻¹⁸ erg cm⁻² s⁻¹, extending previous results to fainter fluxes, especially at z > 4. We put the tightest constraints to date on the low-luminosity end of AGN luminosity function at high redshift. The space density, in particular, declines at z > 3 at all luminosities, with only a marginally steeper slope for low-luminosity AGN. By comparing the evolution of the AGN and galaxy densities, we suggest that such a decline at high luminosities is mainly driven by the underlying galaxy population, while at low luminosities there are hints of an intrinsic evolution of the parameters driving nuclear activity. Also, the black hole accretion rate density and star formation rate density, which are usually found to evolve similarly at z ≲ 3, appear to diverge at higher redshifts.

1 INTRODUCTION

Supermassive black holes (SMBHs) and their hosting galaxies are broadly recognized to influence the evolution of each other over cosmic time. This ‘co-evolution’ is reflected by the tight relations between the masses of SMBH and the properties of host galaxies in the nearby Universe, such as masses and velocity dispersions of the bulges (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Marconi & Hunt 2003) and the broadly similar evolution of the star formation and black hole accretion densities in the last ∼10 Gyr (e.g. Aird et al. 2015), although the details of this interplay are still not well known (see e.g. Kormendy & Ho 2013, and references therein). Studying galaxies and SMBH in the early Universe, where these relations could be set in place, would boost our knowledge of how SMBH and galaxies formed and evolved. However, while galaxy properties have been traced as far back in time as z ∼ 8–10 (e.g. Bouwens et al. 2015), our knowledge of SMBH is limited to later times.

Only ∼90 accreting SMBHs, shining as active galactic nuclei (AGN), have been identified at z > 6 (e.g. Bañados et al. 2016; Wang et al. 2017) and are usually found to have masses of the order of 1–10 billion solar masses (e.g. Mortlock et al. 2011; Wu et al. 2015). The presence of such massive black holes a few 10⁸ yr after the big bang challenges our understanding of SMBH formation and growth in the early Universe, one of the major issues in modern astrophysics (e.g. Reines & Comastri 2016, and references therein). Different classes of theories have been proposed to explain the formation of the BH seeds that eventually became SMBH, the two most promising ones involving ‘light seeds’ (M ∼ 10² M_⊙), as remnants of the first Pop III stars, and ‘heavy seeds’ (M ∼ 10⁴⁻⁶ M_⊙), perhaps formed during the direct collapse of giant pristine gas clouds (e.g. Haiman 2013; Johnson & Haardt 2016; Volonteri et al. 2016b, and references therein). To match the masses of SMBH discovered at z > 6, all such models require continuous nearly Eddington-limited or even super-Eddington accretion phases during which the growing SMBH is plausibly buried in material with large column densities, even exceeding the Compton-thick level (e.g. Pacucci et al. 2015). However, these objects represent the extreme tail of the underlying distribution (in terms of both mass and luminosity) and are not representative of the overall population.

X-ray surveys are the most suitable tools for investigating the evolution of the bulk of the AGN population up to high redshift: being less affected by absorption and galaxy dilution, they provide cleaner and more complete AGN identification with respect to optical/IR surveys (Brandt & Alexander 2015, and references therein). Over the last two decades, several works have focused on the properties and evolution of X-ray selected, z > 3 AGN in wide (e.g. Brusa et al. 2009; Civano et al. 2011; Hiroi et al. 2012; Marchesi et al. 2016) and deep (e.g. Vignali et al. 2002; Fiore et al. 2012; Vito et al. 2013; Giallongo et al. 2015; Weigel et al. 2015; Cappelluti et al. 2016) surveys performed with Chandra and XMM–Newton, or using combinations of different surveys (e.g. Kalfountzou et al. 2014; Vito et al. 2014; Georgakakis et al. 2015). Common findings amongst such works are (1) a decline of the space density of luminous (logL_X ≳ 44) AGN proportional to (1 + z)^d with d ∼ −6 (similar to the exponential decline of the space density of optically selected quasars, e.g. McGreer et al. 2013), and (2) a larger fraction of obscured AGN than that usually derived at lower redshifts, particularly at moderate-to-high luminosities (e.g. Aird et al. 2015; Buchner et al. 2015).

However, most of the low-luminosity (logL_X ≲ 43) and z ≳ 4 AGN are missed even by the deepest surveys, leading to discrepant results amongst different studies. For instance, the evolution of the space density of low-luminosity, X-ray-detected AGN is largely unconstrained, while Georgakakis et al. (2015) reported an apparent strong flattening of the faint end of the AGN X-ray luminosity function (XLF) at z > 3, Vito et al. (2014) found that the decline of the space density of low-luminosity AGN is consistent with that of AGN with higher luminosities. Moreover, Giallongo et al. (2015), using detection techniques that search for clustering of photons in energy, space and time, reported the detection of several faint AGN, resulting in a very steep XLF faint end (see also Fiore et al. 2012). These results also have strong impact on the determination of the AGN contribution to cosmic reionization (e.g. Madau & Haardt 2015). Moreover, the typical obscuration levels in these faint sources remain unknown, although hints of a decrease of the obscured AGN fraction with decreasing luminosity (for logL_X ≲ 44) at high-redshift have been found (e.g. Aird et al. 2015; Buchner et al. 2015; Georgakakis et al. 2015). This relation is the opposite trend to that found at low redshift (e.g. Aird et al. 2015; Buchner et al. 2015), where the obscured AGN fraction shows a clear anticorrelation with AGN luminosity. Finally, the very detection of faint z > 5 AGN in deep X-ray surveys is debated (e.g. Vignali et al. 2002; Giallongo et al. 2015; Weigel et al. 2015; Cappelluti et al. 2016; Parsa et al. 2017).

The recently completed 7 Ms Chandra Deep Field-South (CDF-S; Luo et al. 2017) observations provide the deepest X-ray view of the early Universe, reaching a flux limit of F_{0.5−2 keV} = 6.4 × 10⁻¹⁸ erg cm⁻² s⁻¹. Moreover, the catalogue of X-ray sources in the second deepest X-ray survey to date (limiting flux F_{0.5−2 keV} = 1.2 × 10⁻¹⁷ erg cm⁻² s⁻¹), the 2 Ms Chandra Deep Field-North (CDF-N; Alexander et al. 2003), was recently re-analysed by Xue et al. (2016) with the same detection procedure applied to the CDF-S, which provides detections for more real sources. Therefore, the two deepest Chandra fields allow us now to study high-redshift, faint AGN using homogeneous data sets. In Vito et al. (2016), we applied a stacking technique to CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) selected galaxies to study the X-ray emission from individually undetected sources in the 7 Ms CDF-S, finding that the emission is probably mostly due to X-ray binaries (XRBs) rather than nuclear accretion, and concluding that most of the SMBH growth at 3.5 < z < 6.5 occurred during bright AGN phases. In this paper, we combine the 7 Ms CDF-S and 2 Ms CDF-N data to study the X-ray properties of X-ray-detected AGN at z > 3, with a particular focus on low-luminosity sources (logL_X ≲ 44), which are best sampled by deep, pencil-beam surveys. Taking fully into account the probability distribution functions (PDFs) of the photometric redshifts for sources lacking spectroscopic identifications (see e.g. Marchesi et al. 2016 for a similar use of the photometric redshifts), the final sample consists of ≈102 X-ray-detected AGN at 3 ≤ z < 6. The number of sources contributing to this sample with their PDF(z) is 118. We performed a spectral analysis on our sample assuming the X-ray spectra are well represented by power-law emission subjected to Galactic and intrinsic absorption. The spectral analysis allowed us to take into account the full probability distribution of the intrinsic column densities. We also considered the probability distribution of the count rates of X-ray-detected sources and applied a correction to mitigate the Eddington bias. The flux (and hence luminosity) probability distributions were derived by applying for each source the proper response matrices and the conversion factors between count rate and flux, which depend on the observed spectral shape. We present the trends of the obscured AGN fraction with redshift and luminosity, the number counts, and the space-density evolution of 3 < z < 6 AGN. Throughout this paper, we will use a H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and |$\Omega _\Lambda =0.7$| cosmology and we will assume Galactic column densities of logN_H = 0.9 × 10²⁰ and 1.6 × 10²⁰ cm⁻² along the line of sight of CDF-S and CDF-N, respectively. Errors and upper limits are quoted at the 68 per cent confidence level, unless otherwise noted.

2 THE SAMPLE

2.1 AGN parent sample and redshifts in the 7 Ms CDF-S

We selected a sample of X-ray detected, z > 3 AGN in the 7 Ms CDF-S,¹ the deepest X-ray survey to date, from the Luo et al. (2017, hereafter L17) catalogue, which also provides multiwavelength identifications and spectroscopic and photometric redshifts for the X-ray sources. In particular, photometric redshifts were collected from Luo et al. (2010), Rafferty et al. (2011), Hsu et al. (2014, hereafter H14), Skelton et al. (2014, hereafter S14), Santini et al. (2015) and Straatman et al. (2016, hereafter S16). Each X-ray source can therefore be associated with up to six different photometric redshifts. We considered only sources located in the area (∼330 arcmin², red region in the left-hand panel of Fig. 1) of the survey where the effective exposure is ≥1 Ms, in order to exclude the outskirts of the field, where the point spread function (PSF) distortions and the effects of the vignetting affect the quality of the X-ray data and the optical identification rate and accuracy. Moreover, the inner region of the CDF-S is covered by the deepest optical/IR observations (green region), which are essential to derive highly reliable spectroscopic and photometric redshifts. With this selection, our parent sample in the CDF-S consists of 952 out of the 1008 X-ray sources in L17.

We adopted the L17 definitions for the spectroscopic redshift quality. Following L17, we associated with each X-ray source a spectroscopic redshift if it is defined as ‘secure’ or ‘insecure’ but in agreement within |$\frac{|z_{\mathrm{phot}}-z_{\mathrm{spec}}|}{1+z_{\mathrm{spec}}}<15{\rm per \, cent}$| with at least one photometric redshift. Using more conservative criteria, such as requiring that the spectroscopic redshift agrees with at least 2 or > 50 per cent of the available photometric redshifts, would have no effect on the final sample of z > 3 AGN. The photometric redshifts used to validate the ‘insecure’ spectroscopic redshifts at z > 3 are of good quality, with a 68 per cent confidence level Δz = 0.03–0.2.

We also estimated the accuracy of the confidence intervals provided by the PDFs by computing the fraction of sources whose spectroscopic redshift is included in the 68 per cent confidence interval provided by its PDF (defined as the narrowest redshift interval where the integrated redshift probability is 0.68). If the PDFs provided accurate confidence intervals, that fraction would be 0.68, while we found 0.49, 0.50 and 0.63 for H14, S14 and S16, respectively, reflecting a general mild underestimate of the confidence intervals, hence of the photometric-redshift errors. This effect could be due to underestimating the errors of the fitted photometric data (e.g. see section 5.3 in Yang et al. 2014). We found indeed that the most accurate confidence intervals are provided by S16, who addressed in detail this issue by employing an empirical technique to derive more accurate photometric errors than those usually provided by detection software such as SExtractor. The reported fractions refer to the particular comparative spectroscopic sample, i.e. X-ray-selected galaxies, and are expected to be different considering the entire galaxy samples in those works. The PDFs are usually derived by fitting the observed spectral energy distribution (SED) with models M of galactic emission varying the redshift as |$\mathrm{PDF}(z)\propto \mathrm{\mathrm{exp}}(-0.5\chi ^2(z))$|⁠, where |$\chi ^2(z)=\sum _i\frac{(M_i(z)-\text{SED}_i)^2}{\sigma _i^2}$| is the test statistic of the fit, and the index i represents the different photometric bands. If the photometric errors σ_i are underestimated, the resulting PDFs will be too sharp and their confidence intervals will be underestimated as well. In this case, more accurate confidence intervals can be obtained by multiplying the photometric errors by a factor of α, which represents the average underestimating factor of the photometric errors amongst the used bands, or, equivalently, by using the ‘corrected’ distribution |$\mathrm{PDF^{corr}}(z)=\mathrm{PDF^{input}}(z)^{\frac{1}{\alpha ^2}}$|⁠, where α² is computed empirically such that the 68 per cent considered interval provided by their PDFs encompasses the associated spectroscopic redshift in 68 per cent of the sample. This procedure is equivalent to empirically ‘correcting’ (i.e. increasing) the photometric errors following σ^corr = α × σ^input. We obtained α² = 5.2, 4.4 and 1.5 for H14, S14 and S16, respectively, and use the ‘corrected’ PDFs hereafter.

All of these tests led us to adopt the photometric redshift from S16 as first choice. The photometric redshifts from H14 and S14 have similar accuracy, but H14 provide the redshifts for the entire Extended CDF-S (E-CDF-S), while S14 (as well as S16) is limited to the GOODS-S/CANDELS region. We therefore used the photometric redshifts from H14 and S14 as second and third choices, respectively.² Amongst the 21 remaining sources with no entries in the above-considered photometric-redshift catalogues, two sources have a photometric redshift from Rafferty et al. (2011). We approximated their PDF(z), not provided by that work, as normalized Gaussian functions centred on the nominal redshift and with σ equal to the 1σ error, separately for the positive and negative sides. Given the very low number of sources for which we adopted the Rafferty et al. (2011) redshifts (two, and only one will be included in the final sample), the approximation we used to describe their PDF(z) does not have any effect on the final results.

Nineteen sources out of the 952 X-ray sources in the parent sample (1.9 per cent) remain with no redshift information. Most of them (14/19) are not associated with any counterpart from catalogues at different wavelengths. The five sources with an optical counterpart but no redshift information have been assigned a flat PDF. We preferred not to use redshift priors based on optical magnitudes, as these sources are part of a very peculiar class of objects (i.e. extremely X-ray faint sources), whose redshift distribution is not well represented by any magnitude-based distribution derived for different classes of galaxies. We did not include the 14 X-ray sources with no optical counterparts in this work, as a significant number of them are expected to be spurious detections. In fact, L17 estimated the number of spurious detections in the entire 7 Ms CDF-S main catalogue to be ∼19 and argued that, given the superb multiwavelength coverage of the CDF-S and sharp Chandra PSF, resulting in high-confidence multiwavelength identifications, X-ray sources with multiwavelength counterparts are extremely likely to be real detections. Therefore, most of the spurious detections are expected to be accounted for by X-ray sources with no optical counterpart. This conclusion is especially true considering that most of the unmatched sources lie in the CANDELS/GOODS-S field, where the extremely deep multiwavelength observations would likely have detected their optical counterpart if they were true sources. We also checked their binomial no-source probability (P_B, i.e. the probability that the detected counts of a source are due to a background fluctuation), provided by L17 for all the sources in their catalogue. Thirteen out of the 14 sources with no counterparts have P_B ≳ 10⁻⁴, close to the detection threshold of P_B = 0.007 (L17). For comparison, ∼ 80 per cent of the sources with optical counterparts have P_B < 10⁻⁴. The notable exception is XID 912, which has no optical counterpart, but P_B ∼ 10⁻¹³ and ≈77 net counts. L17 suggested that XID 912 is a off-nuclear X-ray source associated with a nearby low-redshift galaxy, detected in optical/IR observations.

2.2 AGN parent sample and redshifts in the 2 Ms CDF-N

Xue et al. (2016, hereafter X16) presented an updated version (683 sources in total) of the original 2 Ms CDF-N catalogue (Alexander et al. 2003), applying the same detection procedure used in the CDF-S (Xue et al. 2011, L17) and providing multiwavelength identifications and spectroscopic and photometric redshifts from the literature. In particular, photometric redshifts were collected from S14 and Yang et al. (2014, hereafter Y14). Both of these works provide the PDF of the photometric redshifts. We adopted the spectroscopic redshifts collected from X16, as they considered only those marked as ‘secure’ in the original works.

For X-ray sources lacking spectroscopic redshift, we followed the procedure described in Section 2.1 to define a priority order amongst the two used photometric catalogues. In particular, we derived median values of |Δz|/(1 + z_spec) = 0.026 both for S14 and Y14, σ_NMAD = 0.025 and 0.035 and α² = 9.2 and 3.0 for S14 and Y14, respectively. The α² values mean that the photometric redshift PDFs from Y14 account for the redshift uncertainties more realistically than the S14 ones: this behaviour can be explained again by considering that Y14 applied an empirical method to estimate the photometric errors more reliably than those provided by standard detection software. We therefore adopted the photometric redshifts from Y14 and S14 as first and second choice, respectively.

As in Section 2.1, we considered only the sources in the area covered by ≥1 Ms effective exposure (∼215 arcmin²; red contour in the right-hand panel of Fig. 1, almost coincident with the CANDELS survey in that field) as the parent sample (485 sources). Amongst them, 35 sources ( ∼ 7 per cent) have no redshift information. This non-identification rate is mostly due to the method used in X16 to match the X-ray sources with the entries in the photometric-redshift catalogues. First, for each X-ray source they identified a primary multiwavelength counterpart using a likelihood-ratio procedure (Luo et al. 2010). They then matched the coordinates of the primary counterparts with the photometric-redshift catalogues using a 0.5 arcsec radius. However, the positional uncertainties and offsets amongst the different photometric catalogues can be comparable to the utilized matching radius. We therefore visually inspected the positions of all the X-ray sources on the CANDELS/GOODS-N images³ (Grogin et al. 2011; Koekemoer et al. 2011) and could match most of the X-ray sources with no redshift information in X16 with an entry in one of the considered photometric-redshift catalogues. The visual multiwavelength matching allowed us also to associate a spectroscopic redshift from Y14 (which included only high-quality spectroscopic redshifts) with several sources with only photometric or no redshifts in X16.

The resulting number of X-ray sources with no redshift information is 12 ( ∼ 2.5 per cent). As in Section 2, we excluded sources with no multiwavelength counterpart (8 sources) from this analysis, as most of them are expected to be spurious. Two out of the remaining four sources with counterparts but no redshifts in the considered catalogues have photometric-redshift entries in the CANDELS/GOODS-N catalogue (Kodra et al. in preparation). Their PDF(z) lies entirely at z < 3, thus we excluded these sources from the high-redshift sample. Finally, we assigned flat PDF(z) over the range of z = 0–10 to the only two remaining sources with multiwavelength counterparts but no redshifts. Table 1 summarizes the number of X-ray sources in the CDF-N that are associated with spectroscopic, photometric or no redshifts.

2.3 The sample of high-redshift AGN in the Chandra Deep Fields

We checked the photometric redshift accuracy by plotting in Fig. 2 the Δz/(1 + z_spec) and the σ_NMAD for sources with secure spectroscopic redshifts (see also Table 1). The σ_NMAD is computed in a shifting interval of redshift with variable width such to include 10 sources (separately for the positive and negative sides). The photometric redshift for each AGN is chosen following the priority order described in Sections 2.1 and 2.2. We considered only sources within the area covered by ≥1 Ms exposure. The scatter increases slightly at z ∼ 2.5 (by a factor of ∼2–3), but the photometric redshift accuracy does not appear to deteriorate dramatically at high redshift.

Fig. 3 presents the redshift distributions of the sources at z > 3 in the two deep fields, considering their PDF(z). At z > 6, the source statistics are poor, and sources with no redshift information (red line), which carry little information, constitute a significant fraction of the total number of sources. We therefore chose to limit our investigation to the redshift range of 3 ≤ z < 6. To prevent the inclusion of several sources whose PDFs show long tails of extremely low probability at high redshift, we also required that the probability of a source to be in the considered redshift range (computed as in equation 4) is P(3 ≤ z < 6) > 0.2. The number of sources satisfying this criterion is N_hz = 118. Integrating their PDFs in the redshift range of 3 ≤ z < 6, the number of ‘effective’ sources included in the sample is N^eff_hz ≃ 101.6. The cut P(3 ≤ z < 6) > 0.2 results in the rejection of ∼2 effective sources. Table 1 reports the breakdown of the redshift type for sources at 3 ≤ z < 6. Table 2 provides the redshift along with the 68 per cent confidence interval, its origin and the probability of the source to be at 3 ≤ z < 6 according to its PDF for each source in the high-redshift sample.

In principle, equation (4) measures only the uncertainty of the photometric redshift. The probability for a source to be at z₁ ≤ z < z₂ can be computed by weighting equation (4) with the XLF, in order to take into account the intrinsic variation of the AGN space density with redshift, that can introduce systematic errors (e.g. Georgakakis et al. 2015). However, the general high quality of the photometric redshifts we used (see Sections 2.1 and 2.2, and Table 1), due to the availability of deep multiband data in the fields we considered, is reflected in narrow PDF(z). In fact, the PDF(z) of the spectroscopically unidentified sources included in our sample have a median 1σ uncertainty of 0.23. Such narrow PDF(z) would be at most sensitive to the ‘local’ XLF. For continuity reasons, the XLF cannot vary dramatically in a narrow redshift range and the applied weights are therefore expected to be negligible. This statement may not be true in cases of broader PDF(z), which would be therefore sensitive to the space-density variations over larger redshift intervals.

To quantify this effect in our case, we repeated the procedure followed to derive the redshift distribution of the sample, after weighting the PDF(z) with the Ueda et al. (2014) XLF. We also considered additional cases in which the Ueda et al. (2014) XLF is modified at high redshift to match the results of Vito et al. (2014) and Georgakakis et al. (2015). The new redshift distributions are barely distinguishable from Fig. 3, with |$1.7\text{--}5.9{\rm per \, cent}$| (depending on the particular XLF used) of the sample shifting below z = 3. Including in this computation sources with no redshift information (i.e. for which a flat PDF(z) is assumed), we found that |$2.5\text{--}7.7{\rm per \, cent}$| of the samples shift to low redshift. Moreover, the weighting is sensitive in particular to the assumed shape and evolution of the XLF faint end at high redshift, which is currently highly uncertain. For instance, using models with a steeper XLF faint end (e.g. Giallongo et al. 2015) would move a fractional number (<1) of sources from z ≲ 3 to z > 3, slightly increasing our sample size. Given the small effect obtained applying the XLF weighting and the source of additional uncertainty represented by the particular XLF assumed at high redshift, we will use equation (4) to measure the intrinsic redshift probability distribution.

We also made a basic assessment of the effects of the expected presence of redshift outliers in our high-redshift sample (i.e. intrinsically low-redshift sources scattered to high redshift because of the catastrophic failure of their photometric redshifts, and hence PDF(z)). We considered only X-ray sources with a high-quality spectroscopic redshift z_spec < 3, assigned to them a photometric redshift following the priority order described in Sections 2.1 and 2.2, and computed the fraction of redshift outliers, defined as sources for which |z_spec − z_phot|/(1 + z_spec) > 0.15 (see e.g. Luo et al. 2010), that would be counted in the high-redshift sample according to their photometric redshift (i.e. PDF(z)). Considering this fraction and the number of X-ray-selected, spectroscopically unidentified sources with photometric redshift z_phot < 3, we estimated that ≈4 sources in our high-redshift sample ( ≈ 4 per cent) are such low-redshift interlopers. The outlier fraction may be higher for spectroscopically unidentified sources, since these are typically fainter than the control sample we used, but our results are not qualitatively affected even considering an outlier fraction of ≈ 10 per cent.

Our sample includes a few sources studied in detail in previous, dedicated works, like the two Compton-thick candidate AGN in the CDF-S at z = 4.762 (Gilli et al. 2014, ID 539 in Table 2) and z = 3.70 (Comastri et al. 2011, ID 551). We also include the three spectroscopically-identified z > 4 AGN in the CDF-N analysed by Vignali et al. (2002): ID 290 at z = 4.424, ID 330 at z = 5.186, which is also the X-ray source with the highest spectroscopic redshift in the Chandra deep fields, and ID 526 at z = 4.137.⁴

3 SPECTRAL ANALYSIS AND PARAMETER DISTRIBUTIONS

X-ray spectra of the sources in the z > 3 sample were extracted with acis extract (Broos et al. 2010), as described in L17 and X16. Following Vito et al. (2013), we analysed the individual X-ray spectra with xspec⁵ v12.9.0 (Arnaud 1996) assuming an absorbed power-law model with fixed Γ = 1.8 (a typical value for X-ray selected AGN, e.g. Liu et al. 2017), Galactic absorption (Kalberla et al. 2005) and intrinsic absorption (xspec model wabs × zwabs × powerlaw), for a total of two free parameters (the intrinsic column density, N_H, which assumes solar metallicity, and the power-law normalization), in the energy range of 0.5–7 keV. More complex parametrizations are precluded by the general paucity of net counts characterizing the sources in the sample, which, being at high redshift, are typically faint. Appendix A describes our check through spectral simulations that the simplicity of the assumed spectral model does not affect the best-fitting column densities at logN_H ≳ 23 (see Section 5 for discussion). In Appendix B, we study the photon index of a subsample of unobscured sources. The W statistic,⁶ based on the Cash statistic (Cash 1979) and suitable for Poisson data with Poisson background, was used to fit the model to the unbinned spectra. If a spectroscopic redshift is assigned to a source, the fit was performed at that redshift, otherwise we used a grid of redshifts (with step Δz = 0.01), obtaining for each source a set of best-fitting models as a function of redshift. In Section 3.1, we derive the probability distribution of net count rate for each source, including a correction for Eddington bias. In Section 3.2, we convolve the set of best-fitting models with the redshift PDF for each source to derive the probability distribution of best-fitting column density. The fluxes of the sources are derived from the observed count rates in Section 3.3, where we also present the logN-logS of the sample. In Section 3.4, we combine the fluxes and column densities to derive the intrinsic X-ray luminosities. Results from spectral analysis are reported in Table 2.

3.1 Count rates, sky-coverage and Eddington bias

The fundamental quantity measured by an X-ray detector is the count rate of a source, which can then be converted into energy flux assuming a spectral shape and using the detector response matrices. Absorption affects strongly the spectral shape, in particular in the observed soft band (0.5–2 keV). Therefore, different levels of absorption affecting sources with the same count rate result in different measured fluxes. Usually, when performing statistical studies of populations of astrophysical objects detected in X-ray surveys, weights are assigned to each source to take into account the different sensitive areas of the survey at different fluxes. The curves of the area sensitive to a given flux (namely, the sky-coverage) for the 7 Ms CDF-S and 2 Ms CDF-N are shown in Fig. 4 (where we considered only the area covered by ≥1 Ms effective exposure, as per Section 2.1) for our two detection bands (see below). Since the source-detection algorithm used by L17 and X16 operates in the count rate space, sources with the same fluxes in the detection energy band can be detected over different areas, if their spectral shapes are different (i.e. the same count rate can correspond to different fluxes). This is usually not taken into account, and an average spectral shape (usually a power-law with Γ = 1.4) is assumed for all the sources, which can introduce errors up to a factor of several in the weights associated with each source (Vito et al. 2014).

3.2 Intrinsic column-density distribution

For each source, we report in Table 2 the logN_H that maximizes equation (10) (i.e. our best estimate of the intrinsic column density) and the 68 per cent confidence level uncertainties, corresponding to the narrowest logN_H interval in which the integrated P(logN_H) is 68 per cent. Fig. 5 presents the PDF(z) and P(logN_H) for a source as an example. The PDF of redshift, column density, flux (see Section 3.3) and luminosity (see Section 3.4) for every source in the sample are made available with the online version of this paper. Fig. 6 presents the best-fitting column density plotted against the nominal redshift of each source (i.e. the spectroscopic redshift or the redshift corresponding to the peak of PDF(z)).

Fig. 7 shows the binned N^tot(logN_H). Errors have been computed through a bootstrap procedure: we created a list of N_hz sources randomly chosen amongst the N_hz sources in our sample, allowing for repetition (i.e. one source can appear in the new list 0, 1 or multiple times). We then recomputed N^tot(logN_H). Repeating this procedure 1000 times, we derived 1000 distributions of N^tot(logN_H) or, equivalently, at each logN_H we derived 1000 values of probability. The 68 per cent confidence interval at each logN_H has been computed by sorting the corresponding 1000 probabilities and selecting the 16th and 84th quantile values. Liu et al. (2017) presented an X-ray spectral analysis of the 7 Ms CDF-S sources with >80 net counts. Fig. 7 is consistent with the column-density distribution presented in that work in an overlapping redshift bin (z = 2.4–5), in spite of the different assumed spectral models and performed analysis. Also, the priority we assigned to the different photometric-redshift catalogues is different than that used by Liu et al. (2017), resulting in some sources having different redshifts.

The flattening at logN_H ≲ 23 is due to the photoelectric cut-off detection limit. Indeed, the determination of the best-fitting N_H for an X-ray spectrum is largely driven by the detection of the photoelectric cut-off, which is a function of both redshift and column density. In particular, the photoelectric cut-off shifts to higher energies for increasing N_H and to lower energies for increasing redshift. Thus, at high redshift, unless a source is heavily obscured, the photoelectric cut-off is shifted outside the Chandra bandpass and, for small-to-moderate numbers of net counts, it is virtually impossible to determine the intrinsic column density. In Vito et al. (2013), we estimated the minimum column density that could be constrained for a typical (∼100 net counts) AGN in the 4 Ms CDF-S at z ∼ 3.5 to be logN_H ≈ 23.

We stress that the flattening at logN_H ≲ 23, due to the lack of information in X-ray spectra useful to constrain such low values of obscuration at high redshift, should not be considered real. We can derive the probability of a source having logN_H ≲ 23, but not the exact probability distribution in column density bins below this threshold, resulting in a flat behaviour. According to these considerations, we will adopt logN_H = 23 as the threshold between obscured and unobscured sources. This choice is further reinforced by the spectral simulations we performed in Appendix A.

3.3 Flux distribution and number counts

Fig. 9 shows the binned N^tot(F_{0.5–2 keV}) with 68 per cent confidence intervals derived with a bootstrapping procedure, similarly to what was done in Section 3.2 and Fig. 7. We note that, as we remained in the count rate space, the distributions extend slightly below the nominal flux-limit, since the assumed count rate limit corresponds to slightly different fluxes for different observed spectral shapes and response matrices (see Section 3.1).

3.3.1 AGN logN-logS at 3 ≤ z < 6

The cumulative version of equation (17), divided by the total area in deg², is the logN-logS relation of our sample. The 68 per cent confidence region derived through the bootstrap procedure is compared with results from previous works at z > 3 and z > 4 in Fig. 10. We derived the z > 4 logN-logS by repeating the procedures described in this section, using z = 4 as the lower limit for the considered redshift range. Bright and rare sources are not sampled properly by the pencil-beam, deep surveys such as those used in this work. Therefore, Fig. 10 displays the curves up to the flux at which the median expected number of sources from the bootstrap analysis is 1 in the area considered in this work (∼550 arcmin², see Fig. 4, corresponding to ∼6.5 deg⁻² sources).

Our results are in good agreement with previous measurements at logF_{0.5–2 keV} ≳ −16, which is the flux regime better sampled by wide surveys. At fainter fluxes, where we can exploit the excellent sensitivity of the Chandra deep fields, our results are consistent with previous results that used data from the 4 Ms CDF-S (Lehmer et al. 2012; Vito et al. 2014) and with the Ueda et al. (2014) curve. Lehmer et al. (2012) presented the number counts derived in the 4 Ms CDF-S down to slightly fainter fluxes than those reached in this work. This is due to the less-conservative approach they used to compute the sensitivity of the 4 Ms CDF-S compared to that used by L17 for the 7 Ms CDF-S, which we adopted here.

At z > 4, we can push the AGN number counts down to logF_{0.5−2 keV} ≈ −17. While our curve agrees with previous results at bright fluxes and with the Gilli, Comastri & Hasinger (2007) logN-logS reasonably well, it excludes very steep number counts such as those reported by Fiore et al. (2012) at faint fluxes. The selection in Fiore et al. (2012) made use of pre-determined optical positions of galaxies and adaptive X-ray detection bands, and it is therefore less conservative than a blind X-ray detection. However, about half of the 4 Ms CDF-S sources reported in that work at z > 4 were not detected even in the deeper 7 Ms exposures, leaving doubts about their detection significance (especially considering that the Fiore et al. (2012) analysis is plausibly affected by Eddington bias). This issue, together with the different, more recent photometric redshifts we used, which shift some of the Fiore et al. (2012) to lower redshifts, accounts completely for the discrepancies in Fig. 10.

By changing the limits of the integral of column density in equation (17), we derived separately the logN-logS for obscured (23 ≤ logN_H < 25) and unobscured (20 ≤ logN_H < 23) AGN (Fig. 11). The Gilli et al. (2007) X-ray background synthesis model underestimates the number of obscured AGN at high redshift and the steepening of our number counts of unobscured sources at logF_{0.5–2 keV} < −16.5 (see Fig. 11). Such steepening could be due to the population of star-forming galaxies, which begin providing a significant contribution to the number counts at the faintest fluxes probed by this work (e.g. Lehmer et al. 2012).

3.4 Luminosity distribution

P(L_X, logN_H, z) is normalized to unity. Equation (20) includes the correction for Eddington-bias applied in equation (5). However, we checked that neglecting this correction results in a similar luminosity distribution, with a slight decrease in the number of logL_X < 42.5 sources, balanced by an increase at higher luminosities, as expected. The luminosity probability distribution (P(L_X)) of each source is derived by integrating equation (20) over the considered redshift and column-density ranges. Fig. 5 presents P(L_X) for a source as an example. Table 2 reports for each source the luminosity value that maximizes P(L_X) (i.e. our best estimate of the luminosity of each source) and the 68 per cent confidence level uncertainties, corresponding to the narrowest interval containing 68 per cent of the total probability. If less than 68 per cent of P(CR) of a source lies above the count rate limit, we report the upper limit of P(L_X).

Fig. 12 shows the binned N^tot(L_X) with 68 per cent confidence intervals derived from a bootstrapping procedure, similar to what was done in Section 3.2 and Fig. 7.

4 OBSCURATION-DEPENDENT INCOMPLETENESS

The flux (or count rate) limit of a survey is effectively a limit in luminosity, redshift, and column density, as these are the physical parameters defining the flux (assuming our simple spectral model). Therefore, a low-luminosity AGN at a given redshift can be detected only if its column density is below a certain level, causing an obscuration-dependent incompleteness that affects the low-luminosity regime by favouring the detection of unobscured AGN (see also section 5.1.2 of Liu et al. 2017). Such incompleteness usually appears in X-ray studies of flux-limited AGN samples as a lack of detected sources in the low-luminosity and heavy-obscuration region of the L_X − N_H parameter space.

It is worth noting that this effect is not corrected by the Ω(CR) factor, which compensates for the different areas sensitive to different observed fluxes for detected sources, irrespective of their intrinsic properties (e.g. obscuration). In fact, while representing a potentially significant fraction of the population, low-luminosity and heavily obscured AGN are not detected at all, and are therefore not accounted for when deriving intrinsic properties (e.g. space density) of the low-luminosity population. This effect is particularly important when two obscuration-based subsamples are defined to compute the obscured AGN fraction, F_obsc, which would be biased towards lower values, and its trends with redshift and luminosity (e.g. Gilli et al. 2010). Therefore, before deriving F_obsc as a function of redshift and luminosity in Section 5, in this section we assess the effect of the obscuration-dependent incompleteness and derive suitable corrections.

The completeness level in the z–L_X plane is shown in Fig. 13 by computing the mean over the N_H axis separately for logN_H = 23–25 and logN_H = 20–23, weighted over an assumed intrinsic column-density distribution P(logN_H), which we here considered to be represented by Fig. 7. At logL_X ≲ 43, the two classes of objects are characterized by very different completeness levels: in the upper panel of Fig. 13, the detection of unobscured sources is complete up to z ∼ 6 for logL_X ≳ 42.5, while lower-luminosity AGN can be detected only at lower redshifts. Since the photoelectric cut-off for logN_H ≲ 23 shifts close to or below the lower-limit of the Chandra bandpass, for unobscured sources the completeness is almost independent of the column-density distribution. For this reason, the transition between 0 per cent and 100 per cent completeness is sharp, as it depends only on redshift and luminosity.⁸ The lower panel shows that for logN_H = 23–25 at a given redshift the transition is smoother and occurs at larger luminosities. In fact, contrary to the case of unobscured sources, the flux of an obscured source depends strongly on the particular N_H value: more heavily obscured sources can be detected only at higher luminosities than sources affected by milder obscuration. The darker stripe at logL_X ∼ 43–43.5 and z ∼ 3.5–4.5 is due to the inclusion of the hard band; heavily-obscured sources in those redshift and luminosity intervals are more easily detected in the hard band than in the soft band.

By projecting Fig. 13 over the luminosity axis (i.e. averaging over the redshift range), we derived the completeness curve as a function of luminosity (solid lines in Fig. 14). The averaging is weighted assuming an intrinsic redshift distribution characterized by a decline in the space density of high-redshift AGN proportional to (1 + z)⁻⁶ (e.g. Hiroi et al. 2012; Vito et al. 2014, solid lines). The unweighted projection corresponds to a flat distribution in redshift and is shown for completeness as dashed lines in Fig. 14. The specific redshift distribution has a small effect on the displayed curves: this behaviour is due to the negative curvature of absorbed X-ray spectra, which causes inverted K-corrections at increasing redshift.

To assess the dependence on the particular choice of column-density distribution, we also assumed the one from Gilli et al. (2007), in which, in particular, Compton-thick AGN are more numerous (see Fig. 7). While the completeness curve for unobscured sources is not sensitive to the particular choice of P(logN_H) (we do not show the curve corresponding to the Gilli et al. 2007 case for clarity), the completeness curve for obscured AGN is more severely affected by the Gilli et al. (2007) P(logN_H).

The bias affecting the obscured AGN fraction is linked to the relative completeness levels of the obscured and unobscured subsamples as a function of redshift and luminosity. To visualize this effect, the upper panel in Fig. 14 shows the ratio of the completeness characterizing obscured and unobscured subsamples as a function of luminosity. The ratio varies rapidly from logL_X = 42–43, rising from zero to ≈0.8. This is the regime in which the incompleteness effects are strongest: at low intrinsic luminosities, the presence or absence of even moderate obscuration is crucial for observing a flux above or below the survey sensitivity limit. At higher luminosities, only the most-obscured systems cannot be detected, and incompleteness is less severe.

Similarly, we derived the completeness as a function of redshift (Fig. 15) by projecting Fig. 13 over the redshift axis and assuming an intrinsic luminosity distribution (i.e. luminosity function). We used the pure density evolution model of Vito et al. (2014). The upper panel of Fig. 15 shows that the relative strength of incompleteness for obscured and unobscured sources does not significantly vary with redshift, but the evolution is slightly stronger for unobscured sources: this behaviour is due to the inverse K-correction characterizing obscured X-ray sources, which helps the detection of these systems at higher redshifts.

5 OBSCURED AGN FRACTION

As discussed in Section 3.2 and following Vito et al. (2013, 2014), throughout this work we used logN_H = 23 as the column-density threshold dividing obscured and unobscured AGN at z > 3. This value is higher than the commonly adopted threshold of logN_H = 22 in X-ray studies of the AGN population (e.g. Ueda et al. 2014; Aird et al. 2015; Buchner et al. 2015), which is also more in agreement with optical classification (e.g. Merloni et al. 2014). However, such low levels of obscuration are extremely difficult to detect at high redshift, where the photoelectric cut-off in X-ray spectra shifts below the energy band probed by X-ray observatories. The redshifted cut-off, together with the typical limited spectral quality of high-redshift sources, leads to a general overestimate of obscuration for low values of N_H, as can be seen in Appendix A. This effect is especially important when considering the obscured AGN fraction and its trends with other quantities, such as redshift and luminosity. We therefore prefer to define logN_H = 23, as the minimum column density of (heavily) obscured AGN at high redshift and discuss the dependence of (heavily) obscured AGN on redshift and luminosity. Unfortunately, this choice complicates the comparison with previous works utilizing more standard definitions. In Appendix C, we show that neglecting the full probability distributions of redshift and spectral parameters would lead to different results than those presented in the following subsections.

5.1 Obscured AGN fraction versus redshift

To derive an independent and parametric estimate of the dependency of F_obsc on redshift, we performed a Bayesian analysis⁹ of F_obsc(z), assuming a linear model of the form F^m_obsc = q_z + m_z(z − 3). We applied a flat prior to q_z and a minimally informative prior to m_z, such that |$P(m_{z})\propto (1+m_{z}^2)^{3/2}$| (e.g. VanderPlas 2014 and references therein). The latter is preferred over a flat prior on β, which would weight steep slopes more.¹⁰ We also associated a null probability to parameter pairs resulting in unphysical F_obsc values (i.e. negative or larger than unity) in the redshift interval probed by the data.

The resulting best-fitting values are |$q_{z}=0.50_{-0.05}^{+0.05}$| and |$m_{z}=0.00_{-0.06}^{+0.06}$|⁠. The red dashed lines in Fig. 16 enclose the 68 per cent confidence interval of 10⁶ MCMC realizations derived using this parametric approach, and the solid line is the median value. The parametric and non-parametric representations, derived independently, agree well, showing the robustness of our approach. The observed obscured AGN fraction is found in Fig. 16 to be flat at z = 3–5 at F_obsc ≈ 0.5. F_obsc is usually found to increase with redshift up to z ≈ 2–2.5, and then saturates at higher redshifts (e.g. Treister & Urry 2006; Hasinger 2008; Iwasawa et al. 2012; Buchner et al. 2015; Liu et al. 2017). We probe such behaviour up to z ∼ 6. Black circles in Fig. 16 are derived by correcting the red points using the solid lines in Fig. 15, causing F_obsc(z) generally to increase slightly, with no strong dependence on redshift, to values ≈0.65, consistent with the results by Buchner et al. (2015) at high luminosities (logL_X > 43.2) up to z ∼ 4. We could push the investigation of F_obsc(z) down to lower luminosities.

5.2 Obscured AGN fraction versus soft-band flux

Considering only the soft-band detected sources, we defined N_unobs(F_X) and N_obsc(F_X) from equation (17) by limiting the integral over column density to the ranges logN_H = 20–23 and 23–25, respectively, and summing only the soft-band-detected sources.

Similar to the process described in Section 5.1, we derived the parametric (assuming a linear model |$F_{\mathrm{obsc}}=q_{F_{\mathrm{X}}} + m_{F_{\mathrm{X}}}(\mathrm{log}F_{\mathrm{X}}-(-17))$| and non-parametric obscured AGN fraction as a function of flux. The results are shown in Fig. 17 and compared with the predictions of the Gilli et al. (2007) X-ray background synthesis model. The resulting best-fitting values are |$q_{F_{\mathrm{X}}}=0.43_{-0.04}^{+0.04}$| and |$m_{F_{\mathrm{X}}}=-0.04_{-0.06}^{+0.06}$|⁠.

5.3 Obscured AGN fraction versus luminosity

We compare the high fraction of heavily obscured AGN (F_obsc ≈ 0.7–0.8) we derived at high luminosities (logL_X > 43) with previous findings by Aird et al. (2015) and Buchner et al. (2015), where obscured AGN are defined in a more ‘standard’ way as those characterized by logN_H = 22–24. The column-density threshold used to separate obscured and unobscured AGN in their works, logN_H = 22, while useful for comparison purposes with results at lower redshift, is not suitable at high redshift, where the photoelectric cut-off shifts at low energies, close to or even below the lower energy boundary of the Chandra bandpass and is therefore poorly constrained, as discussed in Section 7 and Appendix A. For completeness, we also show the results of Georgakakis et al. (2015), where the obscured AGN fraction is derived by comparing X-ray and UV luminosity functions, and therefore is not directly comparable to the other purely X-ray defined curves.

While at logL_X > 43 incompleteness effects are negligible, they dominate the observed trend of F_obsc at lower luminosities in Fig. 18. Applying the corrections discussed in Section 4 results in values, represented by the black circles in Fig. 18, consistent with those at logL_X > 43 (i.e. F_obsc ≈ 0.7–0.8), although at logL_X < 43 a slightly decreasing trend of F_obsc, down to ≈0.6, with decreasing luminosity remains visible. However, we caution that the applied corrections depend quite strongly on the particular intrinsic column density and luminosity distributions. Also, the possible detection of star-forming galaxies, which have typically steep spectra, at the lowest luminosities probed by this work (logL_X ≈ 42–42.5) may decrease F_obsc in such regime. We therefore do not consider such trend to be significant.

The curve of Buchner et al. (2015) appears to suffer less from similar issues at low luminosities, although they reported a decrease of F_obsc(L_X) at low luminosities. The difference with our observed result is probably due to the different procedure used; Buchner et al. (2015) applied a Bayesian procedure that disfavours strong variations of F_obsc over close redshift and luminosity bins. At high redshift, the faint regime is not well sampled, probably causing the obscured AGN fraction at low luminosities to be dominated by the priors, i.e. to be similar to the value at lower redshift and/or higher luminosity, where incompleteness issues are less severe, alleviating this issue. The slight decrease of F_obsc towards low luminosity, which they ascribed to different possible physical mechanisms, may be at least partly due to incompleteness effects.

A strong anti-correlation between F_obsc and luminosity is usually found at low redshifts. For instance, grey empty symbols in Fig. 18 represent the fraction of AGN with logN_H > 23 derived by Burlon et al. 2011 at z < 0.3. This behaviour appears not to hold at high redshift, or at least to be much less evident. Comparing our points with the Burlon et al. (2011), we note the positive evolution of F_obsc from the local Universe to high redshift, which is stronger at high luminosities. In Vito et al. (2014), where we derived similar results from a combination of different X-ray surveys, we ascribed the larger fraction of luminous obscured AGN at high redshift than at low redshift to the larger gas fractions of galaxies at earlier cosmic epochs (e.g. Carilli & Walter 2013), which can cause larger covering factors and/or longer obscuration phases. This description is especially true if luminous AGN are preferentially triggered by wet-merger episodes (e.g. Di Matteo, Springel & Hernquist 2005; Menci et al. 2008; Treister et al. 2012), whose rate is expected to be higher in the early Universe. In this case, chaotic accretion of gas on to the SMBH can produce large covering factors.

Figs 16 and 18 do not account for the covariance between the F_obsc − z and F_obsc − L_X trends, due to our data set being flux limited. This could in principle bias the results when investigating F_obsc separately as a function of one of the two parameters. To check for this possible effect, we derived the trend of F_obsc with redshift separately for two luminosity bins, logL_X = 43–43.5 and 43.5–45, chosen not to suffer significantly from incompleteness and to include approximately the same number of sources. For both subsamples, F_obsc is consistent with being flat at similar values, reassuring us that this result does not strongly depend on the intrinsic trend with redshift. We also repeated the procedure dividing the sample in two redshift bins and deriving the trend of F_obsc with luminosity in each of them. The results are consistent with those reported above in this section within the (larger) uncertainties.

6 EVOLUTION OF THE HIGH-REDSHIFT AGN POPULATION

6.1 The AGN XLF at high redshift

The AGN XLF can be derived from equation (22) by integrating over narrower redshift intervals and dividing the results by the volume sampled by the surveys in each redshift bin. Fig. 19 presents the XLFs in two redshift bins, chosen to include approximately the same numbers of sources. As stated in Section 3.4, ignoring the correction for Eddington bias does not change significantly the luminosity distribution of our sample, and, as a consequence, the derived XLFs. In particular, the slope of the faint end does not significantly vary when applying or neglecting such a correction. Our XLFs are fairly consistent with previous observational results from Aird et al. (2010), Ueda et al. (2014) and Vito et al. (2014, 2016). They are also in agreement with Buchner et al. (2015) at z < 3.6, which do not provide an analytical expression for the XLF, but rather confidence regions in redshift intervals. At higher redshifts, the Buchner et al. (2015) region exceeds our results. Similarly, the simulated XLFs by Habouzit, Volonteri & Dubois (2017) appear to overpredict the density of low-luminosity AGN.

Giallongo et al. (2015) applied a detection procedure based on searching for clustering of photons in energy, space and time, at the pre-determined optical positions of CANDELS-detected galaxies (see also Fiore et al. 2012), which allowed the authors to push the Chandra sensitivity beyond the limit reachable by blind detection methods (i.e. with no previous knowledge of the positions of galaxies). The large number of low-luminosity AGN at high redshift detected by Giallongo et al. (2015) suggests that AGN could have an important role in cosmic reionization (e.g. Madau & Haardt 2015; but see also Ricci et al. 2017).

Giallongo et al. (2015) derived the UV LF of their sample of X-ray detected, high-redshift AGN candidates by deriving the absolute UV magnitude from the apparent optical magnitude in the filter closest to rest-frame |$1450\, {{^{\circ}_{\rm A}} }$| at the redshift of each source. We transformed their UV LF into an X-ray LF (see Fig. 19) by assuming a SED shape F_ν ∝ ν^−0.5 between 1450 and |$2500\, {{^{\circ}_{\rm A}} }$| (as in Georgakakis et al. 2015), the Lusso et al. (2010) α_ox and Γ = 1.8 for the X-ray spectrum (see also Vito et al. 2016). Our results do not support the very steep LFs derived by Giallongo et al. (2015). Parsa et al. (2017) recently disputed the very detection of the faint z > 4 AGN in Giallongo et al. (2015). Moreover, Cappelluti et al. (2016), applying a similar procedure which makes use of pre-determined optical positions, did not find such a large number of z > 4 sources. Other issues plausibly affect the assessment of the population of faint, high-redshift AGN, such as the uncertainties in the photometric redshifts, the expected presence of Eddington bias (which leads to overestimating the number of detected faint sources), and the XRBs contribution to X-ray emission at faint fluxes (e.g. Lehmer et al. 2012). Moreover, the uncertainties related to the assumed UV/X-ray spectral slope (discussed by Giallongo et al. 2015) may affect the conversion between UV magnitude and X-ray luminosity, although reasonable values for that conversion cannot produce the tension between their slope of the AGN XLF faint end and our results.

At the low luminosities probed by the 7 Ms CDF-S (logL_X = 42 − 43), X-ray emission from XRB may provide a non-negligible contribution (e.g. Bauer et al. 2002; Ranalli, Comastri & Setti 2003, 2005; Lehmer et al. 2010; Mineo, Gilfanov & Sunyaev 2012). The green dashed line in Fig. 19 (right-hand panel) is the XLF of star-forming galaxies at z = 4, where the emission is mostly due to high-mass XRB, from the model by Volonteri et al. (2017, see in particular their section 2.2 for assumptions and caveats). Briefly, they mapped the stellar-mass functions of Song et al. (2016) into an XLF through the main sequence of star-forming galaxies (Salmon et al. 2015) and the SFR-L_X scaling relation by Fragos et al. (2013, model 269), which Lehmer et al. (2016) demonstrated to be the best representation at z = 0–2.5. At logL_X ≳ 42.5, the luminosity limit we applied to our luminosity functions, as the 42 ≲ logL_X ≲ 42.5 regime, was strongly affected by incompleteness; the density of star-forming galaxies is a factor of ≳10 lower than the AGN density. We conclude that at such luminosities, our sample is not significantly contaminated by star-forming galaxies. At lower luminosities, the XLF of star-forming galaxies is comparable to, or even exceeds, the AGN XLF, in agreement with Vito et al. (2016), who found that the X-ray emission in galaxies individually undetected in the 7 Ms CDF-S is mostly due to XRB. Similar conclusions are reached by scaling the star-formation rate functions (SFRF; i.e. the comoving space densities of galaxies per unit star-formation rate) from Smit et al. (2012) and Gruppioni et al. (2015), at z ≈ 3.5 and 4, respectively, into XLF using the best-fitting, redshift-dependent Fragos et al. (2013) and Lehmer et al. (2016) relations.

6.2 The evolution of the AGN space density and comparison with the galaxy population

Integrating the XLF over luminosity, we derived the comoving space density of AGN in three luminosity bins (Fig. 20). The decline of the space density in Fig. 20 appears to be steeper for low-luminosity AGN than for luminous sources. However, incompleteness significantly affects the detection of low-luminosity AGN, as can be seen in Fig. 13. In fact, applying the incompleteness correction discussed in Section 4 to the density of low-luminosity AGN (empty circles in Fig. 20), a slightly flatter behaviour is derived. To quantify the decline of the space density, we fitted a model Φ(z) = A((1 + z)/(1 + z₀))^p, where z₀ = 3, to the points through a simple χ² minimization procedure, before (dashed lines in Fig. 20) and after the correction for incompleteness was applied. Fig. 20 displays the best-fitting model corrected for incompleteness only for the low-luminosity bin (dotted line), as at higher luminosities incompleteness is negligible. Table 3 summarizes the best-fitting values for the total sample and the three luminosity-based subsamples.

The evolution of the space density of luminous (logL_X ≳ 44) AGN in Fig. 20 is consistent with previous results by, e.g., Hiroi et al. (2012), Kalfountzou et al. (2014), Vito et al. (2014) and Marchesi et al. (2016) (i.e. Φ ∝ (1 + z)^p with p ∼ −6, see Table 3). Luminous (logL_X > 44) AGN are usually hosted by massive (M_* ≳ 10¹⁰M_⊙) galaxies (e.g. Xue et al. 2010; Yang et al. 2017). Fig. 20 shows that the space densities of luminous AGN and massive galaxies (from Davidzon et al. 2017) evolve with similar slopes, suggesting that the decline of the space density of logL_X > 44 AGN is merely driven by the evolution of the galaxy number density. This interpretation is further supported by the recent findings of Yang et al. (submitted), who derive a positive evolution of the average black hole accretion rate in massive galaxies from low redshift up to z ≈ 3, followed by a steady behaviour up to z ≈ 4. This trend would be in contrast with the well-established decline of the space density of luminous AGN if the number of massive galaxies itself did not decrease.

Thanks to the unmatched sensitivity of the Chandra deep fields, in particular the 7 Ms CDF-S, coupled with the deep multiwavelength coverage, which allows a nearly complete multiwavelength identification of the X-ray sources, we could push the investigation of the evolution of the AGN space density down to logL_X = 42.5 and up to z = 6. We found a slightly steeper decline for low-luminosity sources than for high-luminosity AGN, even after applying our fiducial corrections for incompleteness (Fig. 20 and Table 3). However, such corrections depend on the assumed intrinsic distribution in redshift, column density and luminosity (see Section 4). Moreover, the fit we performed is an approximation (although widely used in literature), as the error bars of the points in Fig. 20, are derived through a bootstrap procedure, which formally returns a measurement of the dispersion of the bootstraped distribution. For these reasons, we do not consider the best-fitting values reported in Table 3 as definitive evidence, but rather a suggestion, for a differential evolution of the AGN space density in different luminosity regimes. However, the strong flattening of the XLF faint-end found by Georgakakis et al. (2015) is ruled out by Fig. 19. Future X-ray missions, discussed in Section 7.2, will shed light on the evolution of the population of low-luminosity AGN at higher redshifts than those probed in this work.

Intriguingly, the evolution of the space density of low-luminosity AGN is also steeper than the density evolution of both low-mass and high-mass galaxies (see Fig. 20), as derived by Davidzon et al. (2017). This discrepancy requires an evolution of one or more of the physical parameters driving black hole accretion for low-luminosity AGN, other than the mere evolution of the number of galaxies. For instance, a steeper evolution of the AGN density than galaxy density can be expected if the duty cycle and/or occupation fraction decrease at increasing redshift, especially for low-mass galaxies.

6.3 XLF and space density of obscured and unobscured AGN

Fig. 21 (left-hand panel) shows the luminosity functions separately for obscured and unobscured sources, obtained by performing the integral over column density in equation (22) over the ranges logN_H = 23–25 and 20–23, respectively, compared with the XLFs of Hasinger, Miyaji & Schmidt (2005, used by Gilli et al. 2007 for their X-ray background synthesis model) in the same N_H bins. We also report the results from Aird et al. (2015), where obscured and unobscured AGN are defined as those obscured by column densities logN_H = 22–24 and 20–22, respectively. Points in the first luminosity bin (empty circles) are heavily affected by the obscuration-dependent incompleteness discussed in Section 4, as is evident by the drop of the density of obscured AGN with respect to the unobscured subsample. Fig. 21 (right-hand panel) displays the evolution of the space density of obscured and unobscured AGN. Table 3 summarizes the best-fitting values for the obscured and unobscured subsamples. We found no significant difference in the evolution of the space density of obscured and unobscured AGN. This conclusion is consistent with Fig. 16, where the obscured AGN fraction is flat with redshift. After correcting for incompleteness, a slightly shallower decline of the space density of unobscured AGN is derived, because of the faster decrease of detection completeness with increasing redshift for unobscured sources than for obscured sources (see Fig. 15). However, the uncertainties due to the procedure of correcting for incompleteness and the relatively small size of the sample prevent us from drawing strong conclusions.

6.4 Black hole accretion rate density

Fig. 22 displays the BHAD corresponding to the XLFs of Ranalli et al. (2016), Aird et al. (2015), Georgakakis et al. (2015), Ueda et al. (2014) and Vito et al. (2014), converted into BHAD as described above,¹² and that derived by Delvecchio et al. (2014) exploiting Herschel data. We also add the BHAD predicted by simulations (Lodato & Natarajan 2006; Volonteri 2010; Shankar, Weinberg & Miralda-Escudé 2013; Bonoli, Mayer & Callegari 2014; Sijacki et al. 2015; Volonteri et al. 2016a) under different assumptions of BH seed mass and growth mechanisms. Except for Sijacki et al. (2015), who directly reported the BHAD, we derived the curves by differentiating the BH mass densities presented in those works.

Simulations tend to derive BHAD exceeding by about one order of magnitude the observed values at z > 3 (see Fig. 22). This discrepancy is usually ascribed to the contribution of continuous, low-rate accretion, difficult to detect even in deep surveys. In Vito et al. (2016), we used a stacking procedure to demonstrate that such contribution is negligible compared to the BHAD due to X-ray-detected AGN (see also e.g. Pezzulli et al. 2017). With the refined sample of high-redshift AGN of this work, we confirmed this result up to z ∼ 6 and suggest that simulations tend to overestimate the BHAD at high redshift.

A similar issue is known to affect the faint end of the AGN LF at high redshift, with simulations finding steeper slopes than the observed ones. Volonteri et al. (2016a) suggested that this discrepancy could be due to supernovae feedback (i.e. the removal of gas in the galaxy responsible of both nuclear accretion and star formation as a consequence of supernovae activity) being underestimated in low-mass haloes in simulations. Considering only massive (M > 5 × 10¹¹M_⊙) dark-matter haloes, where supernovae feedback has little effect due to the deeper galactic gravitational potential well, improves the agreement with observations. This statement is also true for the BHAD, as can be seen in Fig. 22, where the results of Volonteri et al. (2016a) are presented without any mass cut (continuous blue curve) and for massive haloes only (cyan curve).

The common evolution of star formation rate density (SFRD) and the BHAD at z ≲ 2.5 is generally considered to be a manifestation of the BH–galaxy co-evolution. However, these two quantities are found to evolve with different slopes at higher redshifts (e.g. Silverman et al. 2008; Aird et al. 2015, see e.g. Haiman, Ciotti & Ostriker 2004 for a theoretical interpretation), mainly through extrapolations of the AGN XLF. In this work, we can probe this behaviour, whose physical causes are difficult to identify, down to the faint flux limits of the Chandra deep fields. Most of the BHAD we derived at high redshift is due to luminous (logL_X > 44) AGN, which contribute 70-90 per cent to the black circles in Fig. 22; therefore, the uncertainties on the density of low-luminosity AGN have no strong impacts in this respect. In particular, the contribution of accretion in individually undetected galaxies is negligible (empty black symbols from Vito et al. 2016 in Fig. 22). Indeed, the evolution the BHAD is similar to that of the space density of luminous AGN, both declining by a factor of ≳10 from z = 3 to z = 6. In contrast, the evolution of the SFRD is more similar to that of low-mass galaxies than massive systems, both declining by a factor of ∼3 in the same redshift range (cf. Figs 22 and 20). This discrepancy may mark a breakdown of the BH–galaxy co-evolution at high redshift.

7 FUTURE PROSPECTS

7.1 Spectroscopic follow-up of high-redshift AGN candidates

One of the main sources of uncertainty affecting our knowledge of the demography of high-redshift AGN is the persistently low spectroscopic identification rates of the candidates. The deep Chandra surveys allow sampling of the bulk of the population at high redshift, which is constituted of low luminosity and/or obscured sources. However, these objects are typically faint in the optical/NIR bands (e.g. ≈ 80 per cent of our sources located in the CANDELS/GOODS-S area have H > 24), and are therefore difficult to follow up spectroscopically with current facilities in order to confirm their high-redshift nature. While photometric redshifts of good quality are available in the fields we considered, and despite the careful treatment of their PDF, they may introduce non-negligible uncertainties in the final results (e.g. see section 4.3 of Georgakakis et al. 2015). This issue is less problematic for shallower X-ray surveys, which sample AGN with brighter optical/NIR counterparts, but miss the sub-population of low-luminosity AGN.

Current facilities that could be used effectively to increase the spectroscopic completeness of our sample include ALMA¹³ (e.g. Walter et al. 2016), the Keck observatory¹⁴ and VLT/MUSE.¹⁵ The last, in particular, with 27 h of integration time in the Hubble Ultra Deep Field, has already been able to obtain spectroscopic redshifts for sources below the Hubble detection limit (Bacon et al. 2015). In the near future, facilities like James Webb Space Telescope (JWST),¹⁶ and, in the somewhat longer term, the Extremely Large Telescopes (with ≳30 m diameter mirrors) will provide the required spectroscopic sensitivity to perform ultradeep spectroscopic follow-up campaigns aimed at identifying high-redshift AGN candidates with extremely high efficiency (see e.g. the discussion in Brandt & Alexander 2015).

7.2 The case for Athena and Lynx

The understanding of the SMBH formation and growth mechanisms in the early Universe requires the collection of much larger samples of X-ray-selected AGN and galaxies at higher redshifts and lower luminosities than those probed by current X-ray surveys. Moreover, while Chandra and XMM–Newton, jointly with large-area optical/IR surveys, are successfully improving our knowledge of the bright quasar population up to z ≈ 7 (e.g. Brandt & Vito 2017), the bulk of the AGN population, constituted by low luminosity and possibly obscured systems, are currently completely missed beyond z ≈ 5–6 even by the deepest X-ray surveys available, and Chandra observations substantially deeper than the 7 Ms CDF-S are not foreseeable. The Athena X-ray Observatory (Nandra et al. 2013)¹⁷, which will be launched in ∼2028, is expected to detect hundreds of ≲L_* AGN at z > 6 (Aird et al. 2013), thanks to its survey power a factor of ∼100 faster than Chandra or XMM–Newton, boosting our knowledge of the high-rate accretion phases, which likely dominate SMBH growth at high redshift (Vito et al. 2016).

Simulations suggest that different combinations of the physical parameters governing SMBH formation and early growth (e.g. seed-mass and Eddington ratio distributions, occupation fraction, feedback mode, and their dependence on the host-galaxy and environment properties) produce different shapes of the faint-end (logL_X ≲ 42) of the AGN LF. The green-striped region in Fig. 19 (right-hand panel) is the locus of the XLF models of Volonteri et al. (2017) at z = 4, under different assumptions. In particular, the upper bound is derived using a minimum initial BH mass of 10² M_⊙ and occupation fraction equal to unity, while for the lower bound a minimum initial mass of 10⁵M_⊙ and a modified occupation fraction from Habouzit et al. (2017), flagged as ‘H17low1’ in Volonteri et al. (2017). Different combinations of parameters are mostly encompassed by these boundaries. The models are normalized to the observed points, which explain the small spread at logL_X ≳ 43. However, at lower luminosities, which even the currently deepest surveys cannot sample, different combinations of parameters produce significantly different XLF shapes. The spread of the XLF shapes in this faint regime is expected to increase at higher redshifts (see fig. 12 of Volonteri et al. 2017). A tremendous step forward towards understanding the early growth of SMBH will be provided by Lynx (also known as X-ray Surveyor; Weisskopf et al. 2015), which, with its Chandra-like angular resolution and ≈25 times higher sensitivity, will reach an expected flux limit of logF ≈ −18.5 in a blind 4 Ms survey, sufficient to probe the faint end of the XLF down to logL_X ≈ 40–42 at z ≲ 10, providing unprecedented constraints on the combinations of parameters driving SMBH formation and early growth.

The constraints derived by X-ray telescopes will be complemented by deep optical/IR observations provided by future facilities such as JWST and WFIRST¹⁸, necessary to identify the X-ray sources. Recently, several works have focused on the multiwavelength properties of emission from BHs during the early accretion phases. For instance, Volonteri et al. (2017) proposed a colour-colour selection based on JWST bands to separate AGN from normal galaxies at high redshift. In particular, the characteristic shape of the SED of accreting direct-collapse BH (with initial M_BH = 10⁴–10⁵ M_⊙) has been used by Natarajan et al. (2017) to propose a JWST colour-colour method to select such objects. In contrast, optical/IR emission from accreting light-seeds BH (M_BH = 10²–10³ M_⊙) probably will be outshined by stellar emission in the young host galaxies. The joint effort of the next generation X-ray observatories and optical/IR telescopes therefore will open a completely new window on the cosmic epoch where SMBH formed.

8 CONCLUSIONS

We have investigated the population of 3 ≤ z < 6 AGN detected in the two deepest X-ray fields, the 7 Ms CDF-S and the 2 Ms CDF-N, selected on the basis of a careful assessment of the photometric-redshift uncertainties for sources lacking spectroscopic identification. For each source, uncertainties on redshift, column density, flux and luminosity were accounted for by using the full probability distributions of such parameters, and a statistical correction for Eddington bias was applied. Our main results are the following:

• We derived the best-constrained column-density distribution of moderate-luminosity obscured AGN at high redshift, finding it to peak at logN_H ∼ 23.5–24 (see Section 3.2).

• Comparing the mean photon index of X-ray spectra of high-redshift, L_X < L_* AGN with results derived for luminous quasars, we found no significant variation with redshift and luminosity (see Appendix Section B).

• The obscured AGN fraction (F_obsc), where the column-density threshold to define obscured versus unobscured AGN is set to logN_H = 23, is consistent with being flat with redshift at observed values of ≈0.5–0.6. Once corrected for incompleteness, we derive F_obsc ≈ 0.6–0.7. This result agrees with previous works finding a saturation of this quantity at z ≳ 2, after a positive evolution at z ≈ 0–2 (see Section 5.1).

• Thanks to the sensitivity of the Chandra deep fields, we could probe the number counts of high-redshift AGN down to fainter fluxes, especially at z > 4, than those reached by wider and shallower surveys. Comparing the logN-logS of our sample with the predictions from the Gilli et al. (2007) X-ray background synthesis model, which was constrained at lower redshifts, separately for obscured and unobscured AGN, we found an excess of obscured sources, likely caused by the positive evolution of F_obsc from the local Universe to high redshift (see Section 3.3.1).

• We found F_obsc ≈ 0.7–0.8 at logL_X ≳ 43, somewhat larger, but still in agreement with other works performed at similar redshifts. In contrast with previous results, derived especially at low redshift, there is no significant anticorrelation between F_obsc and luminosity at logL_X ≳ 43, suggesting a stronger positive evolution of F_obsc towards earlier cosmic times for luminous AGN than for moderate-luminosity sources (see Section 5.3).

• At low luminosities (logL_X ≲ 43), the observed behaviour of F_obsc is strongly affected by sample incompleteness: low luminosity, obscured sources are more difficult to detect than unobscured sources, biasing the observed result towards low values of F_obsc. We suggest that recent claims of a turn-down of F_obsc towards low luminosities at high redshift may be at least partly due to similar incompleteness issues. By assuming intrinsic distributions of sources in N_H and L_X space, we applied tentative corrections that led F_obsc to be fairly consistent with the values at higher luminosities, although a possible slight decrease down to values of ≈0.6 is still visible (see Sections 5.3 and 4).

• The space density of luminous (logL_X ≳ 44) AGN declines from z = 3 to z = 6, as (1 + z)^d, with d ≈ −6, in agreement with previous results. The XLF in the redshift range probed by this work does not show significant steepening of the faint end, which has been recently linked with a possible significant contribution of low-luminosity AGN to cosmic reionization. Our results tend to exclude this hypothesis. By contrast, the space density of AGN with low-to-moderate luminosity appears to decline slightly faster than for high-luminosity AGN. However, uncertainties due to incompleteness might affect this result. The evolution of the AGN XLF faint end at high redshift is particularly important for placing constraints on the combination of the physical parameters driving the formation and early growth of SMBHs, one of the goals of future X-ray observatories (see Sections 6 and 7.2).

• The space densities of luminous AGN and massive galaxies show a similar trend with redshift, suggesting that the evolution of luminous AGN is driven by the underlying galaxy population. By contrast, the space density of low-luminosity AGN declines with increasing redshift faster than the space density of both high-mass and low-mass galaxies. This difference suggests an evolution in the physical parameters regulating nuclear activity, such as duty cycle and occupation fraction. We also derived the BHAD at z = 3–6, and noted that, while the BHAD and SFRD track each other at z ≈ 0–2.5, at higher redshifts the BHAD decreases with a faster rate than the SFRD (see Section 6).

Acknowledgements

We are grateful to the referee, A. Georgakakis, for useful comments that helped in clarifying several aspects of this work. FV, WNB and GY acknowledge support from Chandra X-ray Center grant GO4-15130A, the Penn State ACIS Instrument Team Contract SV4-74018 (issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060), and the V.M. Willaman Endowment. BL acknowledge support from the National Natural Science Foundation of China grant 11673010 and the Ministry of Science and Technology of China grant 2016YFA0400702. YQX acknowledges support from the 973 Program (2015CB857004), NSFC-11473026, NSFC-11421303, the CAS Frontier Science Key Research Program (QYZDJ-SSW-SLH006), and the Fundamental Research Funds for the Central Universities. We warmly thank Harry Ferguson and Dritan Kodra for providing PDFs of CANDELS photometric-redshifts, and Iary Davidzon, Carlotta Gruppioni, and Melanie Habouzit for providing useful data. The Guaranteed Time Observations for the CDF-N included here were selected by the ACIS Instrument Principal Investigator, Gordon P. Garmire, currently of the Huntingdon Institute for X-ray Astronomy, LLC, which is under contract to the Smithsonian Astrophysical Observatory; Contract SV2-82024.

REFERENCES

SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

online_figures

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

APPENDIX A: CHECKING THE ACCURACY OF THE X-RAY SPECTRAL ANALYSIS

In Section 3, we assumed a simple spectral model, an absorbed power law with photon index fixed to Γ = 1.8, as a representation of the X-ray spectra of the sources in our high-redshift sample, and performed a spectral analysis to derive the probability distribution of the intrinsic column density. Although the use of more complex spectral models is generally precluded by the limited number of counts in the X-ray spectra of the typically faint AGN in our sample, we check in this appendix for potential biases arising from this assumption using X-ray spectral simulations.

We assumed that the intrinsic X-ray spectra in cases of column densities logN_H > 22 are well represented by the MYTorus spectral model¹⁹ (Murphy & Yaqoob 2009), which fully accounts for the transmission and scattered components, as well as the production of the iron K edge and K_α and K_β emission lines. We fixed the Galactic absorption to logN_H = 20 (the characteristic value towards the CDF-S; Kalberla et al. 2005), Γ = 1.8, the inclination angle Θ = 90° (i.e. edge-on configuration), and the normalization of the scattered and line components to unity with respect to the transmitted component. The MYTorus model does not allow lower values of column densities, and therefore for logN_H < 22 we used the xspecplcabs model, which also takes into account Compton scattering. We assumed three values of observed soft-band flux, logF = −17, −16 and −15, corresponding to sources close to the flux limit of the 7 Ms CDF-S, faint sources and bright sources (in deep X-ray surveys), respectively. We also considered three different redshifts (z = 3, 4 and 5) and eight different values of column density (from logN_H = 21 to 24.5 in steps of 0.5 dex). For each of the resulting 72 combinations of these parameters, we simulated 10 X-ray source and background spectra with a 6 Ms effective exposure (typical of the CDF-S), allowing for Poissonian fluctuations, using the fakeit command in xspec. The background spectrum is sampled from the real CDF-S background. We then fitted the simulated spectra with our simple model and derived the best-fitting column density and its errors at the 90 per cent confidence level.

Fig. A1 shows the best-fitting column densities against the simulated ones for the different redshift-flux combinations. The output values in cases of large simulated column densities are consistent with the input ones. When decreasing the simulated flux and increasing the redshift, low values of column density (logN_H ≲ 22.5 and, less severely, logN_H ∼ 23) become increasingly difficult to constrain. In particular, the output column densities of sources simulated with logN_H ≤ 22 are never constrained to be lower than that value, which is the widely used threshold separating obscured and unobscured AGN at low redshift. This result is due to the photoelectric cut-off shifting close to, and even below, the low-energy spectral limit of the Chandra bandpass. Small differences in the intrinsic photon index (within the expected range of Γ = 1.7–2.0) do not affect significantly the output column densities. We therefore conclude that the spectral model we assumed, despite its simplicity, is a useful and appropriate approximation, at least for the redshift and flux ranges probed by this work.

APPENDIX B: The X-ray power-law photon index of high-redshift, sub-L_* AGN

The photon index (Γ) of the primary power-law emission in AGN X-ray spectra can be used to estimate the Eddington ratio (λ_Edd = L/L_Edd) of a source through the Γ–λ_Edd relation (e.g. Shemmer et al. 2006, 2008; Risaliti, Young & Elvis 2009; Brightman et al. 2013, 2016; Fanali et al. 2013, but see the caveats discussed in Trakhtenbrot et al. 2017), according to which higher λ_Edd would produce steeper Γ. This relation is likely produced by the coupling between the accretion disc and the hot corona, and can therefore be used to investigate whether the accretion mode depends on the Eddington ratio at different redshifts.

The spectral analysis presented by Liu et al. (2017) on bright (>80 net counts) X-ray-selected AGN in the 7 Ms CDF-S reveals a median value of Γ = 1.82 ± 0.15, typical of Seyfert galaxies. Evidence for steep photon indexes in a few luminous Type-I QSOs at z > 6 have been found by Farrah et al. (2004), Page et al. (2014) (whose results are disputed by Moretti et al. 2014) and Ai et al. (2016). However, studies of samples of z = 4–7 Type-I QSOs did not find a significant evolution of the average Γ in luminous QSOs with redshift (e.g. Shemmer et al. 2005; Vignali et al. 2005; Nanni et al. 2017). Flatter photon indexes (Γ ≈ 1.1–2.0) have been derived by Vignali et al. (2002) for three spectroscopically confirmed z > 4, moderate-luminosity AGN in the CDF-N, probably due to the presence of intrinsic absorption.

Here, we compare the photon indexes of L ≲ L_* AGN with results derived for much more luminous objects. We selected a subsample of bona fide unobscured sources at z = 3–6 from Table 2, requiring that the best-fitting column density is logN_H < 23 at the 68 per cent confidence level and a signal-to-noise ratio (SNR) > 3 in the 0.5–7 keV band. We performed a more detailed spectral analysis than that reported in Section 3, by allowing both the intrinsic column density and the photon index to vary. Fig. B1 (blue squares) shows the resulting best-fitting values of Γ as a function of intrinsic luminosity for sources with >110 net counts in the 0.5–7 keV band. This threshold was chosen in order to obtain reasonable uncertainties on the best-fitting Γ for individual AGN.

For the fainter sources, we derived the average photon indexes by performing a joint spectral analysis in luminosity bins chosen to include approximately the same number of sources. During the fitting, the column density was left free to vary for each spectrum, while the photon index parameters were linked amongst all sources, in order to derive an average value. Red circles in Fig. B1 present the results of this joint analysis for sources with <110 net counts. Finally, we performed the joint spectral analysis in three luminosity bins, including both bright and faint objects (black points). Fit results are reported in Table B1 for the individual fit of bright sources and for the joint fit of faint and all (bright and faint) sources.

Fig. B1 presents as continuous and dashed grey lines the median Γ and the standard deviation of the Γ distribution (Γ = 1.82 ± 0.15), respectively, derived by Liu et al. (2017) for sources with >80 net counts in the 7 Ms CDF-S. Our sources are on average broadly consistent with such values, justifying the use of Γ = 1.8 for all sources in Section 3. The photon indexes of the luminous sources in our sample appear to be on average slightly flatter. This result is probably driven by the small number of luminous AGN considered here, which is due to the limited area covered by the deep surveys we used. Moreover, the best-fitting average Γ in the highest luminosity bin is quite strongly driven by the 1–2 brightest sources, which have Γ ≈ 1.6 and alone provide the 38–59 per cent of the net counts for the joint fitting in this luminosity bin. Neither of these two sources are detected in radio bands, and, in general, discarding the minority of sources (2 out of 23) detected in radio surveys does not change the results.

Our findings at 42 < logL_X ≲ 44 are compared in Fig. B1 with the average Γ of optically selected, luminous (logL_X > 45) Type-I QSOs at low-to-intermediate redshift (z = 1–5, Just et al. 2007) and z = 4–7 (Vignali et al. 2005; Shemmer et al. 2006; Nanni et al. 2017). Our results are consistent with both the median value of Γ derived by Liu et al. (2017) for similar luminosities at all redshifts and the average values of luminous QSOs up to z ≈ 7.

APPENDIX C: THE EFFECT OF NEGLECTING THE FULL PROBABILITY DISTRIBUTION OF SPECTRAL PARAMETERS

In Section 3, we derived the full probability distribution of flux, and luminosity by convolving the probability distribution of count rate and redshift, and using the proper count-rate-to-flux conversion for each source. Such information has been used in Section 5 to derive the dependence of the obscured AGN fraction on redshift, flux and luminosity.

To check the importance of applying this careful but computationally expensive procedure, we repeat the analysis in Section 5 by considering only the nominal value of parameters for each source, according to the following steps. (1) A single redshift (z^peak), corresponding to the spectroscopic redshift or the peak of PDF(z), was assigned to each source. We did not include objects with no redshift information. (2) We adopted the values corresponding to the peak of the count rate (CR^peak) and column-density (N_H^peak) distributions. (3) Applying the proper conversion factor for z^peak and CR^peak, we derived single values of flux (F^peak) and luminosity (L^peak). Finally, following the procedures described in Section 5, we derive F_obsc as a function of redshift, flux and luminosity for our sample (Fig. C1). Larger fractions of obscured AGN are in general derived by this procedure, affecting in particular the trends with redshift and flux. This may be at least in part due to an overestimation of the best-fitting column density affecting intrinsically unobscured sources, due to statistical fluctuations in low-quality X-ray spectra discussed in section 3.3 of Vito et al. (2013). When considering the full probability distribution of column density, this issue is greatly alleviated. We also checked that, considering the F^peak and L^peak values, the number counts resulted not to be strongly different, while the XLF tends to be flatter at low luminosities. We again ascribe this behaviour to the tendency of overestimating the column densities, which has the effect of shifting the obscuration-corrected luminosities to larger values.