Space Densities and Emissivities of Active Galactic Nuclei at z > 4

As in G15 the AGN selection from the parent z > 4 catalog is based on the detection of significant X-ray emission in the source position measured in the HST/WFC3 H-band image. In this updated analysis we benefit from the new 7 Ms GDS Chandra image as well as from the shallower 2 Ms GDN and 0.7 Ms EGS Chandra images. The data reduction was done reprocessing all the observations using the Chandra Interactive Analysis Observations (CIAO) software (v4.8; Fruscione et al. 2006) and CALDB v4.8. Intervals of high background were determined by creating 0.3–10 keV background light curves in intervals of 100 s. We rejected time intervals where the background count rate was 3σ above the mean value of the background count rate in the observation. For each observation, we produced energy-filtered events files and exposure maps in several energy bands. We registered and refined X-ray astrometry of each observation adopting a sample of bright pointlike X-ray sources selected from the Chandra catalogs available for the three fields (Alexander et al. 2003; Xue et al. 2011, 2016, and Nandra et al. 2015). The positions of these bright sources have been recovered in each single observation, providing the needed rototranslation. The cleaned and astrometrically corrected event files were then generated for each frame and coadded with a residual positional error of 002. The astrometric solutions found for the full mosaics following this procedure may not be consistent with the CANDELS astrometric solution. Shifts of ≈01–02 are seen between the X-ray and the CANDELS bright sources (see also Luo et al. 2017). Since we are interested in collecting X-ray counts at the position of the CANDELS sources, we realigned the mosaics to the CANDELS reference frames using bright X-ray and optical sources.

We searched for X-ray counts in the final coadded X-ray images at the position of each CANDELS source in the NIR H-band. This is the key to reach the lowest possible flux limit, since analyzing the X-ray emission at the position of known sources allows one to use a less conservative threshold for source detection than in a blind search. We do not correct for possible offsets between the X-ray centoid and the CANDELS position. The X-ray detection strategy and photometry are based on the ephot software and are extensively discussed in Fiore et al. (2012). To reach the faintest X-ray flux limits we chose the (circular) photometric apertures and energy bands that minimize the background and maximize the S/N. We use apertures from 1–7 arcsec (diameter). Because the Chandra point-spread function (PSF) varies strongly with the off-axis angle (and also with the energy band, although less strongly), we use apertures allowing to collect at least 35%–40% of the counts at each off-axis angle. This is important to obtain reliable fluxes, limiting the uncertainty on the PSF aperture correction, which is obviously larger for larger PSF shapes. On the other hand, large apertures may be affected by contamination of foreground sources close to the position of the main target. To limit the contribution of contaminants, we carefully checked all detections based on large apertures and excluded from our sample all the cases where the spurious X-ray flux comes from adjacent bright sources. As we will discuss in the following, this is the main difference with respect to the analysis reported in G15.

The sky coverage and associated incompleteness correction is estimated by using Monte Carlo simulations (see Fiore et al. 2012 for details). The study of faint X-ray source population requires the most careful possible characterization of the background. Following Fiore et al. (2012), we extracted average backgrounds by excluding regions of 10, 15, and 20 arcsec of radii around bright sources. The corresponding background spectra were very similar, and the background obtained with a 10 arcsec exclusion region has been adopted to estimate the background at the position of each CANDELS source. We first extracted spectra at the position of each CANDELS source. We then normalized the average background to the counts detected in the X-ray spectra at the position of the CANDELS sources in the 7–11 keV energy band, where the contribution of the X-ray sources with respect to the internal Chandra ACIS background is negligible (for faint sources) or small (for bright sources), due to the sharp decrease in the Chandra effective area. This procedure allows better count statistics and minimizes systematic errors due to, e.g., source crowding, varying exposure times etc. It also allows a direct estimate of the Poisson probability for background fluctuations at the position of the CANDELS sources.

We estimated the Poisson probability that the counts extracted from a given energy band and a given aperture at the position of each CANDELS source were a fluctuation of the background (estimated as described above). We finally chose the aperture and energy bands producing the smallest probability.

