ACTIVE GALACTIC NUCLEI AS MAIN CONTRIBUTORS TO THE ULTRAVIOLET IONIZING EMISSIVITY AT HIGH REDSHIFTS: PREDICTIONS FROM A Λ-CDM MODEL WITH LINKED AGN/GALAXY EVOLUTION

The assessment of the thermal history of the intergalactic medium (IGM) at high redshifts is fundamental to understanding the physical processes involved in galaxy formation and evolution, including triggering and fueling of active galactic nuclei (AGNs).

The ionization state of the IGM is a function of the hydrogen ionizing UV background (UVB) conceivably produced by sources like high-redshift galaxies and AGNs. Star-forming galaxies are more common in the high-redshift universe, and, in principle, can be responsible for the reionization of the IGM starting at z  >  7. However, even a first-order estimate of their contribution to the ionizing UVB has two main sources of uncertainty. The first is connected with the increasing difficulty in evaluating the faint end slope of the UV luminosity function (LF) especially at the highest redshifts (Bouwens et al. 2011; Grazian et al. 2011). The second is even less understood and is connected with the poor knowledge of the average escape fraction of ionizing Lyman continuum photons from the interstellar medium of each galaxy. On the other hand, the AGN contribution is mainly affected by the poor knowledge of the faint end slope of their LF, especially at high redshifts due to the lack of deep AGN surveys at various wavelengths (see, e.g., Shankar & Mathur 2007). The observed escape fraction of ionizing photons in the majority of the bright AGNs is high and of the order of unity at the Lyman continuum, although the value for the faintest AGNs is unknown.

Since the apparent number density of bright QSOs and AGNs is rapidly decreasing outside the redshift interval 2  <  z  <  3, it is usually assumed that star-forming galaxies can give a considerable contribution to the ionizing flux that could become dominant at z  <  1 and at z  >  3 (Madau 1991; Giallongo et al. 1997; Haardt & Madau 1996, 2012). In this context several attempts have been made to derive the escape fraction of UV ionizing photons both at low redshifts from space and at high redshifts z = 2–4 from ground-based observations. However, there is little evidence supporting a scenario where enough ionizing photons escape from galaxies at intermediate redshifts, and most attempts gave only upper limits on the escape fractions.

Indeed, recent measurements by Cowie et al. (2009) provided a stringent 2σ upper limit of 0.8% for the relative escape fraction of galaxies at 0.9  <  z  <  1.4, suggesting that galaxies are not the dominant contributors to the ionizing UVB at z  <  2. At high redshifts the situation is more difficult from the observational point of view. Direct measures of the Lyman continuum flux from galaxies at z  >  4 are increasingly difficult because of the sharp increase in the number density of intervening Lyman limit absorptions due to optically thick clouds in the IGM. For this reason the search has been focused on the range z ∼ 3–4, but a discrepancy is found there by different teams using different methods. Steidel and collaborators support the detection of appreciable Lyman limit flux escaping at least from the UV brightest galaxies at z ∼ 3 (Steidel et al. 2001). More recent analyses suggest average escape fractions of the order of 10%–15% (Shapley et al. 2006; Nestor et al. 2011), although their sample can be contaminated by foreground objects (Vanzella et al. 2010b; Boutsia et al. 2011). Other teams using both spectroscopic and very deep broad or narrowband imaging gave only upper limits in the range <5%–15% (Giallongo et al. 2002; Fernandez-Soto et al. 2003; Vanzella et al. 2010a; Boutsia et al. 2011). To summarize, the recent trends appear to favor low escape fractions from relatively bright galaxies, where the recent limits are below 1% at low redshifts and below 10% at z ∼ 3–4, and cast serious doubts on any redshift evolution in the galaxy ionizing escape fraction. This leaves open the possibility that AGNs could still play a leading role in contributing to the cosmological ionizing background.

The traditional view based on relatively bright QSO surveys in optical and X-ray bands leads to the standard view of an evolutionary ionizing QSO background increasing from the local value out to z = 2–3 and then quickly decreasing (see e.g., Madau et al. 1999). Previous estimates were based on simple parametric extrapolations both in luminosity and redshifts of the AGN LF without any link to specific physical processes (Madau et al. 1999). However, recent deep optical QSO surveys at z = 3 (Siana et al. 2008; Fontanot et al. 2007) and z = 4 (Glikman et al. 2011) are showing the presence of a considerable number of previously unknown faint AGNs producing a rather steep LF. The inclusion of X-ray detection in the selection methods enables us to extend the knowledge of the LF to even fainter limits (Fiore et al. 2012) of the order of L(2–10 keV) >10⁴³ erg s⁻¹ out to z  >  4.

The presence of this faint population is changing our estimate of the AGN contribution to the ionizing UVB, although the selection of faint AGNs at the highest redshifts z ∼ 6 becomes difficult with the current instrumentation. Given the current observational limits, an insight into the origin of the ionizing flux at high redshifts can be gained from theoretical modeling. The latter should account for the formation and evolution of galaxies in a cosmological context, and for the co-evolution of the AGN population. Given the complex sub-grid physics involved in the AGN activity over a wide range of cosmic times and of galactic masses, semi-analytic models (SAMs) constitute a powerful tool to provide a statistically relevant sample of simulated galaxies and AGNs. Such models connect the physical processes involving gas and star formation to the merging histories of dark matter (DM) halos collapsed from the primordial density field. While several SAMs, including the growth of supermassive black holes (SMBHs), have been presented in the literature (starting from Kauffmann & Haenhelt 2000; see also Croton et al. 2006; Bower et al. 2006; Marulli et al. 2008; Menci et al. 2006, 2008), addressing the problem of the ionizing background requires a physical description of the escape fraction in galaxies. Such a problem is deeply connected with that concerning the feedback of the AGNs onto the interstellar gas. In fact, since even a tiny amount (column densities N_H ∼ 10¹⁷ cm⁻²) of galactic gas would suppress the radiation at the Lyman limit escaping from the galaxy, estimates of the escape fraction translate into computing the fraction of photons that are emitted along the directions where the interstellar gas of the galaxy has been depleted. In our previous work (Menci et al. 2008), we have developed our SAM model (Rome-SAM, hereafter R-SAM)) to describe the expansion of the blast wave produced in the galaxy interstellar medium by AGN-driven outflows. Here we use the same model to derive the escape fraction of photons in active galaxies as the fraction of photons emitted along the directions where the AGN-induced blast wave has expelled the galactic gas. When coupled with the integrated LF of AGN self-consistently derived in the same model, this allows us to compute the ionizing background produced by AGNs within a cosmological model of galaxy formation. We stress that the basic quantities predicted by our computation have been extensively tested against observations in our previous works. On one hand, the predictions concerning the growth of SMBHs and the corresponding evolution of the LF have been compared against a wide range of different observations, from optical to X-rays, for SMBHs masses and AGN luminosities in a wide range of cosmic times extending up to z = 6 for the brightest QSOs (Menci et al. 2003, 2004, 2006, 2008; see also Fiore et al. (2012)); on the other hand, the blast-wave model for the AGN feedback allowed us to derive—within a full cosmological context and in connection with the galaxy properties—the distributions of column densities in active galaxies, as a function of redshift and of the AGN luminosity, which we successfully tested against observations (Menci et al. 2008). In this paper we extend the analysis of the UV AGN LF to faint objects and adopt the blast-wave mechanism to predict the amount of UV photons escaping from the AGN host galaxies.

The paper is organized as follows. In Section 2.1 we recall the basic properties of our SAM for the evolution of galaxies and AGNs. In Section 2.2 we describe how our blast-wave model for the AGN feedback (developed in previous papers) allows us to compute the luminosity-dependent escape fraction of ionizing photons emitted by AGNs. The evolution of the LF of AGNs predicted by our model is compared with data in Section 3; the comparison with the current data set is discussed taking into consideration possible biases in the model and in the data. In Section 4 we convolve the predicted intrinsic LF of the AGN with the escape fraction derived in our model to derive the emissivity and the photoionization rate predicted by the model at various redshifts, and compare the latter with existing data; the contribution of AGNs with different luminosities to such observables is also shown and discussed. Section 5 is devoted to a summary and to conclusions.

The R-SAM we develop and use connects the AGN evolution with the evolution of galaxies in a cosmological framework (see Menci et al. 2008 and references therein).

2.1. Hierarchical Galaxy/AGN Evolution in the R-SAM

Galaxy formation and evolution is driven by the collapse and growth of DM halos, which originate by gravitational instability of overdense regions in the primordial DM density field. This is taken to be a random, Gaussian density field with cold dark matter power spectrum expected within the "concordance cosmology" (Spergel et al. 2007) for which we adopt round parameters Ω_Λ = 0.7, Ω₀ = 0.3, baryonic density Ω_b = 0.04, and Hubble constant h = 0.7 in units of 100 km s⁻¹ Mpc⁻¹. The normalization of the spectrum is taken to be σ₈ = 0.9 in terms of the variance of the field smoothed over regions of 8 h⁻¹ Mpc. The merging histories of DM halos and sub-halos are followed through a Monte Carlo simulation of the collapse and subsequent merging history of the peaks of the primordial density field, which enables us to generate a synthetic catalog of model galaxies and of their past merging history. The physics of baryons (gas radiative cooling, disk formation, quiescent and impulsive star formation, stellar and AGN feedback) is connected to the merging histories of DM halos as described in detail in our previous papers (e.g., Menci et al. 2005). Note that the impulsive star formation implemented in the model is triggered not only by (major and minor) merging events but also by fly-by events between galaxies. The R-SAM also includes a treatment of SMBHs growing at the center of galaxies by interaction-triggered inflow of cold gas, following the physical model proposed by Cavaliere & Vittorini (2000) and implemented in Menci et al. (2003, 2006, 2008). The accretion of cold gas is triggered by galaxy encounters (both of fly-by and of merging kind), which destabilize part of the galactic cold gas mass m_c by inducing loss of angular momentum. The fraction of cold gas accreted by the BH in an interaction event is computed in terms of the variation Δj of the specific angular momentum j ≈ Gm/v_d of the gas, to read (Menci et al. 2003) f_acc ≈ 10⁻¹ |Δj/j|. For minor merging events and for the encounters among galaxies with very unequal mass ratios m'/m ≪ 1 that dominate the statistics in all hierarchical models of galaxy formation, the accreted fraction takes on values 10⁻³ ≲ f_acc ≲ 10⁻². The average amount of cold gas accreted during an accretion episode is thus Δm_acc = f_acc m_c, and the duration of an accretion episode, which provides the timescale for the QSO or AGN to shine, is assumed to be the crossing time τ of the destabilized cold gas component. The time-averaged bolometric luminosity so produced by a QSO hosted in a given galaxy is then provided by