To associate a reliable probability to each X-ray search we resorted again to simulations. We first simulated between 10⁵ and 10⁶ background X-ray spectra at the position of each CANDELS source, to use exactly the same exposure time, vignetting, and PSF. We then applied the ephot detection tool to these simulations and studied the number of spurious sources as a function of the parameters used: (1) the Poisson probability that the simulated counts are indeed a background fluctuation; (2) the size of the source extraction region; (3) the background subtracted counts in the broad 0.3–4 keV band; and (4) the energy bandwidth E_max/E_min. In this multidimensional space we chose a combination of parameters providing a number of spurious detections smaller than one every 5000 spectra. The number of candidates with z > 4 in the three fields analyzed is 4084 (1489 in GDS, 1341 in GDN, and 1254 in the EGS), and we expect about one spurious detection in the overall sample. The estimates of spurious detection probabilities are given in Table 1.

ephot was also run on the galaxy samples with fixed energy bands (0.5–2 keV), thus optimizing only for the size of the source extraction region. The X-ray fluxes in the band 0.5–2 keV were estimated from ephot count rates in the 0.5–2 keV band, after PSF correction, if the S/N in this band is higher than 2.5 or from the count rates in the band, which optimizes the detection otherwise. To convert count rates into fluxes, we used a spectrum with a photon index estimated from the ratio of the counts in the 0.5–2 keV and 2–4 keV bands when a source was detected in both bands, or with a fixed photon index of 1.4 otherwise. To test our photometry, we compared our 0.5–2 keV flux with those published by Luo et al. (2017) for the Chandra Deep Field South, Xue et al. (2016) for the Chandra Deep Field North, and Nandra et al. (2015) for the EGS. The agreement is generally good, with the median of our samples in agreement with those of the published samples within 2% for the CDFS sources, 6% for the CDFN sources, and 4% for the EGS sources.

We have detected in the X-ray images 32 AGN candidates with 4 < z < 6.1 whose positions, with relative astrometric accuracy of ∼1 arcsec, redshifts, and H magnitudes are given in Table 1. This table also lists the results from the X-ray detection procedure: X-ray best detection energy band, photometric aperture (diameter), probability of spurious detection, and 0.5–2 keV flux.

Of the 32 AGN candidates, 19 come from the GDS field (11 in common with Luo et al. 2017), 8 from GDN (6 in common with Xue et al. 2016), and 5 from the EGS field (1 in common with Nandra et al. 2015). The X-ray contours of the 32 sources overlaid with the WFC3 H-band images are shown in the Appendix (Figures 12–14).

Of the 19 AGN candidates in the GDS, 15 are in common with G15. In 4 of the 7 candidates removed from the original G15 catalog (GDS4285, GDS4952, GDS5501, and GDS31334) the X-ray detection remains confirmed at the chosen probability threshold, using a typical detection cell of 2–3 arcsec around the CANDELS sources. However, the comparison of the Chandra map with the HST image suggests that most of the X-ray emission is actually produced by contaminating sources within 1–2 arcsec from the CANDELS high-z target. This contamination was underestimated in G15 because of the shallower X-ray images with respect to the present Chandra 7 Ms images. Thus the X-ray association in the present sample is now more robust. The statistical significance associated with the remaining 3 of the 7 removed candidates (GDS12130, GDS9713, and GDS33073) was just above the threshold in the 4 Ms analysis (G15), but it is just below the threshold in the 7 Ms analysis. This may be due to either source variability or background fluctuation. These 7 sources were also removed by the Parsa et al. (2018) analysis of GDS field. On the other hand there are 4 new candidates at z > 4 in the new 7Ms GDS image, which increase the total number to 19.

We have also produced an X-ray stack of all the new 14 sources not included in the previous X-ray selected catalogs (see Figure 15 in the Appendix). The X-ray stack image in the fixed 0.8–3 keV energy band shows a significant S/N ∼ 10, implying that the bulk of our new sources are X-ray emitters. More details are provided in the Appendix.

Here we remind that the associated X-ray luminosities are in general L_X ≳ 10⁴³ erg s⁻¹ in the 2–10 keV band. These luminosities are more typical of Seyfert-like and QSO sources rather than starburst galaxies for which these high X-ray luminosities would correspond to very high SFRs ≳2000 M_⊙ yr⁻¹ (Ranalli et al. 2003) well above the average star formation rates (SFRs) derived from stellar SED fits to our composite sample. Moreover selecting high-redshift X-ray sources would imply sampling high rest-frame spectral energies >7 keV, where the galaxy contribution is progressively decreasing. We cannot exclude however a significant contribution from stars to the X-ray luminosity of few peculiar sources in our sample.

As in G15 we note that the possible presence of significant X-ray absorption does not imply a similar absorption shortward of the Lyman continuum, because the physical regions responsible for the UV and X-ray emissions are different in size and ionization state. Indeed X-ray absorption is in general close to the emitting region and originated by metals associated to a highly ionized hydrogen.

Most sources in the HST H-band image are compact or too faint for any morphological classification, and only a few spectroscopic redshifts are available in the literature (see Table 1).

Thus, a critical issue for these sources is related to the photometric estimate of redshifts. In spite of the fact that we have adopted a combination of different redshift estimates provided by different authors, the redshifts of a few sources remain uncertain due to their faintness, the almost power-law shape of the SED, and the contamination of their UV/blue ground-based photometry by nearby sources. Few AGN candidates suffer from these large uncertainties, namely, GDS2527, GDS11847, and GDS33160. It is worth noting however that GDS2527 has M₁₄₅₀ > −18, and it is not used for the computation of the z = 4–5 LF. GDS11847 and GDS33160 are at z > 5, and their broad PDF(z) could bias the high-z, LF estimate. We quantify this bias when computing volume densities in Section 3.1.

The SEDs and PDFs(z) of all the AGN candidates together with their X-ray, optical, and IR multiwavelength images are given in the Appendix. In particular, the PDFs(z) show the uncertainties associated with the redshift estimates, which however still depend on the SED templates adopted. From the images it is possible to check the faintness of the z > 4 AGN candidates in the optical bands and the need for deep IR and X-ray detections for any effective multiwavelength selection.

A recent reevaluation of the G15 sample by Parsa et al. (2018) resulted in lower photometric redshifts for a significant number of sources and consequently lower LFs, especially at z > 5. In all the discrepant cases (except one) this is due to low photometric redshifts obtained by means of the adoption by Parsa et al. (2018) of AGN templates coupled with large dust reddening (see also the Appendix). First of all we note that they adopt the Calzetti extinction curve, which is more suitable for starburst galaxies rather than for AGNs, which typically show an SMC extinction curve especially at z < 4 (e.g., Richards et al. 2003). Second, they remove from the fit IRAC data at 5.8 and 8 μm, which are important especially for red sources. This introduces further degeneracy favoring in some cases low-redshift, dusty solutions.

To check the robustness of their criticism we have also included pure AGN templates in our library following the recipe adopted by Hsu et al. (2014) in the extended GOODS-S field, adding possible dust reddening. In Figure 8 in the Appendix we show the difference between the two redshift estimates. The average difference is very small with essentially no offset between the two estimates. We also note there is not significant difference between the average dust reddening derived by redshift best-fit solutions obtained using galaxy or AGN templates (E(B − V) ∼ 0.27 versus E(B − V) ∼ 0.15, respectively). However, there is a fraction of 20% of objects with low-redshift solutions by AGN templates with strong reddening. These solutions however would imply pure AGN SEDs in sources with low luminosities, which in the extreme cases are more typical of dwarf galaxies M_B ∼ −12, −15 without assuming any contamination by the UV light of the host galaxy. We consider unlikely these low-z solutions.

Dusty and reddened AGN templates from the UV to the NIR giving low-redshift solutions appear moreover in contrast with those resulting from recent multicomponent analyses performed on a sample of reddened AGNs by Bongiorno et al. (2014) and LaMassa et al. (2017). Indeed their multicomponent analyses performed on reddened broad-line AGNs with spectroscopic redshifts show optical spectra dominated by the SEDs of the host galaxy. The AGN continuum appears to mainly shape the NIR rest-frame region (Bongiorno et al. 2014; LaMassa et al. 2017) and not the rest-frame optical region as assumed by Parsa et al. (2018). Moreover, the fainter high-redshift AGNs at z ∼ 5–6 seem less dusty with respect to the bright QSOs observed in the same redshift interval. Indeed recent Atacama Large Millimeter/submillimeter Array observations seem to support a scenario where faint AGNs appear to inhabit normal main-sequence or quiescent galaxies (Izumi et al. 2018) in contrast with the brightest QSOs at similar redshifts.

For all these reasons we keep the photometric redshift solutions obtained by a combination of the probability redshift distributions derived by galaxy templates as described above.