$$\begin{equation} L={\eta \,c^2\Delta m_{{\rm acc}}\over \tau }. \end{equation} \tag{ 1 }$$

We adopt an energy-conversion efficiency η = 0.1 (see Yu & Tremaine 2002), and derive luminosities in the various bands on adopting standard spectral energy distributions as in Marconi et al. (2004) and A. Marconi et al. (in preparation). The SMBH mass m_BH grows mainly through accretion episodes as described above, besides coalescence with other SMBHs during galaxy merging. As an initial condition, we assume small seed BHs of mass 10² M_☉ (Madau & Rees 2001) to be initially present in all galaxy progenitors; our results are insensitive to the specific value as long as it is smaller than some 10⁵ M_☉. Due to the crucial role of AGN feedback in affecting the star formation history of massive galaxies, this part of the model has been extensively tested. The predicted QSO LF and X-ray LF and emissivity have been checked against observation over a wide range of redshift 0  <  z  <  5 and bolometric luminosities 10⁴³  <  L  <  10⁴⁶ erg s⁻¹ (Menci et al. 2003, 2004). The consistency of the modeling of the AGN feedback with the color distribution of galaxies, including a description of the bimodal appearance, has been discussed in Menci et al. (2006).

2.2. Blast-wave Model for AGN Feedback and the Escape Fraction of UV Ionizing Radiation

Fast winds with velocities up to 10⁻¹c are observed in the central regions of AGNs. They likely originate from the acceleration of disk outflows by the AGN radiation field (see Begelman 2003 for a review), and affect the environment in the host galaxy and beyond, leaving imprints out to large scales of some 10² kpc in the intracluster medium (ICM). Clear examples of outflows in local AGNs are those observed on galactic scales in the infrared and optical bands in Mrk 231 and Mrk 573 (Feruglio et al. 2010; Fischer et al. 2010; Rupke & Veilleux 2011). A detailed model for the transport of energy from the inner, outflow region to larger scales has been developed by Cavaliere et al. (2002) and Lapi et al. (2005). Central, highly supersonic outflows compress the gas into a blast wave terminated by a leading shock front, which moves outward with a lower but still supersonic speed and sweeps out the surrounding medium. The key quantity determining all shock properties is the total energy ΔE injected by AGNs into the surrounding gas. This is computed as

$$\begin{equation} \Delta E=\epsilon _{{\rm AGN}} \eta c^2 \Delta m_{{\rm acc}}=\epsilon _{{\rm AGN}}\,L\,\tau \end{equation} \tag{ 2 }$$

for each SMBH accretion episode in our Monte Carlo simulation; the value of the energy feedback efficiency for coupling with the surrounding gas is taken as _AGN = 5 × 10⁻², which is consistent with the values required to match the X-ray properties of the ICM in clusters of galaxies (see Cavaliere et al. 2002). This is also consistent with the observations of wind speeds up to v_w ∼ 0.1c in the central regions, which yield _AGN ∼ v_w/2c ∼ 0.05 by momentum conservation between photons and particles (see Chartas et al. 2002; Pounds et al. 2003); this value has also been adopted in a number of simulations (e.g., Di Matteo et al. 2005) and in our R-SAM (e.g., Menci et al. 2006).

The blast expands into the ISM or the ICM as described by hydrodynamical equations that include the effects of not only an initial density gradient, but also those of upstream pressure and DM gravity, clearly important quantities in galaxies, as discussed by Lapi et al. (2005). The solutions show in detail how the perturbed gas is confined to an expanding shell bounded by an outer shock at radius R_s(t), which sweeps out the gas surrounding the AGN. An analytic expression for R_s(t) has been derived by the above authors in the case of shock expansion in a gas with power-law density profile ρ∝r^−w; in the relevant central regions of the galaxies (within one-half the scale length r_d of the exponential disk) the expansion of the shock is given by (see Menci et al. 2008)

$$\begin{equation} R_s(t)\approx r_d\,{\mathcal M} \,t/t_d, \end{equation} \tag{ 3 }$$

where ${\mathcal M}$ is the Mach number. This is related to the energy ΔE (Equation (2)) injected by the AGN (relative to the initial thermal energy content E∝m_c of the galactic gas) by the relation ${\mathcal M}^2=1+\Delta E / E$ (Lapi et al. 2005). Thus the expansion velocity of the blast wave is directly related to the AGN luminosity.

The above blast-wave model for the AGN feedback allows us to self-consistently compute the escape fraction of ionizing photons. In fact, the typical galactic column densities are in the range N_H ≈ 10²⁰÷10²⁴ cm⁻² and effectively suppress the escape of ionizing radiation along the lines of sight where the galactic gas has not been swept out by the blast wave produced by the AGN. Thus, the escape fraction can be computed as f(L) = ω/2 π, where ω is the solid angle subtending the intercept of the expanding blast wave with the galactic disk. This is illustrated in Figure 1, which shows how only ionizing photons propagating along the directions (indicated with ℓ₁ in the figure) within the solid angle ω can travel unabsorbed and escape the galactic disk, while directions outside the solid angle ω (indicated with ℓ₂ in the figure) intercept a fraction of the gas disk, thus blocking the escape of ionizing photons. Note that the solid angle ω increases as the shock radius R(t) propagates outward.