We have extracted Multi-Unit Spectroscopic Explorer (MUSE) spectra at the position of our 12 candidates falling in the MUSE-Wide survey in GDS (Herenz et al. 2017; Urrutia et al. 2018). The MUSE-Wide survey is a Guaranteed Time program covering a relatively large area in GDS with a relatively shallow exposure time of 1 hr. MUSE is the integral field unit at the ESO-Very Large Telescope with a field of view of 1 arcmin² covered with a spatial sampling at 0.2 arcsec and a spectral resolution of 2.5 Å from 4750–9350 Å (Bacon et al. 2009). Details on the data reduction analysis are provided in Herenz et al. (2017), where the detection algorithm LSDCat has been used to select emission lines by means of a matched filtering procedure. This analysis has provided the detection of GDS9945 as a clear Lyα emitter at z = 4.50 with a possible AGN signature of weak NV emission. This source was already known as Lyα emitter by means of previous unpublished spectroscopic information (see G15), however since its X-ray position could be contaminated by a close brighter galaxy (see Figure 9 in the Appendix) the source was not used for the computation of the LF. We have also visually checked for weak emission line features, finding another tentative detection in GDS8687 at z = 4.40 in broad agreement with our photometric redshift but in contrast with the Parsa et al. (2018) estimate, which puts GDS8687 at z = 3.57, so out of the z > 4 sample. The emission lines of the two sources are shown in Figure 1. The lines have a S/N ∼ 34, 4, fluxes f ∼ 40, 9 × 10⁻¹⁸ erg s⁻¹ cm⁻², and luminosities L ∼ 7, 1.8 × 10⁴² erg s⁻¹, which are typical of Lyα emitters at z ∼ 5–6. In particular luminosities are confined between 0.1L* and L* at z ∼ 4–4.5 at the faint end of the luminosity interval where AGN activity could be present in Lyα emitters (e.g., Sobral et al. 2019).

Zoom In Zoom Out Reset image size

Although the photometric redshifts have been estimated from galaxy templates, we are assuming that the AGN luminosity is driving the UV SED of our sources with a minor contamination by the host galaxy. To check the validity of our assumption we show in Figure 2 the 2 keV luminosity versus monochromatic UV luminosity at 2500 Å for our 32 sources. Luminosities at 2500 Å have been derived from the observed apparent magnitudes closer to 2500 × (1 + z) Å, namely, the J-band for 4 < z < 5 sources and the H-band for 5 < z < 6. Luminosities at 2 keV have been derived from the 0.3–2 keV flux, assuming a photon spectral index Γ = 1.4 for these faint sources, which is more similar to the background spectrum and takes into account a possible absorption.

Zoom In Zoom Out Reset image size

We have compared the distribution of our sources in the L_X−L_UV plane to the relation measured by Lusso et al. (2010) for the brighter, type 1, COSMOS AGNs. It appears that most of the objects are within 1σ from the correlation without the presence of any significant bias. We conclude that in our sample the AGN contribution to the UV luminosity is on average prevailing on the host galaxy contribution. This is consistent with the SED analysis performed by Bongiorno et al. (2012) on the type 1 COSMOS sample (their Figure 4, right panel), where a two-component model (AGN+galaxy templates) has been adopted. An AGN contribution at least by 50%–60% of the total UV flux has been found within the 25th percentile of the unobscured AGN sample for λ < 3000 Å. Thus our sample in Figure 2 simply represents the fainter end of the COSMOS L_X−L_UV correlation but at higher redshifts.

Also, for the computation of the UV LF in the next section, rest-frame UV 1450 Å absolute magnitudes M₁₄₅₀ are derived directly from apparent magnitudes in those filters with effective rest-frame wavelengths closer to 1450 × (1 + z) Å. This minimizes uncertainties in k-corrections.

Our magnitude-selected galaxy sample in the UV rest-frame band is then used to estimate the 1450 Å LF following G15. The extended version of the standard 1/V_max algorithm (Schmidt 1968) is adopted where different regions with different magnitude limits are combined together in the volume estimate of each object (e.g., Avni & Bahcall 1980).

For a given redshift interval (z_low, z_up), these volumes are confined between z_low and z_lim (j), the latter being defined as the minimum between z_up and the maximum redshift at which the object could have been observed within the magnitude limit of the jth region. Indeed, the complexity of the exposure map in the H-band image of the CANDELS fields requires a composition of several subareas with associated magnitude limits (see G15).

As in G15 the galaxy volume density ϕ(M, z) in a given (Δz, ΔM) bin is

$$\begin{eqnarray}&&\phi (M,z)=\displaystyle \frac{1}{{\rm{\Delta }}M}\displaystyle \sum _{i=1}^{n}{\left[\sum _{j}\omega (j){\int }_{{z}_{\mathrm{low}}}^{{z}_{\mathrm{lim}}(i,j)}\displaystyle \frac{{dV}}{{dz}}{dz}\right]}^{-1},\end{eqnarray} \tag{ 1 }$$

where ω(j) is the area in units of steradians corresponding to the different regions, n is the number of objects in the redshift-magnitude bin, and dV/dz is the comoving volume element per steradian.

The estimates of the UV LFs of AGN candidates selected by their X-ray flux are subject to significant correction for incompleteness, since at faint NIR (rest-frame UV) fluxes only sources that are relatively bright in the X-ray band can be detected even in the deepest GDS field. Figure 3 shows the X/H flux ratio as a function of the H-band magnitude for our AGN candidates in the three fields (GDS9945 left out). The straight lines indicate the locus of the constant X-ray flux limits adopted for the three fields. It is clear, for example, that at H ≳ 26 sources with F_X/(λ f_H) ≲ 0.1 cannot be detected at the X-ray flux limit of the EGS survey. Overall our AGN sample is thus biased against relatively weak X-ray emitters as pointed out in G15. This aspect has been taken into account when computing the faint-end UV LF.

Zoom In Zoom Out Reset image size

The resulting volume densities as a function of M₁₄₅₀, z are shown in Table 2 and Figure 4 in two redshift intervals Δz = 4–5 and Δz = 5–6.1. Volume densities have been first corrected for incompleteness in the H-band galaxy counts at the faintest limits showing, e.g., a 50% drop at H ∼ 27 in GDS and at H ∼ 26.5 in the EGS field. The second correction takes into account the loss of AGN candidates with X-ray flux below the X-ray flux limit for AGNs with a given H magnitude. The incompleteness fraction is derived from the same X/H distribution observed above the X-ray flux threshold for a given H-band flux of each source. The objects have been weighted considering the shape of the corrected H counts. It is clear that the adopted X/H distribution could be biased by selection effects (e.g., volume effects in small area surveys) and could be different from the intrinsic X/H distribution of faint AGNs at z > 4, which is essentially unknown at the redshifts and luminosities probed here. Nevertheless, this is an attempt to estimate any correction for incompleteness due to the X/H properties of the faint AGN population. Again, as in G15, corrections by 10%–20% have also been applied to volume densities due to spatial fluctuations of the X-ray flux limits for each X-ray position.

Zoom In Zoom Out Reset image size

Table 2.  AGN LFs from 1/V_max Analysis

Δz	M1450	ϕobs	ϕcorr	Nobj	ϕMC
4–5	−19	6.81	${14.54}_{-5.81}^{+8.72}$	6	18.05 ± 4.27
	−20	4.74	${11.47}_{-4.59}^{+6.88}$	6	8.03 ± 3.34
	−21	3.29	${5.08}_{-2.21}^{+3.45}$	5	4.52 ± 1.15
	−22	1.24	${1.31}_{-0.87}^{+1.74}$	2	1.33 ± 0.11
5–6.1	−19	3.62	${7.27}_{-4.02}^{+7.12}$	3	6.27 ± 3.42
	−20	3.12	${4.77}_{-2.31}^{+3.79}$	4	2.91 ± 1.84
	−21	0.65	${0.69}_{-0.60}^{+1.61}$	1	1.13 ± 0.70
	−22	0.61	${0.62}_{-0.54}^{+1.44}$	1	0.80 ± 0.33

Note. ϕ_corr is ϕ_obs volume corrected for incompleteness in the X/H distribution. ϕ, ϕ_corr, and ϕ_MC are in units of 10⁻⁶ Mpc⁻³ mag⁻¹. ϕ_MC are average volume densities derived from 1000 simulated catalogs where random photometric redshifts have been extracted from the PDF(z) of each source. See details in Section 3.1.

Download table as:  ASCIITypeset image

Poisson errors in LF bins have been derived adopting the recipe by Gehrels (1986), which is also valid for a small number of sources. Systematic errors as for example due to cosmic variance in small area surveys are expected to be smaller than in G15 due to the three fields used in the present analysis.