Zoom In Zoom Out Reset image size

We can compute the solid angle ω = 2π (1 − cos  θ/2) in terms of the opening angle θ (the aperture of the inner cones in Figure 1), which in turn is related to the vertical scale height h of the disk and to the radius of the expanding shell R(t) by the relation cos  θ/2 = h/R(t). Thus we obtain the time-dependent escape fraction f(L, t) = ω/2π = 1 − h/R(t), where the shock radius R(t) depends on the AGN luminosity (see Equation (3) and text below). Since we are interested in the statistical, effective escape fraction, we compute the average over the duration τ of the AGN activity

$$\begin{equation} \langle f(L)\rangle ={1\over \tau }\,\int _{t_S}^{\tau }\,dt\,\Bigg [1-{h\over R(t)}\Bigg ]. \end{equation} \tag{ 4 }$$

Here t_S is the time (computed from the start of the AGN phase) at which the shock radius first encompasses the width of the disk (for earlier times escape fraction is null), so that R(t_S) = h. Using Equation (3) we express t_S as a function of the disk dynamical time t_d as $\tilde{t}_S\equiv t_S/t_d=(h/r_d)/\mathcal M$. Performing the integral in Equation (4) then yields

$$\begin{equation} \langle f(L)\rangle =1- { \tilde{t}_S \over \tilde{\tau } }\Bigg [1+ { \ln {\tilde{t}_S \over \tilde{\tau }}}\Bigg ], \end{equation} \tag{ 5 }$$

where $\tilde{\tau }=\tau /t_d\approx 1$ (for interaction-driven AGNs). Note that for shocks with $\mathcal M\gg 1$ (i.e., large AGN-injected energies ΔE ≫ E, see Equations (2) and (3)) imply $\tilde{t}_S\ll \tilde{\tau }$, yielding large average escape fractions 〈f(L)〉 ≈ 1, while slowly expanding shocks yield 〈f(L)〉 ≈ 0. Thus Equation (5) results in an effective escape fraction increasing with the AGN luminosity L. At fixed L the escape fractions in Equation (5) increase with redshift since—at given ΔE—the smaller host galaxy masses yield larger ΔE/E at higher z.

For a ratio h/r_d = 1/15 (see Narayan & Jog 2002) determining the timescale t_S for each AGN luminosity (corresponding to h ≈ 200 pc for a typical L galaxy), we obtain low-redshift (z ≃ 0.5) escape fractions increasing from 〈f〉 ≈ 0.1 for M₁₄₅₀ = −22 to 〈f〉 ≈ 0.4 at M₁₄₅₀ = −24, while at higher redshifts (z ≃ 4) the values range from 〈f〉 ≈ 0.3 for M₁₄₅₀ = −22 to 〈f〉 ≈ 0.8 at M₁₄₅₀ = −24; quantitative information on the dependence of the escape fraction on the AGN luminosity and redshift can be found in Figure 2. Note that the value of the parameter h/r_d = 1/15 adopted here to compute 〈f〉 from Equation (5) is the same adopted in Menci et al. (2008), where we showed that the blast-wave model for AGN feedback provides column densities N_H (computed after averaging over all the line of sights schematically shown in Figure 1) depending inversely on the AGN luminosity. Such behavior is due to the fact faster shocks correspond to a lower fraction of still unperturbed gas outward of r = R_s(t), and is in quantitative agreement with observations, as we showed in the paper mentioned above.

Zoom In Zoom Out Reset image size

To predict the UV ionizing emissivity produced by our AGN evolution model, we first evaluate the evolution of the AGN LF as a function of redshift at λ = 1450 Å, where most of the UV rest frame data on AGNs are collected and where the average AGN emission shows the peak in the spectral energy distribution, representative of the overall bolometric emission. The contribution of each luminosity to the ionizing background will be the computed by convolving such LFs with the escape fraction in Equation (5).

In Figure 2 we show the evolution of the UV LF N(L₁₄₅₀) computed in our R-SAM (binned in magnitude, dashed line), and the fraction N(L₁₄₅₀) 〈f〉 that will contribute to the ionizing background (solid line), where 〈f〉 is the luminosity-dependent escape fraction given in Equation (5). The bolometric luminosity L (computed for each AGN in our R-SAM after Equation (2)) is related to the emission L₁₄₅₀ at 1450 Å after the bolometric correction given in Marconi et al. (2004) and A. Marconi et al. (in preparation).

In the figure we compare our predictions with data from literature. Note that the data should be compared with the global UV LF (dashed line). At redshifts z ⩽ 2.5 we have compared our predictions with the LF provided by Richards et al. (2009) and Croom et al. (2009). The data come from the combined 2dF-SDSS sample of quasars photometrically selected from Sloan Digital Sky Survey, and then spectroscopically confirmed using the 2DF instrument on the Anglo Australian Telescope. The LFs are complete down to M₁₄₅₀ = −21.5 and show a gradual flattening below M₁₄₅₀ = −24. The shaded region brackets the determinations derived in the X-ray band by La Franca et al. (2005), Ebrero et al. (2009), and Aird et al. (2010), and thus accounts for systematic errors only. Conversions to 1450 Å have been obtained using the Marconi et al. (2004) relation. Note that La Franca et al. (2005) and Ebrero et al. (2009) determinations account for the obscured sources whose flux is reduced below the flux limit by intervening column densities, they provide densities slightly higher than those of Aird et al. (2010), which are not corrected for this effect.

In this redshift range, the agreement with the model predictions is good, if we consider that the errors shown are only Poissonian and that various systematics (e.g., in the flux corrections) probably would widen the uncertainties involved. In any case, the LFs derived from X-ray surveys (see shaded regions in Figure 2) constitute upper boundaries that are fully consistent with our predictions.

At higher redshifts z  >  3 we compare with the rest-frame UV LFs from Fontanot et al. (2007), Bongiorno et al. (2007), Siana et al. (2008), and Glikman et al. (2011). Here the comparison with data at the faint end is more critical, since the corrections for incompleteness of the observations play a more relevant role. For example, the data by Siana et al. (2008) have been derived combining IR SPITZER SWIRE data and optical through a standard morphological criterion and a color selection based on the UV dropout expected at these redshifts coupled with an IR color selection aimed at removing low-redshift interlopers. These authors attain a completeness level at the faintest survey limits of about 75% which is taken into account in their effective volumes. However, their completeness estimates are based on detailed simulations involving both templates of quasar spectra and a model for IGM absorption, which is crucial for any UV dropout selection. In their analysis they refer to the Madau (1995) model, whose analysis has been superseded since larger statistical information about the IGM absorption became available (e.g., Meiksin 2006; Inoue & Iwata 2008), so that now the Madau model represents an upper limit to the average absorption. More recent and reliable estimates by Prochaska et al. (2009) show average Lyman limit opacities by the IGM observed in 1800 Sloan quasars at z ∼ 3.5–4.5 about a factor of two lower. This should allow the presence of a significant number of z ∼ 3 quasars bluer than expected by the UV dropout threshold adopted by Siana et al. (2008), implying a lower completeness level of the LF difficult to quantify. A similar bias could be present in the Fontanot et al. (2007) data which go two magnitudes fainter at a slightly higher average redshift, z ∼ 3.7. The LF has been derived from combined standard color selection and X-ray matching in the GOODS fields. The color selection was effective in the interval 3.5  <  z  <  5.2. In Figure 2 we show the LF derived in the 3.5  <  z  <  4 interval where the completeness of the survey is higher. Even in this survey the average IGM absorption adopted to derive the color threshold is similar to that predicted by the Madau model at least up to z ∼ 4.5, although they apply some correction factor to the Lyα forest. The important point in this respect is that at z  >  3 the incidence of IGM absorption systems which dominate the Lyman limit opacity has a slower evolution compared to the Lyα lines (Prochaska et al. 2009). A further important source of incompleteness at magnitudes M₁₄₅₀ ≳ −24 is constituted by the morphological selection criterion adopted in the spectroscopic follow-up to reduce the contamination from Lyman break galaxies at the same redshift as described in Glikman et al. (2010, 2011).

In Figure 2 we include two faint points derived from a recent very deep search for AGNs in the GOODS south field by Fiore et al. (2012). The latter authors adopt different criteria starting from the standard color dropout technique but including faint detection in the 4 Ms Chandra images in the 2–10 keV X-ray band and relaxing in this way any morphological restriction. They use spectroscopic redshifts where available or accurate photometric redshifts from the GOODS MUSIC and GOODS ERS multicolor catalogs (Grazian et al. 2006, 2011). In this way they can select faint AGN activity among the numerous Lyman break galaxies in a wide redshift interval 3  <  z  <  7 putting interesting constraints on the faint end of the LF at z = 3–4, 4–5, and z  >  5.8 dominated by Seyfert-like objects. Although a fraction of the order of 20% of the detected AGNs appears strongly absorbed in X-ray or Compton thick, the same objects do not show any reddening in the optical–UV part of the spectrum which shows emission lines with various widths. Indeed, a subsequent spectroscopic optical follow-up revealed even in the case of Compton thick sources, typical highly ionized UV emission lines like C iv and N v. In summary, the limited spectroscopic information available coupled with the analysis of the overall spectral energy distributions suggest the absence of any significant link between the presence of strong absorption in X-ray and emission of UV photons from the central nucleus, at least in cases where the X-ray flux is measured for optically UV-selected AGNs. In this respect the sample is more complete with respect to the optical surveys and this is reflected in the higher densities attained by the Fiore et al. (2012) LFs. Nevertheless, the sample could be affected by several uncertainties and significant incompleteness, especially at the faintest limits where the objects have been selected from very deep Hubble Space Telescope NIR images at H ∼ 27 and measured at the faintest X-ray fluxes of F_X ∼ 2 × 10⁻¹⁷ erg s⁻¹ cm⁻². At the lowest X-ray luminosities, objects with relatively high X/optical ratio, and consequently with H > 27, are missed in their survey. Most of the redshifts have been evaluated photometrically and this could introduce low-redshift interlopers in the sample. On the other hand, even in this case the color selection based on the standard UV dropout technique adopts a predicted IGM opacity higher than measured by Prochaska et al. (2009) although lower than expected in the Madau model, allowing for some incompleteness. As a consequence, uncertainties in the Fiore et al. LF data due to systematic errors could be larger than shown in the figure based on number statistics. Keeping in mind these limitations we have adapted their LF data after having recast the X-ray luminosities in UV (1450 Å) magnitudes following the standard Marconi et al. (2004) corrections. The same standard corrections have been applied to shallower X-ray data from the COSMOS field by Brusa et al. (2010) and Civano et al. (2011).