The LFs are shown in Table 2 and Figure 4 with (red-filled squares) or without (green crosses) corrections for incompleteness with respect to the expected X/H distribution. They amount to a factor of ∼2 at the faintest magnitude bins. It is to note that the small corrections applied to the brightest points at M₁₄₅₀ ∼ −22 in our sample are probably lower limits since even at these relatively bright magnitudes there could be AGNs not detected in X-ray. We know for example there is a relatively bright AGN (R ∼ 24.8) with spectroscopic redshift at z ∼ 2 in GDN, which lacks X-ray detection in the 2 Ms X-ray image (see, e.g., Figure 2 in Steidel et al. 2002). Figure 4 shows a rather flat LF in the UV luminosity interval typical of the local Seyfert population −21 ≲ M₁₄₅₀ ≲ −18.5 with corrected densities ∼10⁻⁵ Mpc⁻³ mag⁻¹ at z ∼ 4.5.

We note that the number of candidates in the higher redshift bin selected Δz = 5–6.1 is small and prevents any preliminary estimate of a specific shape. In this context the X-ray detection of sources at z ≳ 6.5 is strongly disfavored by sampling the SED at progressively higher rest-frame energies going from 10–16 keV in the redshift interval Δz = 4–7 for observations taken at ∼2 keV. In this range of energies a typical AGN spectrum drops by a factor of ≲2. Thus a typical AGN with a fixed bolometric luminosity (fixed H magnitude), which would be barely detected at 2 keV at z ∼ 4, will be undetected if at z ∼ 7.

As already mentioned, a few AGN candidates suffer from large uncertainties in the estimate of the photometric redshift. Thus to check in general how the uncertainties in the photometric redshifts affect the LF estimates we have extracted randomly by a Monte Carlo technique the photometric redshift for each source according to its PDF(z) shown in Figure 7. In this procedure we have also included AGN candidates at 3 < z < 4, some of which having a significant probability of being at z > 4. These lower z sources have been selected adopting the same procedure used for the z > 4 AGN candidates. In case a spectroscopic redshift is available for an AGN we fix the photometric redshift to the spectroscopic one. We then processed the simulated catalog with the same software adopted for the observed catalog. We did not vary the H-band magnitudes and the weight assigned to each object to take into account the X/optical flux ratio incompleteness. We produced 1000 simulated catalogs used to compute the mean values of the simulated ϕ_MC for the same magnitude bins adopted for the LF derived by the best-fit photometric z. The scatter of the LF has been derived for each absolute magnitude and redshift interval as half of the difference between the 84th and 16th percentiles.

Form Table 2 it is possible to conclude that the resulting ϕ_MC mean values are consistent with the volume densities derived from the best estimates of the photometric redshifts, indicating that the few broad PDF(z) distributions present in our sample do not significantly bias the estimate of the LF. The resulting scatter of the simulated ϕ_MC results is lower than the errors computed according to the Gehrels (1986) recipe, probably due to the small number of sources per bin. Thus the Poissonian errors shown in Table 2 and Figure 4 do not represent a significant underestimate of the true errors.

For these reasons the volume densities based on the redshift best estimates are therefore used in the next sections to estimate the global shape of the LF and its associated emissivity.

Our new corrected volume densities at z ∼ 4.5 are, e.g., a factor of ∼2 lower with respect to those derived in G15 at M₁₄₅₀ = −20. These lower values are due to fluctuations of number densities of AGN candidates in the three fields, to the different incompleteness correction derived from the new X/H distribution and to the cleaning of uncertain X-ray sources, as described in Section 3.1. Our volume densities are consistent with those derived by Parsa et al. (2018) at z < 5 and M₁₄₅₀ ∼ −20, while they are a factor of 3 higher at z = 5.6. This discrepancy is mainly due to the lower redshifts derived by Parsa et al. (2018) for a significant fraction of our GDS sample, probably due to the reddened AGN templates adopted for the estimate of the photometric redshifts, as discussed in the previous section.

The prediction of the ionizing AGN emissivity critically depends on the global shape of the LF and on the escape fraction of ionizing photons from the AGN host galaxy.

In this context a homogeneous UV sample of z > 4 AGNs selected on the basis of the source X-ray detection and extended on a sufficiently large magnitude interval (−29 ≲ M₁₄₅₀ ≲ −18) is not available with the present instrumentation. For this reason we have connected our data points of the LF to that derived by different optical selected surveys under various assumptions about their possible incompleteness.