At redshifts z  >  4 the above observational biases and possible incompleteness could explain the model over prediction appearing at fainter magnitudes M₁₄₅₀ > −24 respect to the present data set. The discrepancy is a factor of 2–3 in the magnitude interval −24 < M₁₄₅₀ < −22 up to z ∼ 6 and increases to fainter magnitudes. Since the estimate of the observational incompleteness is difficult to assess without the data input from deep X-ray surveys, we shall rely on our model to get an insight on the contribution of AGNs to the ionizing background. To this aim, we shall couple our predicted LFs with the physical description of the escape fraction based on the blast-wave model.

We now proceed to compute the contribution of our predicted LFs to the AGN emissivity at the Lyman limit. This is given by

$$\begin{equation} \epsilon =\int \langle f(L)\rangle N(L)\,L\,\,{E_{912}\over E}\,{\it dL}, \end{equation} \tag{ 6 }$$

where the 〈f(L)〉 is the fraction of AGNs of a given luminosity which shows Lyman continuum emission given by Equation (5), and E/E₉₁₂ is the (luminosity dependent) bolometric correction at 912 Å. In the integral we have considered all the AGNs produced in our Monte Carlo code down to very faint magnitudes; however, the total emissivity at a given redshift is independent of the faintest luminosity cutoff. This is clearly shown in Figure 3, which shows the contribution to the total emissivity from AGNs with different bolometric luminosity L at different redshifts. Since the bolometric correction from the UV M₁₄₅₀ magnitude is not linear, we have indicated on top of the first panel the UV magnitudes at 1450 Å corresponding to the peak bolometric luminosities. Note that while at z ∼ 3 the main contribution comes from AGNs in the relatively wide magnitude interval −24 < M₁₄₅₀ < −18.5, at higher redshifts the main contribution comes from brighter AGNs with −24 < M₁₄₅₀ < −22. In this magnitude interval the discrepancy between data and model prediction is within a factor of two to three, which is of the same order as the differences among surveys where AGNs are selected with different methods in various bands.

Zoom In Zoom Out Reset image size

In this respect it is just in this magnitude interval that a more careful analysis should be performed to investigate whether data and model prediction can be reconciled. Unfortunately, this magnitude interval is sparsely sampled by present surveys; more statistics with much more controlled systematics should be produced in the future. In particular, multiwavelength surveys of medium depth (e.g., I ∼ 22) in an adequate sky area should be planned for spectroscopic and X-ray follow-up. Meanwhile, we explore the ability of our linked AGN/galaxy model in providing sufficient UV photons to keep the IGM ionized up to z ∼ 6. In Figure 4 and Table 1 the predicted ionizing UV emissivity ₂₄ in units of 10²⁴ erg s⁻¹ Hz⁻¹ Mpc⁻³ is shown as a function of redshift. The sharp increase up to z ∼ 2 is due to the luminosity evolution of the AGN LF. At z  >  4.5 the emissivity gradually declines up to z ∼ 8 and then drops quickly. The rather flat behavior at high redshifts is due to the frequent galaxy merging rate which drives the AGN activity coupled with the short gas cooling time. Both guarantee a high AGN duty cycle and consequently an average AGN emissivity weakly evolving in time.

Zoom In Zoom Out Reset image size

Table 1. Lyman Continuum Emissivity as a Function of Redshift Produced by the AGN Population as Shown in Figure 4

Redshift	log 24	Redshift	log 24
0.25	−0.66	4.75	1.16
0.50	−0.23	5.00	1.14
0.75	0.06	5.25	1.12
1.00	0.30	5.50	1.10
1.25	0.46	5.75	1.08
1.50	0.59	6.00	1.05
1.75	0.66	6.25	1.03
2.00	0.78	6.50	1.00
2.25	0.86	6.75	0.97
2.50	0.91	7.00	0.94
2.75	0.96	7.25	0.91
3.00	1.01	7.50	0.87
3.25	1.06	7.75	0.81
3.50	1.10	8.00	0.75
3.75	1.14	8.25	0.69
4.00	1.17	8.50	0.58
4.25	1.19	8.75	0.46
4.50	1.18	9.00	0.32

Download table as:  ASCIITypeset image

The photoionization rate per hydrogen atom Γ₋₁₂ in units of 10⁻¹² s⁻¹ is the standard quantity used to evaluate the ionization status of the IGM. We have computed this quantity following Madau et al. (1999) and Schirber & Bullock (2003)

$$\begin{eqnarray} &&\Gamma _{-12}(z)\simeq 0.46 {\epsilon _{24}(z)\over 3+{\mid} \alpha _{{\rm UV}}{\mid} } \left(\Delta l\over 50\ {\rm Mpc}\right){\left(1+z\over 4.5\right)}^{3-(2.5+\gamma)}, \nonumber \\ \end{eqnarray} \tag{ 7 }$$

where α_UV = −1.76 is the average UV spectral index of AGNs derived from the average SED we adopt (Marconi et al. 2004). The power-law index −(2.5 + γ) describes the decrease in redshift of the mean free path (mfp) Δl of ionizing photons in the IGM due to the increase in redshift of the Lyman limit absorption systems with slope γ. We adopt γ = 1.94 and the normalization of the mfp to 50 Mpc at z = 3.5 as in Songaila & Cowie (2010). At z  >  4.2 we adopt a steeper redshift evolution with γ = 5 to fit the mfp data point at z ≃ 5.8 (see Songaila & Cowie 2010 for details). This equation represents a first-order approximation where only ionizing sources within one absorption length (where τ_IGM = 1) contribute to the ionizing background. The resulting curve from our Monte Carlo model is shown in Figure 5. The dotted one represents the intrinsic production of UV ionizing photons by the overall AGN population. The ratio between the two curves depends on the ratio between the Lyman limit optically thin AGNs and the overall population.

Zoom In Zoom Out Reset image size

Data points with different symbols are also shown for comparison. They represent different estimates derived in a model-dependent way from the IGM statistics. Indeed most values come from methods of estimating the UVB; and so Γ based on the mean flux decrement in the Lyα forest of QSO spectra in combination with numerical simulations. The UVB value is adjusted until the mean flux in simulated Lyα forest spectra is equal to that in real data (Faucher-Giguere et al. 2008; Bolton et al. 2005; Becker et al. 2007; Fan et al. 2006; Wyithe & Bolton 2011). An alternative method is based on the proximity effect which has recently been proposed at the highest redshifts z = 5–6 by Calverley et al. (2011). However, both methods are subject to uncertainties and systematics which are responsible for the scatter among the data points reported in Figure 5. The data shown in the figure derived from flux decrement analysis in the IGM Lyman forest of QSO spectra, have been scaled to the same IGM temperature–density relation as in Calverley et al. (2011).

Indeed, converting the mean flux decrement into an ionization rate depends on modeling the gas density and temperature distribution in the IGM at very low densities (Miralda-Escudé et al. 2000) in various cosmological scenarios. On the other hand, estimates from the proximity effect are sensitive to the assumption concerning the presence of overdensities of matter close to the measured quasar and to the assumption that gas temperature within the proximity region is similar to that in the general IGM (Calverley et al. 2011).

The continuous curve represents the hydrogen ionization rate as a function of redshift predicted by our model, where the oscillating behavior is related to the stochastic nature of the AGN trigger, and rendered by our Monte Carlo code. The model is able to reach an almost constant value of the order of Γ₋₁₂ ≃ 1 consistent with the data points, in a broad redshift interval from z ∼ 2 up to z ∼ 4. A decline by a factor of four is then appearing at higher redshifts up to z ∼ 6 mainly due to the significant decrease of the mfp of ionizing photons.

The gradual decline at high redshifts is a strong intrinsic prediction of our model based on a linked AGN/galaxy evolution and is in line with recent data analyses at z ∼ 5–6 (Calverley et al. 2011; Wyithe & Bolton 2011). As we have shown in the previous section, the model overestimates the number of relatively (M₁₄₅₀ ∼ −23) faint AGNs at z  >  4 by a factor of ∼2–3; on the other hand, we argued the possible presence of significant incompleteness in the scanty present faint data set.

Since the predicted emissivity by hard spectrum sources like AGNs keeps high up to z ∼ 6 it is interesting to evaluate the impact of this evolution on the ionization history of He ii in the IGM. Indeed, the higher He ii ionization threshold respect to hydrogen and its smaller photoionization cross section imply a complete He ii reionization only when numerous hard spectrum sources like AGNs are distributed within the IGM. There is some consensus about a late He ii reionization which could complete at z ∼ 3 as suggested by observations of He ii Lyα Gunn–Peterson trough and patchy absorption in the He ii Lyα forest at z ∼ 3 (Shull et al. 2010 and references therein). Other indirect evidence comes, e.g., from an apparent increase of the average IGM temperature at z ∼ 3 measured from metal line ratios (e.g., Songaila 1998).

Following early analysis by Madau et al. (1999) and more recently by Haardt & Madau (2012) we have computed the volume filling factor of He iii in the IGM as a function of redshift. We have also considered the case where Lyman limit absorption by optically thick clouds (LLSs) limits the AGN ionization, as proposed by Bolton et al. (2009). In the latter case each AGN is able to ionize only a limited region due to the fact that, at variance with hydrogen ionization, the distance between He ii Lyman limit optically thick clouds is smaller than the typical AGN separation at z ≳ 3. Thus, He ii-ionizing photons will have intersected several helium LLSs being subject to absorption. The equation describing the evolution of the He iii volume filling factor Q is, according to Bolton et al. (2009): $\dot{Q} =\dot{n}-RQ-\dot{n}Q$, where $\dot{n}$ is the He ii ionization rate per helium atom and R(t) is the average recombination rate which includes an effective clumping factor. The last term describes absorption by LLSs clumps. In Figure 6 we show the evolution of the volume filling factor assuming only recombination with moderate clumping factors. We also show the evolution including LLSs absorption as a further photon sink as in Bolton et al. (2009). The Q evolution adopting the Haardt & Madau (2012) emissivity is also shown for comparison. It appears that the He ii reionization can be completed down to z ∼ 3 assuming modest clumping especially when LLSs allow each AGN to ionize only a limited volume. The early hydrogen reionization is shown for comparison only in the recombination case since the mfp of hydrogen LLSs is larger than typical AGN separation and the photon sink term can be neglected. Of course, the knowledge of the detailed He ii reionization history requires three-dimensional hydrodynamical models. The adopted simplified description is only suggesting that some degree of inhomogeneity could help to reduce the tension between the possible presence of a high level of ionizing emissivity by AGNs at z ∼ 6 with an He ii reionization ending at z ∼ 3.