More specifically, to derive a first guess on the shape of the UV LF we have compared in Figure 4 our volume densities derived at M₁₄₅₀ > −22 with that of the brightest QSOs selected in the SDSS survey, where selection effects with respect to the morphological appearance and X-ray properties are thought to be small. In other words, at the brightest absolute magnitudes sampled by the SDSS survey M₁₄₅₀ < −27 no X-ray QSOs with strong absorption in the rest-frame optical/UV are expected, and the optically selected sample should be representative of the overall AGN population. Volume densities derived at slightly different redshifts have been scaled to our average redshifts adopting a rescaling in normalization. We adopted a redshift decrease Δlog ϕ = −kΔz, where k = 0.34 in the redshift interval z = 4–5 (Schindler et al. 2018) and k = 0.5 at higher z (Fan et al. 2001). At z ∼ 4.5, QSO volume densities derived from the same SDSS sample by different authors differ by more than a factor of 2 at M₁₄₅₀ > −27 due to different procedures adopted to derive the selection function and consequently the correction for incompleteness (e.g., Richards et al. 2006; Fontanot et al. 2007). For this reason we considered only QSO volume densities derived at M₁₄₅₀ < −27.

We have also connected our CANDELS densities with the ones derived by us in the COSMOS field based on an extensive spectroscopic campaign of candidates selected on a multiwavelength criterion (Boutsia et al. 2018). This includes not only the standard color selection but also the use of photometric redshift catalogs and especially X-ray detection. In this respect it is the most similar selection criterion to the one adopted in G15 and in the present work. For this reason the derived volume densities are particularly high at intermediate absolute magnitudes (M₁₄₅₀ ∼ −24), much higher than found by, e.g., Akiyama et al. (2018) in their color-selected QSO sample but closer to the ones derived by Glikman et al. (2011) in their multicolor survey (NOAO Deep Wide-Field Survey).

We have therefore included both Boutsia et al. (2018) and Glikman et al. (2011) in the best-fit analysis of the LF, scaling to z = 4.5 their published densities estimated at z = 3.9 and z = 4.2, respectively. Model 1 includes our COSMOS sample and the SDSS data. Model 2 adds the NOAO sample to the data of model 1.

To check how the inclusion of brighter surveys can alter the faint-end slope of the LF we have first derived a best-fit power-law slope β = 1.74 for our CANDELS data only. The slope is definitely flatter than found at the bright end, suggesting the presence of a significant break at intermediate magnitudes. For this reason we adopted a two power-law shape for the LF of the type

$$\begin{eqnarray}&&\phi =\displaystyle \frac{{\phi }^{* }}{{10}^{0.4({M}_{\mathrm{break}}-M)(\beta -1)}+{10}^{0.4({M}_{\mathrm{break}}-M)(\gamma -1)}}.\end{eqnarray} \tag{ 2 }$$

Two best-fit solutions are shown in Table 3 and Figure 4 (upper panel) at z ∼ 4.5. The first solution (model 1) shows a flat faint-end slope consistent with that derived by the CANDELS data alone. A sharp break is present at intermediate magnitudes (M₁₄₅₀ ∼ −25.8) followed at the bright end by a steep power law γ ∼ 3.7. Such a steep slope at bright magnitudes (M₁₄₅₀ < −27) seems supported by the recent evaluations based on the North sample of the Extremely Luminous Quasar Survey QSO survey at z ≲ 4 (Schindler et al. 2018).

Table 3.  Parametric AGN LFs, Emissivities, and Photoionization Rates^a

Model	Δz	β	γ	Mbreak	log ϕ*	1450	912	Γ
1	4–5	1.70	3.71	−25.81	−6.68	${6.50}_{-3.77}^{+24.26}$	${3.89}_{-2.26}^{+14.52}$	${0.56}_{-0.33}^{+2.09}$
2	4–5	1.74	3.72	−25.89	−6.76	${6.35}_{-3.51}^{+8.53}$	${3.80}_{-2.10}^{+5.10}$	${0.55}_{-0.30}^{+0.73}$
3	5–6.1	1.92	3.09	−25.06	−7.29	1.33	0.80	0.07
4	5–6.1	1.74b	3.72b	−25.37	−7.05	${1.94}_{-1.59}^{+7.37}$	${1.16}_{-0.95}^{+4.41}$	${0.11}_{-0.09}^{+0.41}$

Notes.