Zoom In Zoom Out Reset image size

A faint AGN population is emerging from recent multiwavelength surveys, in particular from those including X-ray detections. These are changing our perception about the AGN contribution to the ionizing UVB, although the selection of faint AGNs at the highest redshifts z ∼ 5–6 becomes difficult with the present instrumentation. Given the current observational limits, insight into the origin of ionizing flux at high redshifts can be gained from theoretical modeling. We have adopted our SAM R-SAM which successfully links the evolution of the galaxy population to that of the AGN population through the computation of the growth of SMBHs in the nuclei of galaxies. The model was able to reproduce several observables such as luminosity and mass functions of both galaxies and AGNs.

In this context we first concentrated on the comparison between the predicted and observed LFs at the UV rest-frame wavelength 1450 Å. The agreement is satisfactory at low and intermediate redshifts up to z ∼ 3. At higher redshifts the predicted LFs tend to overestimate the observed data by a factor of 2–3 at the fainter magnitudes M₁₄₅₀ ∼ −22. We have discussed the possible reasons for the discrepancy. From the observational point of view it is clear that surveys of optically/UV-selected AGNs which also rely on X-ray detection as a discriminant in the candidate selection process appear more complete on approaching the predicted values by the model. Indeed, without the X-ray discriminant, the point-like morphological criterion is needed to avoid confusion with the more numerous Lyman break galaxy population. But this introduces significant incompleteness at faint magnitudes. Moreover, several uncertainties in the selection process like, e.g., the lack of spectroscopic information for the faintest AGN candidates, make the estimates of the faint end of the AGN LF appreciably uncertain at very high redshifts.

Since faint AGN surveys could be affected by significant bias we have relied on our model to obtain an estimate of the hydrogen photoionization rate produced by the AGNs. To this aim, we have coupled our predicted LFs with the physical description of the escape fraction based on the blast-wave model already introduced in our R-SAM to describe the expansion of the blast wave produced in the interstellar medium of the hosting galaxy by the AGN-driven outflows. The blast wave clears the way for the UV photons emitted by the nuclear AGN which are free to escape outside the galaxy and ionize the IGM. To evaluate whether AGNs can ionize the IGM, we have first computed the AGN emissivity at the Lyman limit as a function of luminosity and redshift. Interestingly, the predicted AGN emissivity has a definite maximum at about M₁₄₅₀ ∼ −23 which is almost independent of redshift. Fainter sources give a negligible contribution because of the flattening of the shape of the LF and because the escape fraction of UV photons decreases on average in fainter AGNs and reaches 〈f〉 ≈ 0.1–0.3 at M₁₄₅₀ ∼ −22 at intermediate and high redshifts, respectively. Note that for AGNs of a given luminosity the escape fraction increases with redshift following the increase of the merging rate which is the main triggering of the QSO activity in our model. Since at z  >  3 AGNs with M₁₄₅₀ ∼ −22 show a duty cycle of about 10% and so represent about 10% of the host galaxy population (Fiore et al. 2012), then the bulk of the Lyman continuum emitters hiding faint AGN activity represents about 1% of the UV-dropout galaxy population similar to recent statistics on large Lyman break galaxy samples (Vanzella et al. 2010b). The resulting photoionization rate evolves strongly from z = 0 to z = 2, and then remains nearly constant to a value Γ₋₁₂ ∼ 1 up to z ∼ 4.5. A gradual decline is apparent at higher redshifts, consistent with the trend derived from most recent analyses of the Lyman forest in QSO spectra. We have also shown that a high level of AGN ionizing emissivity up to z ∼ 6 could still be consistent with a more gradual He ii reionization ending at z ∼ 3, especially if absorption by He ii Lyman limit IGM clouds is taken into account.

In summary, the main prediction of the model is that in a scenario where the AGN activity is linked to the dynamical galaxy evolution triggered by interactions and merging, only a mild high-redshift evolution of the hydrogen photoionization rate is expected in the redshift interval z ∼ 2–5. New deep multiwavelength surveys with accurate selection procedures are needed to build reliable AGN volume densities down to M₁₄₅₀ ∼ −23 in the relevant redshift interval. These values are critical to understand whether AGNs can play a key role in reionizing the high-redshift IGM. Our model suggests that it is time to reconsider the AGNs as the main population driving the ionization history of the universe over the whole known redshift interval (up to z ∼ 6) where the cosmic evolution of structures is observed.

We thank the referee for useful comments which significantly improved our paper. We are grateful to A. Cavaliere for a critical reading of the manuscript. We also thank F. Haardt for useful discussion on the helium reionization history.