^aϕ* in units Mpc⁻³ Mag⁻¹, in units of 10²⁴ erg s⁻¹ Hz⁻¹ Mpc⁻³, Γ in units of 10⁻¹² s⁻¹; model 2 also includes the NOAO data points (Glikman et al. 2011). 1σ errors in β, γ, M_break, and log ϕ* are ∼0.4, 0.7, 1.3, and 1.0 respectively for both models 1 and 2. 1σ errors in M_break, log ϕ* are ∼0.7,0.7 for model 4. β, γ, M_break, and log ϕ* for model 3 are essentially not constrained by the present data. Errors for and Γ are derived computing the 68% joint probability distribution for M_break and log ϕ* having fixed the two slopes. ^bFixed value.

Download table as:  ASCIITypeset image

The second solution (model 2) includes the NOAO survey in the analysis. The resulting densities are lower but closer to our COSMOS data and the LF shape derived from the analysis of the four samples is very similar to the one derived in model 1. Our CANDELS data thus represent the natural extension of the LF after the bending shown in the Boutsia et al. (2018) and Glikman et al. (2011) data.

At z ∼ 5.6 given the poor and uncertain statistics we cannot derive with a similar accuracy a global LF and the uncertainties associated with the LF parameters are much larger. Indeed, in this redshift interval our sample is statistically confined to M₁₄₅₀ > −20 since each of the two brightest bins includes only one object, although the source in the brightest CANDELS bin at M₁₄₅₀ ∼ −22 is a secure AGN (GDN3333) at the spectroscopic redshift z = 5.186. Therefore a best-fit LF derived from the joint analysis of the bright M₁₄₅₀ < −27 SDSS sample and the very faint CANDELS sample resulted in a steeper faint-end slope β ∼ 1.9 and a shallower bright-end slope γ ∼ 3.1, as shown in Table 3 (model 3). If we assume that the shape of the LF does not change appreciably from z = 4.5–5.6 we can fix the two LF slopes found at z = 4.5 in model 2 and adapt normalization and break magnitude to follow the slightly different density evolutions of the bright and faint sides of the LF. For model 4 we derived a break magnitude M₁₄₅₀ ≃ −25.2 and a density decline with respect to the previous redshift bin by factors of 2–5 from the faint to the bright end, respectively. These factors however are sensitive to the uncertain break position of the LF.

The very uncertain model 3 predicts volume densities at intermediate magnitudes −23 ≲ M₁₄₅₀ ≲ −25, which are consistent with or slightly higher than the ones derived from standard optical color-selected QSO surveys by McGreer et al. (2018) at z ∼ 5 and by Onoue et al. (2017) and Matsuoka et al. (2018) at z ∼ 6, all scaled at z = 5.6. Model 4, with its slopes invariance, predicts a number of sources in the same magnitude interval, which is a factor of 3–10 higher. While this appears as a larger factor we note that a factor of 3–4 is also present at z ∼ 4.5 between the Subaru color survey by Akiyama et al. (2018) and the COSMOS spectroscopic sample by Boutsia et al. (2018). This systematic discrepancy between spectroscopic samples of color-selected AGNs and AGNs selected by a broader multiwavelength criterion as in the COSMOS field could depend on the progressive increase of the incompleteness at fainter magnitudes in color and morphological selected AGNs due to various selection effects. For example, in the COSMOS field about half of the X-ray spectroscopically confirmed AGNs lies outside the standard optical color selection. Moreover, the Subaru surveys seem to rely on tighter color selections and appear more conservative to avoid large contamination by Galactic stars; for this reason the final estimated volume densities are the results of a non-straightforward balance between contamination and incompleteness corrections. Thus either a large incompleteness is still present at intermediate absolute magnitudes in optical color-selected, pointlike surveys, or a strong and quick change in the shape of the LF should be envisaged.

In summary, a sharp break of the LF with slope values changing by ∼2 is required at z ∼ 4.5 to follow the volume density decrease from the CANDELS sample down to the COSMOS and SDSS samples. This sharp change implies a major role of AGNs with intermediate luminosity to the ionizing emissivity of the global AGN population up to z ∼ 5 as outlined in the next section. With the present multiwavelength CANDELS data set it is not possible to derive similar constraints in the redshift interval 5 < z < 6.1 and predictions on the global AGN ionizing emissivity rely on extrapolations about the LF shape, as already pointed out in G15.