Constraining the contribution of active galactic nuclei to reionisation

TL;DR

This study evaluates whether active galactic nuclei (AGN) could have significantly contributed to cosmic reionisation, concluding that AGN alone are unlikely to be the dominant source based on modeling and observational constraints.

Contribution

The paper introduces a semi-numerical model incorporating a quasar halo occupation distribution to assess AGN's role in reionisation, highlighting limitations of AGN-only scenarios.

Findings

AGN cannot match observed optical depth and neutral fraction simultaneously.

AGN-dominated models produce larger ionised bubbles and more 21cm power.

A mixed galaxy-AGN model is marginally consistent with current observations.

Abstract

Recent results have suggested that active galactic nuclei (AGN) could provide enough photons to reionise the Universe. We assess the viability of this scenario using a semi-numerical framework for modeling reionisation, to which we add a quasar contribution by constructing a Quasar Halo Occupation Distribution (QHOD) based on Giallongo et al. observations. Assuming a constant QHOD, we find that an AGN-only model cannot simultaneously match observations of the optical depth $\tau_e$, neutral fraction, and ionising emissivity. Such a model predicts $\tau_e$ too low by $\sim 2\sigma$ relative to Planck constraints, and reionises the Universe at $z\lesssim 5$. Arbitrarily increasing the AGN emissivity to match these results yields a strong mismatch with the observed ionising emissivity at $z\sim 5$. If we instead assume a redshift-independent AGN luminosity function yielding an emissivity evolution like that assumed in Madau & Haardt model, then we can match $\tau_e$ albeit with late reionisation, however such evolution is inconsistent with observations at $z\sim 4-6$ and poorly motivated physically. These results arise because AGN are more biased towards massive halos than typical reionising galaxies, resulting in stronger clustering and later formation times. AGN-dominated models produce larger ionising bubbles that are reflected in $\sim\times 2$ more 21cm power on all scales. A model with equal parts galaxies and AGN contribution is still (barely) consistent with observations, but could be distinguished using next-generation 21cm experiments HERA and SKA-low. We conclude that, even with recent claims of more faint AGN than previously thought, AGN are highly unlikely to dominate the ionising photon budget for reionisation.

Figures

Tables

Table 1. Table 1: Summary of models considered in section  3 to compare between the AGN and star-forming galaxies impacts on different reionisation constraints. Columns (from left to right) are: models’ names, the photon escape fractions from star-forming galaxies f esc , ⋆ subscript 𝑓 esc ⋆ f_{\mathrm{esc,\star}} and AGN f esc , AGN subscript 𝑓 esc AGN f_{\mathrm{esc,AGN}} , the ionisation rate used in Equation ( 1 ), the optical depth τ e subscript 𝜏 𝑒 \tau_{e} , and reionisation redshift z reion subscript 𝑧 reion z_{\rm reion} defined at neutral fraction limit x HI < 10 − 3 subscript 𝑥 HI superscript 10 3 x_{\rm HI}<10^{-3} .

Model	$f_{\mathrm{esc},\star }$	$f_{\mathrm{esc},\mathrm{AGN}}$	ionisation rate	${\tau }_e$	$z_{\mathrm{reion}}$
Galaxies	0.25	0.0	$f_{\mathrm{esc},\star }{R}_{\mathrm{ion},\star }$	0.057	7.5
constant QHOD AGN	0.0	1.0	$f_{\mathrm{esc},\mathrm{AGN}}{R}_{\mathrm{ion},\mathrm{AGN}}$	0.036	5.0
50-50	0.125	0.5	$f_{\mathrm{esc},\star }{R}_{\mathrm{ion},\star }+f_{\mathrm{esc},\mathrm{AGN}}{R}_{\mathrm{ion},\mathrm{AGN}}$	0.049	6.5
constant LF AGN	0.0	1.0	$f_{\mathrm{esc},\mathrm{AGN}}{R}_{\mathrm{ion},\mathrm{AGN}}$	0.048	4.0

Table 2. Table 2: Summary of parameters used in 21cmSense package to obtain the thermal noise sensitivity for each experiment. Columns (from left to right) are: experiments’ names, designs, antenna diameter, total collecting area, and the reciver temperature. The sky temperature is given by T sky = 60 λ 2.55 superscript 𝜆 2.55 \lambda^{2.55} K.

Experiment	Design	Diameter [m]	Collecting area [m2]	Receiver temperature [mK]
LOFAR	48 tiles of bow-tie high band antennae	30.75	35,762	140,000
HERA	331 hexagonally packed antennae	14	50,953	100,000
SKA	866 compact core antennae	35	833,189	100 Tsky + 40,000

Equations

$$\int f_{\mathrm{esc}}{R}_{\mathrm{ion}}dt\ge \int x_{\mathrm{HII}}{R}_{\mathrm{rec}}dt+(1-x_{\mathrm{HII}})N_{\mathrm{H}},$$

$$\frac{R_{\mathrm{rec}}}{V}=A_{\mathrm{rec}}(1+z{)}^{D_{\mathrm{rec}}}{\left[\frac{{\left(\Delta /B_{\mathrm{rec}}\right)}^{C_{\mathrm{rec}}}}{1+{\left(\Delta /B_{\mathrm{rec}}\right)}^{C_{\mathrm{rec}}}}\right]}^4,$$

$$\frac{R_{\mathrm{ion},\star }}{M_{\mathrm{h}}}=A_{\mathrm{ion}}(1+z{)}^{D_{\mathrm{ion}}}(M_{\mathrm{h}}/B_{\mathrm{ion}}{)}^{C_{\mathrm{ion}}}\exp \left(-(B_{\mathrm{ion}}/M_{\mathrm{h}}{)}^{3.0}\right),$$

$$\frac{M_{\mathrm{bh}}}{M_{\odot }}=A{\left(\frac{v_{\mathrm{cir}}}{159.4\mathrm{km}{\mathrm{s}}^{-1}}\right)}^5,$$

$$\frac{v_{\mathrm{cir}}}{\mathrm{km}{\mathrm{s}}^{-1}}=0.014{\left(M_h\sqrt{{\Omega }_m(1+z{)}^3+{\Omega }_{\Lambda }}\right)}^{1/3}.$$

$$\frac{L_B}{L_{\odot ,B}}=5.7\times 10^3\frac{M_{\mathrm{bh}}}{M_{\odot }}.$$

$$L_{\nu }=L_{912}{\left(\frac{\nu }{{\nu }_{912}}\right)}^{-1.57},$$

$$\frac{L_{912}}{\mathrm{ergs}{\mathrm{s}}^{-1}\mathrm{H}{\mathrm{z}}^{-1}}=10^{18.05}\frac{L_B}{L_{\odot ,B}}.$$

$$R_{\mathrm{ion},\mathrm{AGN}}={\int }_{{\nu }_{912}}^{\infty }\frac{L_{\nu }}{h\nu }d\nu =\frac{L_{912}}{1.57h}.$$

$$\varphi =\frac{{\varphi }^{\ast }}{10^{0.4(M_{\mathrm{break}}-M)(\beta -1.0)}+10^{0.4(M_{\mathrm{break}}-M)(\gamma -1.0)}},$$

$$\mathrm{N}={\left(\frac{M_h}{2.19\times 10^{12}}\right)}^{0.9}+0.023.$$

$$\delta T_b(\nu )=23x_{\mathrm{HI}}\Delta \left(\frac{{\Omega }_b{h}^2}{0.02}\right)\sqrt{\frac{1+z}{10}\frac{0.15}{{\Omega }_m{h}^2}}\left(\frac{H}{H+dv/dr}\right)\mathrm{mK},$$

Full text

††pubyear: 2017††pagerange: Constraining the contribution of active galactic nuclei to reionisation–References

1 Introduction

The nature of the sources driving the epoch of reionisation (EoR) in the early Universe remains uncertain. It is canonically believed that star-forming galaxies have provided the bulk of the ionising photon budget required to complete reionisation (Barkana & Loeb, 2001; Loeb & Barkana, 2001). This is because there is a significant decrease of observed active galactic nuclei (AGN) candidates at redshifts $z>3$, such that the contribution from star-forming galaxies is expected to well exceed that of AGN at $z>6$ (Shapiro & Giroux, 1987; Shankar & Mathur, 2007; Hopkins et al., 2007; Glikman et al., 2011; Masters et al., 2012; Haardt & Madau, 2012; Micheva et al., 2017; Ricci et al., 2017). However, there remain large uncertainties in the contribution of both star-forming galaxies and AGN to reionisation. Current constraints are now consistent with a minimal contribution from very low metallicity Population III stars (e.g. Robertson et al., 2015), but there is still the issue of the highly uncertain ionising photon escape fraction $f_{\mathrm{esc,\star}}$. Direct observations of $f_{\mathrm{esc,\star}}$ are quite difficult at $z\ga 4$ owing to the ubiquity of strong absorption systems that suppress Lyman continuum flux and the difficulty in removing foreground interlopers, but careful measurements generally indicate $f_{\mathrm{esc,\star}}$ less than a few percent (e.g. Grazian et al., 2016; Vasei et al., 2016), with some evidence for a higher $f_{\mathrm{esc,\star}}$ in lower-mass galaxies (Vanzella et al., 2016; Grazian et al., 2017; Bian et al., 2017).

Theoretical models have tried to constrain $f_{\mathrm{esc,\star}}$ indirectly by matching models to other data, making a variety of assumptions for $f_{\mathrm{esc,\star}}$ such as a constant (e.g. Robertson et al. 2013; Finkelstein et al. 2015; Ma et al. 2015; Hassan et al. 2016), redshift-dependent (e.g. Kuhlen & Faucher-Giguére 2012; Mitra et al. 2013; Finlator et al. 2015; Qin et al. 2017), mass-dependent (e.g. Gnedin 2007; Yajima et al. 2011; Wise et al. 2014; Paardekooper et al. 2015), and recently UV magnitude-dependent $f_{\mathrm{esc,\star}}$ (Anderson et al., 2017; Japelj et al., 2017) in order to match simultaneously various reionisation constraints. The currently-favoured lower value of Thomson scattering integrated optical depth ($\tau=0.058\pm 0.012$) measured by Planck (2016) prefers rather sudden and late reionisation scenarios, which relaxes the previously stringent constraints on the ionising photon budget. In Hassan et al. (2017), we performed a detailed Monte Carlo Markov Chain (MCMC) analysis to constrain our semi-numerical model to several EoR key observables and found that $f_{\mathrm{esc,\star}}$ is highly degenerate with the ionising emissivity amplitude, leading to a best fit value of $f_{\mathrm{esc,\star}}=0.25^{+0.26}_{-0.13}$, which allows a substantial range but is generally higher than available (lower-redshift) observations. Without a firmer understanding or direct measurement of $f_{\mathrm{esc,\star}}$, it is difficult to conclusively argue that star-forming galaxies can provide all the photons required for reionisation.

Recently, there has been renewed interest in assessing the contribution of AGN to the reionising photon budget. Previous estimates of the AGN contribution relied on an extrapolation to faint luminosities based on lower-redshift results. But recent deep observations have enabled a more direct characterisation of the faint end. Giallongo et al. (2015, hereafter G15) identified 22 faint AGN candidates at $z>4$ and inferred a significantly steeper faint-end slope than what is seen at lower redshifts. We note that claims of such a steep faint end remain controversial; for instance Parsa et al. (2017) was unable to confirm a substantial fraction of the G15 candidates based on additional multi-wavelength data. Furthermore, recent spectroscopic surveys (Kim et al., 2015; Jiang et al., 2016) have concluded that the observed quasars population at high redshift might not be enough to fully reionise the Universe. Nonetheless, the differing claims have led to speculation that AGN could provide the primary ionising photon contribution in order to keep the inter-galactic medium (IGM) highly ionised (e.g. Madau & Haardt, 2015). These claims further favor a late reionisation scenario in which the flatness observed in the ionising emissivity measurements by Becker & Bolton (2013) might arise naturally. In addition, they might also support the early and extended Helium reionisation observed by Worseck et al. (2016). Independently, Chardin et al. (2017) argued that the large scale opacity fluctuations in the Ly$\alpha$ forest measured by Becker et al. (2015) could be explained if AGN dominate the ionising UV background at $z\sim 6$ (see also Chardin et al. 2015). Hence the contribution of AGN to reionisation remains uncertain and potentially important or even dominant.

The idea that AGN might have driven cosmic reionisation has so far been investigated mostly in terms of global quantities, such as the ability to match the optical depth or comoving ionising emissivity constraints (e.g. Madau & Haardt, 2015; Mitra et al., 2016; Mao & Kim, 2016; Khaire et al., 2016; Qin et al., 2017). It remains to be demonstrated whether AGN-driven models are able to simultaneously satisfy all the current reionisation-epoch constraints. An important upcoming addition to the pantheon of constraints will be the 21cm EoR power spectrum, which may be substantially different for AGN- versus star formation-driven reionisation, if AGN and star-forming galaxies cluster in different ways as one might naively expect. An early attempt by Geil & Wyithe (2009) to assess the effect of AGN on the 21cm power spectrum using a seminumerical scheme concluded that the effect is likely to be small, but more recent semi-numerical models by Kulkarni et al. (2017) have suggested the opposite, that AGN produce signficantly different 21cm signal. However, Kulkarni et al. (2017) populates AGN only in the most massive halos using abundance matching to the halo velocity, employing the observed velocity-black hole mass relation at lower redshifts (Tremaine et al., 2002; Ferrarese, 2002), which thus effectively adopts a unity duty cycle of AGN for massive halos. However, recent results from Hyper Suprime-Cam suggest that quasars do not necessarily live in the most overdense regions where massive halos are expected to reside, and that their duty cycle is below a few percent (He et al., 2017). Accounting for sub-unity duty cycles inevitably drives black holes into lower-mass halos, altering the implied emissivity associated with haloes and epochs where they are not directly-measured. Moreover, the G15 data suggest that the AGN driving reionisation are rather faint, which may not be associated with the most massive halos. Without a proper treatment for AGN occupancy (duty cycle) and a more comprehensive analysis of all the implications of AGN-driven reionisation, it is difficult to properly assess the viability of this scenario.

In this paper, we build on our semi-numerical framework based on the SimFast21 code to evolve the EoR ionisation field, which allows us to examine a range of EoR observations as we have done in Hassan et al. (2016, 2017). To explore the AGN contribution, we populate AGN into halos with a more physically motivated approach that utilises both the observed luminosity function and abundance matching, thereby generating a Quasar Halo Occupancy Distribution (QHOD); our scheme partially follows the recipe summarized in Choudhury & Ferrara (2005). We constrain our QHOD to match the G15 AGN LF fit at $z=5.75$, and assign AGN randomly into halos. This QHOD predicts a duty cycle that is close to unity for extremely massive halos, but drops to sub-percent values at intermediate halo masses. To obtain the AGN emissivity, we utilise the strong correlation observed between the circular velocity and black hole mass following low redshift observations (Tremaine et al., 2002; Ferrarese, 2002). We account for this additional AGN photon contribution while we evolve our SimFast21 density and ionisation field, including the effects of recombination and time-evolving neutral fractions.

This work improves on previous efforts in several ways. First, using our SimFast21-based framework, we examine a wider variety of simultaneous constraints on the evolution of AGN-driven reionisation, including the Thomson optical depth, the mean cosmic neutral fraction evolution, and the ionising emissivity at the end of reionisation. Second, our model for populating quasars into halos is more realistic than previous works because we apply constraints beyond just abundance matching, allowing us to directly constrain the duty cycle of AGN as a function of halo mass. Third, we forecast upcoming 21cm EoR power spectrum measurements from LOFAR, HERA, and SKA, and illustrate how such future data might be able to constrain the fractional contribution of AGN to reionisation. Our primary conclusion is that it is very difficult to reconcile purely AGN-driven reionisation based on the (optimistic) G15 AGN luminosity function measurements with current global reionisation constraints. Future 21cm data should provide a new avenue to more precisely characterise the contribution of AGN to reionisation.

This paper is organized as follows: In Section 2, we describe our semi-numerical simulation and the AGN model implementation and calibration. We compare AGN with star-forming galaxies models in terms of their EoR observables, present the 21cm predictions, and discuss how future experiments can discriminate between these models in Section 3. We finally conclude in Section 4. Throughout this work, we adopt a $\Lambda$CDM cosmology in which $\Omega_{\rm M}=0.3$, $\Omega_{\rm\Lambda}=0.7$, $h\equiv H_{0}/(100\,\mathrm{km/s/Mpc})=0.7$, a primordial power spectrum index $n=0.96$, an amplitude of the mass fluctuations scaled to $\sigma_{8}=0.8$, and $\Omega_{b}=0.045$. We quote all results in comoving units, unless otherwise stated.

2 Simulations Using SimFast21

We use the recently developed Time-integrated version of our semi-numerical code SimFast21 (Santos et al. 2010) that has been presented in  Hassan et al. (2017). We here briefly review the simulation and defer to  Santos et al. (2010) for full details about the basic algorithm, and to  Hassan et al. (2016, 2017) for more information about our subsequent improvements.

The dark matter density field is generated using a Monte-Carlo Gaussian approach, which is dynamically evolved into the non-linear regime via applying the Zel’Dovich (1970) approximation. The dark matter halos are generated using the well known excursion set formalism (ESF). In the Time-integrated model, the ionised regions are identified using a similar form of the ESF that is based on comparing the time-integrated ionisation rate ${R_{\rm ion}}$ with that of the recombination rate ${R_{\rm rec}}$ and the local neutral Hydrogen density within each spherical volume specified by the ESF. Regions are flagged as ionised if:

[TABLE]

where $f_{\mathrm{esc}}$ is the photon escape fraction, $x_{\rm HII}$ is the ionised fraction, and $N_{\rm H}$ is the total number of hydrogen atoms. This is a generalized form of the ionisation condition in the Time-integrated model, which can be used for any ionising source or sink populations to run the reionisation calculations. With this ionisation condition, reionisation occurs more suddenly compared to our previous Instantaneous model developed in Hassan et al. (2016), in which the ionisation condition was based on an instantaneous comparison of ${R_{\rm ion}}$ and ${R_{\rm rec}}$. The more sudden reionisation is favoured by recent Planck (2016) data, and in Hassan et al. (2017) we showed that the Time-integrated ionisation condition produces larger ionised bubbles, resulting in 21 power spectrum enhancement on large scales.

2.1 Sink model

Reionisation, in short, is an evolving battle between ionising photon sources and sinks. To model sinks, we must account for the clumping effects from small scales below what we can directly evolve using the large-scale SimFast21 code (typically, sub-Mpc scales). We thus parametrize the inhomogeneous recombination rate ${R_{\rm rec}}$ from high-resolution full radiative transfer hydrodynamic simulations (hereafter 6/256-RT) (Finlator et al., 2015) as a function of over-density $\Delta$ and redshift $z$, as follows:

[TABLE]

where $A_{\rm rec}=9.85\times 10^{-24}$cm*-3s-1* (proper units), $B_{\rm rec}=1.76$, $C_{\rm rec}=0.82$, $D_{\rm rec}=5.07$. Consistent with  Sobacchi & Mesinger (2014), our recombination rate ${R_{\rm rec}}$ parametrization suppresses the ionisation and 21cm power spectrum on large scales. Full details about the inhomogeneous recombinations ${R_{\rm rec}}$ parametrizations and impact on the EoR observables can be found in Hassan et al. (2016).

We note that AGN-only reionisation scenarios are found to substantially heat the IGM (D’Aloisio et al., 2016; Oñorbe et al., 2017), which lowers the recombination rate. This may reduce ${R_{\rm rec}}$ by up to a factor of $\sim$ 2 in our AGN-only models, which in turn may slightly advance reionisation by AGN, and hence improving the viability of AGN-only models. We do not account for this effect in our calculation since we expect it to be sub-dominant compared to other effects related to halo growth, and here simply use the same sink model to compare reionisation histories produced by Galaxies versus AGN.

2.2 Source model: Star-forming galaxies

For the stellar contribution, we use a parametrization obtained from combining the 6/256-RT with larger-volume hydrodynamic galaxy formation simulation (Davé et al., 2013) (hereafter 32/512), that have both been shown to match a range of observations including lower redshift data. From these simulations, we parametrize the non-linear ionisation ${R_{\rm ion,\star}}$ rate as a function of halo mass $M_{\rm h}$ and redshift $z$ as follows:

[TABLE]

where $A_{\rm ion}=1.08\times 10^{40}\;{\rm M}_{\odot}^{-1}$s*-1*, $B_{\rm ion}=9.51\times 10^{7}\;{\rm M}_{\odot}$, $C_{\rm ion}=0.41$ and $D_{\rm ion}=2.28$. This ionisation rate is computed directly from the star formation rate (SFR) of these simulations based on stellar population models applied to star formation histories of simulated galaxies.

In Hassan et al. (2017) we considered a more generalized form of this source model, and found that constraining these parameters against several EoR observations using MCMC analysis resulted in best-fit values that matched the above parameters to within uncertainties, thereby validating the extrapolation from the small scales of 6/256-RT and 32/512 simulations to large scales covered by SimFast21 simulations. We further showed that using this non-linear ionisation rate relation boosts the small scales 21cm power spectrum as compared with models assuming a linear relation between the ionisation rate and halo mass; see Hassan et al. (2016, 2017) for more details.

2.3 Source model: AGN

The new aspect of the source model for this work is the AGN ionising photon output. We compute the ionisation rate from AGN ${R_{\rm ion,AGN}}$ following partially the recipe summarized in Choudhury & Ferrara (2005). Motivated by low redshift observations (Tremaine et al., 2002; Ferrarese, 2002), the basic assumption is that the black hole mass $M_{\rm bh}$ is strongly correlated with the hosting halo’s circular velocity $v_{\rm cir}$. We assume that this correlation is independent of redshift and valid during the reionisation redshifts. This correlation can be written as:

[TABLE]

where $A$ may be regarded as the black hole formation efficiency, which is our only free parameter in the AGN source model at fixed $f_{\mathrm{esc,AGN}}$. We then fix $A$ to match the AGN ionising emissivity constraints from G15, as we describe in §2.4. It is worthwhile to mention that this observed correlation (Equation 4) would naturally arise if one applies a self-regulation condition on the black hole growth as previously shown by Wyithe & Loeb (2003).

For our adopted cosmology, the circular velocity $v_{\rm cir}$ of a given halo mass $M_{h}$ is given by:

[TABLE]

Having obtained the black hole mass $\rm M_{\rm bh}$, the Eddington luminosity in the B-band is given by (Choudhury & Ferrara, 2005)

[TABLE]

Given this B-band luminosity, we must now determine the ionising photon output. Following Schirber & Bullock (2003) and Telfer et al. (2002), we assume the spectral energy distribution for AGN takes a power law form:

[TABLE]

where $\rm L_{912}$ is the luminosity at the Lyman limit that is given by:

[TABLE]

We then integrate the above over all frequencies to find the ionisation rate ${R_{\rm ion,AGN}}$ as follows:

[TABLE]

This explains how we compute the AGN ionisation rate ${R_{\rm ion,AGN}}$ given the host halo properties.

Next, we must populate the AGN into our halos. Here is where we make use of the G15 AGN Luminosity Function (LF). G15 evaluated this at $\lambda=1450{\rm\r{A}}$ for several redshifts higher than $z=4$. We then use G15 LF fit at their highest redshift z=5.75 to compute the number of AGN as a function of halo mass in our simulations. G15 LF at z=5.75 can be best fit using a double power law as follows:

[TABLE]

where $\phi$ is the comoving AGN density, $M$ is the absolute magnitude computed via the standard relation $M_{\rm AB}=-2.5\,{\rm log}_{10}(L_{\nu})+51.60$, log10 $\phi^{*}=-5.8$ Mpc*-3*, $M_{\rm break}=-23.4$, $\beta=1.66$, and $\gamma=3.35$.

Putting this together, our procedure to populate AGN into halos is as follows:

We bin the halo catalogs as a function of halo mass, and find the average halo mass in every bin of $\Delta\log_{10}M_{\rm h}=0.34$. 2. 2.

Using Equations 4 to 7, we compute the corresponding AGN $L_{1450}$ and $M_{1450}$ of each halo mass bin. 3. 3.

We then obtain the number of AGN for each halo mass bin using Equation (10), which turns out to always be less than the actual number of halos; the ratio of these numbers is the duty cycle of AGN for that halo mass bin. 4. 4.

We randomly assign the appropriate number of AGN into halos within that mass bin.

Note that in step (iii) one ideally may assume Poisson fluctuations around the number of AGN following McQuinn et al. (2009) QSO Method I, since the luminosity function, in principle, yields the average number of AGN at a given magnitude bin per the simulation volume. We ignore these fluctuations for two complementary reasons. First, the number of AGN obtained from the G15 LF is very large, particularly, at the faint end (N $\sim$ 106) around which the Poisson fluctuations can be neglected ($\sqrt{N}\sim$ 103). Second, as will be seen later, our results are mainly driven by the strong AGN clustering at the faint end, and hence adding Poisson fluctuations at the very bright end (e.g. first few magnitude bins) is unlikely to affect the results since bright sources are rare. In such a situation, the average number of AGN is a very good approximation to the actual number of AGN. We then round off the resulting AGN number in order not to populate halos with fractional AGN, but this in fact is a small correction.

Following the above procedure, we now have plausible AGN population in our simulation box at $z=5.75$. To quantify the evolution of the AGN population, we must make a choice regarding the evolution of AGN relative to that of the halos. Given that theoretical predictions for AGN evolution are relatively uncertain, the simplest assumption is to assume that the relationship between AGN and their host halos do not change at $z\geq 5.75$. In other words, we assume that the Quasar Halo Occupation Distribution (QHOD) is non-evolving. This is a reasonable assumption since the HOD of galaxies have been studied extensively (see Yoshikawa et al. 2001; Berlind et al. 2003, and references therein) and it is found that the HOD is nearly a redshift-independent quantity.

The QHOD as a function of halo mass $M_{h}$ can directly be calculated from G15 LF fit at z=5.75 (Equation 10) as the ratio between the number of AGN to that of their hosting halos for each halo mass bin. Our QHOD can be well fit with a constant plus a power law as follows:

[TABLE]

Note that the QHOD changes with different values of our free parameter $A$ relating black hole mass to circular velocity, which translates into a shift in magnitudes of the AGN. Here we have used $A=5\times 10^{5}$, a value at which our constant QHOD AGN model is calibrated to reproduce the G15 ionising emissivity constraints, as will be discussed next in §2.4.

Figure 1 shows the QHOD computed from G15 observations at z=5.75 (closed circles) and our QHOD fit in Equation (11). The QHOD represents a plausible description of the AGN occupancy in their hosting halos. Indeed, this can also be regarded as an AGN duty cycle, if one (reasonably) postulates that every halo contains a black hole but only some fraction of them are detectably active. The He et al. (2017) observations suggest a duty cycle of $0.001-0.06$ for moderate-mass halos, which is somewhat lower than our model assumes but qualitatively agrees with the trend that the duty cycle is smaller in lower-mass halos. For comparison, we also plot the QHOD at $z=4.25$ that is computed from G15 LF at that redshift bin (open circles in Figure 1). We notice that the QHOD data at $z=5.75$ and $z=4.25$ are fairly similar, differing by $\sim 30\%$ for all $M_{h}\ga 10^{10}M_{\odot}$. This suggests that the QHOD doesn’t evolve strongly with redshift, and motivates us to fiducially assume that the QHOD does not evolve. We will call this the “constant QHOD" case. In this case, we replace Equation (10) with Equation (11) in step (ii) to compute directly the number of AGN in each halo mass bin at higher redshifts. Note that in step (iv) AGN assignment is completely random and redshift-independent. As a result, halos with AGN may or may not have AGN at the next time step. This is realistic since the simulation time step (dz=0.125) is typically larger than the AGN lifetime.

As a counterpoint to this case, we also consider a model where the AGN luminosity function is constant with time. Here, we calculate the AGN number at all redshifts based on G15 LF fit at $z=5.75$ (Equation 10). We will call this the “constant LF" case. This is less realistic because the QHOD here increases strongly with redshift, since there are many fewer halos at higher redshifts but the number of AGN remains fixed following the assumed constant LF. Also, observations of the AGN LF at lower redshifts ($z\la 6$) exhibit significant evolution, so it seems unlikely that this evolution should suddenly cease. Nonetheless, the emissivity evolution in this case turns out to be similar to the AGN comoving ionising emissivity model assumed by Madau & Haardt (2015), hence it represents an interesting contrasting case that we will examine in §3.

2.4 AGN source model calibration

Our next task is to calibrate the relationship between black hole mass and circular velocity via the normalization parameter A in Equation (4). Observationally, this parameter is not tightly constrained and has only been measured at low redshifts.

We first calibrate the constant QHOD and constant LF AGN models to at least match the G15 ionising emissivity constraints at z=5.75 in order to verify the possibility to complete reionisation solely by AGN. This we achieve by tuning the black hole formation efficiency $A$ in our AGN models to match the total ionising emissivity measurements at 912 Å ($\epsilon_{912}$), which is the total escaped $L_{912}$ of all AGN divided by the simulation comoving volume. The simulation configurations of these models are presented in §3 with the rest of our fiducial models. We assume $f_{\mathrm{esc,AGN}}=100\%$ for AGN, which is standard (e.g. Madau & Haardt, 2015).

We find that the constant LF AGN model can match the G15 ionising emissivity constraints with $A=10^{6}$ whereas the constant QHOD AGN model requires $A=5\times 10^{5}$. We note that what is really constrained here is the product of $Af_{\mathrm{esc,AGN}}$, so we have the freedom to keep $A$ fixed and tune $f_{\mathrm{esc,AGN}}$ instead; all our results would be unchanged. In this case, the constant QHOD and constant LF AGN models would require $f_{\mathrm{esc,AGN}}=50,100\%$ at $A=10^{6}$ to match the G15 constraints, respectively. Note that we do make the assumption here that the product $Af_{\mathrm{esc,AGN}}$ does not vary with redshift.

In Figure 2, we show the comoving ionising emissivity evolution obtained with the procedure discussed above. The constant LF AGN model produces a slowly growing emissivity, which is similar to the evolution expected from models in which ionising radiation is dominated by star-forming galaxies, and similar in shape to that assumed in Madau & Haardt (2015). While matching the G15 constraints at $z=5.75$, the constant LF AGN model under-estimates the ionising emissivity by a factor of $\sim$ 3 as compared with G15 constraints at $z=4.75$ and $z=4.25$. In contrast, the emissivity from our fiducial constant QHOD AGN model matches simultaneously G15 data at several redshifts bins due to the rapidly growing emissivity evolution as expected from an AGN dominated model. This further validates our assumption that using constant QHOD is a more physically motivated approach than using the constant LF.

Note that we have intentionally not applied a magnitude cut-off in computing the integral of emissivity ($\epsilon_{912}$); G15 used a cutoff of $M_{1450}=-18$. At $z=5.75$, the total comoving ionising emissivity is $\epsilon_{912}=2.12\times 10^{24}$ erg s*-1* Mpc*-3* Hz*-1*, whereas with a magnitude cut-off of $M_{1450}=-18$ it becomes $\epsilon_{912}=1.58\times 10^{24}$ erg s*-1* Mpc*-3* Hz*-1*. This shows that those fainter AGN contribute $\sim$ 25$\%$ to the total emissivity, which is modest but not negligible. We include this in order to check whether including all faint AGN would allow reionisation completion to be consistent with neutral fraction and optical depth constraints. This means that we are effectively studying an optimal case for reionisation by AGN, since those fainter than $M_{1450}>-18$ might already be a part of the galaxy population as discussed in Chardin et al. (2017), due to the overlap between the galaxy and AGN luminosity functions at this faint limit. Applying a magnitude cut-off would suppress the ionising emissivity and further delay reionisation.

In summary, we have described our procedure to obtain the ionising emissivity of AGN as a function of halo mass, and then populate AGN into halos within SimFast21 via constraining the halo occupancy of AGN (QHOD) using the G15 AGN LF. The total emissivity is calibrated to match that observed by G15 at $z=5.75$, which fixes our free parameter relating black hole mass to halo circular velocity. Our fiducial model assumes a constant QHOD, and we will also consider a constant AGN LF. We now study the predictions of our model for reionisation observables, and compare the results with our previous SimFast21 models where we considered only star-forming galaxies.

3 EoR Observables

3.1 SimFast21 runs

We run all of our EoR realizations using the Time-integrated model (Hassan et al., 2017) to establish a proper comparison between the different source models. Using the same density field and halo catalogs generated in a box size L$=300$ Mpc and N$=560^{3}$ number of cells, we run 4 different EoR models based on different ionisation sources as follows (and summarized in Table 1):

• •

Galaxies: This model only considers ionising photons emitted by star-forming galaxies using Equation (3) with parameters: $f_{\mathrm{esc,\star}}$ = 0.25 , $A_{\rm ion}$ = 4.27$\times 10^{39}$, $C_{\rm ion}$ = 0.44. These parameters are suggested by our recent MCMC analysis in Hassan et al. (2017) to match simultaneously various EoR constraints including the SFR densities at several redshift bins as compiled by Bouwens et al. (2015), integrated ionising emissivity at z $\sim$ 5 by Becker & Bolton (2013) and Planck (2016) optical depth.

• •

constant QHOD: This is our fiducial AGN model in which the AGN are the only source for ionising radiation using $f_{\mathrm{esc,AGN}}$=1.0, and our QHOD fitting function (Equation 11) computed from G15 LF fit at $z=5.75$.

• •

50-50: This model contains an equal contribution from the Galaxies and constant QHOD models, specifically we use $f_{\mathrm{esc,\star}}$ = 0.125 and $f_{\mathrm{esc,AGN}}$ = 0.5.

• •

constant LF: This is our alternative AGN model which uses the actual LF fit of G15 at $z=5.75$ to compute the number of AGN at all redshifts, with $f_{\mathrm{esc,AGN}}$=1.0.

3.2 Ionising emissivity

We begin by comparing the integrated ionising emissivity $\dot{N}_{\rm ion}$ of these models, which is the total number of ionising photons per second per comoving volume. We compare our models with results from  Qin et al. (2017) based on the DRAGONS simulation, which uses the Meraxes semi-analytic galaxy formation model built upon the Tiamat $N$-body simulation. The Qin et al. (2017) model is able to track the growth of central super-massive black holes and reproduce wide range of observations including the observed quasar LF from $z\sim 0.6-6$. Their model predicts that AGN contribution to EoR is minimal, so it is interesting to compare our AGN emissivity with theirs.

Figure 3 shows the redshift evolution of the comoving integrated ionising emissivity $\dot{N}_{\rm ion}$, in units of $10^{51}$ s*-1Mpc-3*, for our four models. For the Qin et al. (2017) models (blue lines), we show their full galaxies+AGN model as well as their AGN-only contribution. At $z\la 5$, we show the observational constraints from Becker & Bolton (2013) inferred from Ly$\alpha$ forest measurements.

Our Galaxies model skirts the upper limit of Becker & Bolton (2013) constraints while simultaneously calibrated to reproduce the Bouwens et al. (2015) SFR and Planck (2016) optical depth, as discussed in Hassan et al. (2017). Adopting a mildly redshift-dependent or even mass-dependent $f_{\mathrm{esc,\star}}$ would permit a better match with the amplitude and flat redshift dependence of the Becker & Bolton (2013) emissivity measurements, as suggested by Mutch et al. (2016), without much altering the Thomson optical depth.

Our fiducial constant QHOD model shows a much more rapid growth of ionising emissivity with time than the Galaxies model, which matches the Becker & Bolton (2013) ionising emissivity at the upper end of the observed redshift range but overshoots the low end. In this model, the AGN contribution overtakes the Galaxies contribution at $z\sim 3$, which is in agreement with what is typically inferred (e.g. Haardt & Madau, 2012).

We further see that our 50-50 model (dotted black in Figure 3) is similar and much closer to the Galaxies model than to the constant QHOD model. This indicates that the contribution from star-forming galaxies dominates the ionising emissivity while AGN contribution is minor.

Finally, the constant LF model shows a relatively shallow evolution, approximately parallel to the galaxies case but significantly lower in amplitude. This falls just below the ionising emissivity data at $z\la 5$.

The blue dotted and dashed lines show the evolution of the emissivities from the Qin et al. (2017) model, with the latter showing only the contribution from AGN. Their AGN contribution shows a similar redshift dependence to our constant QHOD model, which further supports the validity of this assumption, at least down to $z\sim 6$. It is possible that our assumption breaks down at $z\la 5$ as their AGN contribution flattens, and if ours did this then the agreement with the observed evolution of the emissivity would improve. However, DRAGONS also predicts that galaxies dominate the ionising photon budget at all redshifts, which may be contrary to studies of the hardness of the ionising background in the $z\sim 2-4$ IGM (e.g. Schaye et al., 2007; Oppenheimer et al., 2009).

The emissivity evolution in our Galaxies model is similar to that from the Galaxies+AGN model by Qin et al. (2017), whereas there is a difference in the amplitude due to these models’ differences in the $f_{\mathrm{esc}}$ treatment (constant versus redshift-dependent). As previously noted, the constant LF AGN model shows a slowly growing emissivity similar to those of our Galaxies model and Galaxies + AGN model by Qin et al. (2017).

In summary, our fiducial constant QHOD AGN model matches the Becker & Bolton (2013) emissivity measurements reasonably well, at least at $z\sim 5$, but shows a dramatically different redshift evolution compared to our Galaxies and constant LF models. The Qin et al. (2017) DRAGONS model shows an evolution during reionisation that is consistent with our constant QHOD model for their AGN component, and with our Galaxies model for their overall emissivity (which is dominated by galaxies), but somewhat lower in amplitude in each case. All these models are broadly in agreement with the emissivity measures at $z\sim 5$ given current uncertainties.

3.3 Global ionisation history and optical depth

We now explore our model predictions for other current observational constraints on the global evolution of reionisation, particularly the evolution of the volume-weighted neutral fraction $x_{\rm HI}$ and the integrated Thomson optical depth to electron scattering to the CMB $\tau_{e}$.

Figure 4 shows a comparison between our models in terms of their global ionisation history, as characterised by the volume-weighted average neutral fraction $x_{\rm HI}$. We see that our Galaxies and 50-50 models are consistent with several Ly-$\alpha$ forest measurements (shaded areas from  Fan et al. 2006; Becker et al. 2015; Bouwens et al. 2015 and orange upper limits by McGreer et al. 2015). While both are consistent with data, the 50-50 model delays reionisation by $\Delta z\sim 1.0$, owing to the fact that galaxies are the main driver of reionisation as discussed in the previous section and their contribution has been halved in this model.

Turning to our AGN-driven reionisation models, we see that both the constant LF and constant QHOD reionise the Universe very late, with $\rm{x_{\rm HI}}<10^{-3}$ not occurring until $z\sim 4$ and $z\sim 5$, respectively. This is highly inconsistent with the Bouwens et al. (2015) constraints as seen in Figure 4, as well as direct observations of the Ly$\alpha$ forest in quasars at these epochs (Fan et al., 2006; Becker et al., 2015). The constant LF AGN model starts reionisation earlier than constant QHOD AGN model due to its higher emissivity at high redshifts (Figure 3). This too-late reionisation strongly suggests that AGN are unable to drive the bulk of reionisation, when constrained to match the observed ionising emissivity after the end of reionisation ($z\sim 5$).

As mentioned earlier, all models use the Time-integrated ionisation condition (Equation 1) to identify the ionised regions in the excursion set-formalism. We showed in Hassan et al. (2017) that this ionisation condition results in a more sudden reionisation as opposed to our previous Instantaneous model (Hassan et al., 2016) which yields an extended reionisation scenario. From Figure 4, all models yield a fairly sudden reionisation, consistent with our previous results, except the constant LF AGN model that shows an extended reionisation. This is because the fixed LF likely over-estimates the number of AGN at high redshifts, and furthermore the source population does not grow in concert with the growing sink population, which delays reionisation.

Figure 5 shows the evolution of the Thomson scattering optical depth ($\tau_{e}$) as a function of redshift in these models. The Galaxies model is consistent with the recent Planck (2016) measurements (red shaded areas), mostly because it was constrained to do so via MCMC. The 50-50 and constant LF yield a lower optical depth of about $\tau\sim$ 0.049, consistent with the lower limit of 1-$\sigma$ level. However, the optical depth obtained by our fiducial constant QHOD AGN model is very low ($\tau\sim 0.036$) at the lower limit of the 2-$\sigma$ level. Essentially, this model does not produce enough early photons in order to obtain a sufficient ionised path length to the CMB.

One might question why Madau & Haardt (2015) was able to match these constraints based on the G15 model and thereby argue for purely quasar-driven reionisation, whereas we reach an opposite conclusion. There are two main differences. First, Madau & Haardt (2015) made a rough fit to the ionising emissivity, which resulted in a much flatter redshift dependence than we obtain from our constant QHOD model, more like our constant LF model except higher in amplitude by $\sim\times 2-3$ (see Figure 2). Indeed, our constant LF model results in an ionising emissivity evolution much like theirs, and is consistent (albeit marginally) with $\tau_{e}$. Second, Madau & Haardt (2015) purely relied on counting emitted photons, and did not account for an evolving spatially-clustered nature of sinks which absorbs more photons, and hence retards reionisation particularly at early epochs compared to a spatially-homogeneous clumping factor model. These differences result in substantially more neutral gas at early epochs in our model, and lead to strong disagreement with the measured $\tau_{e}$ as well as $x_{\rm HI}$.

Our 50-50 model can plausibly match the constraints within current $1\sigma$ uncertainties. Recall that for our Galaxies model, Hassan et al. (2017) found $f_{\mathrm{esc}}=0.25^{+0.26}_{-0.13}$, which means that the galaxy contribution alone in the 50-50 model corresponds to an escape fraction that is at the $1\sigma$ bound allowed by the MCMC fit. Correspondingly, the predicted $\tau_{e}$ and the redshift of reionisation are near the $1\sigma$-low end of their respective allowed observational ranges.

The Qin et al. (2017) AGN model yields an optical depth of $\tau\sim 0.025$, corresponding to an end of reionisation at $z\sim 3$. Our constant QHOD AGN model would obtain similar results if a magnitude cut-off had been implemented in the ionising emissivity integral to exclude the very faint AGN. Even with our more favorable case of reionisation by AGN, reionisation still occurs very late.

Overall, we find that our constant QHOD model that is constrained to match the AGN emissivity of G15 can adequately match the global ionising emissivity at $z\sim 5$, but strongly fails to reionise the Universe by $z\sim 6$ and produces too low Thomson optical depth. A constant LF model does somewhat better at matching constraints, but the underlying assumption is not physically well-motivated, does not match the observed emissivity evolution from $z\sim 6\rightarrow 4$, and is inconsistent with the self-consistent calculations of Qin et al. (2017). A 50-50 model is still within the allowed bounds of the observations considered here, but any larger contribution from AGN would be disfavoured – and we reiterate that our AGN model is already pushed towards increasing the AGN emissivity as much as possible. Our $x_{\rm HI}$ and $\tau_{e}$ constraints for all these models are listed in Table 1. Our results thus suggest that AGN-dominated reionisation is highly unlikely, and therefore that galaxies dominate the ionising photon budget during the EoR.

3.4 EoR topology

We have shown that our AGN-only models are highly disfavoured given current observational data. However, a 50-50 model is still permissible, if only marginally. Clearly, increasing the precision of current measures should in principle enable more stringent constraints on the relative contribution of AGN and star-forming galaxies to reionisation. But we can also appeal to other aspects such as the topology of reionisation in order to discriminate between models. This will be particularly fruitful in the era of 21cm EoR experiments, which will quantify the power spectrum of neutral hydrogen on large-scales. In this section, we discuss the topology of neutral gas in our various models, and in the following section we will quantify this by forecasting the 21cm power spectrum for upcoming experiments.

We first investigate these models’ differences in terms of their ionisation field maps. We choose to compare these models’ maps at a fixed neutral fraction, since we have shown in Hassan et al. (2017) that the 21cm fluctuations are more sensitive to the topology of the ionisation field while the density field contribution is secondary. This then allows us to compare the topology even though the actual redshift where a given ionisation occurs varies substantially between models.

Figure 6 shows 2D maps of the ionisation field of all models at different neutral fractions (${x_{\rm HI}}\approx 0.25,0.5,0.75$), projected through the entire volume. Note that the redshifts at which this occurs are much later for the AGN-only models, consistent with their late reionisation.

Our Galaxies model shows a range of bubble sizes, with many small bubbles around low-mass galaxies that have low clustering. The largest bubbles towards the end of reionisation span $\sim 100$ Mpc, in agreement with many previous studies (e.g. Barkana & Loeb, 2004; Furlanetto et al., 2004; Mesinger et al., 2011; Zahn et al., 2011; Iliev et al., 2014; Majumdar et al., 2014). It is clear that there will be significant power on all scales owing to this topology.

The constant QHOD and constant LF AGN models show fairly similar H ii bubble sizes and distributions. Compared to the Galaxies case, there are fewer small bubbles owing to the fact that AGN tend to populate more massive halos than the typical galaxy contributing to reionisation in the Galaxies model. However, there are still some small halos hosting small bubbles even in the AGN case. This contrasts with the Kulkarni et al. (2017) AGN model where AGN are assigned only in the massive halos using a circular velocity cut-off, and thus their model does not yield any small-scale H ii regions (see their Figure 2). Nonetheless, because the duty cycle of AGN in our model is very low in low-mass halos, many fewer small bubbles are seen compared to the Galaxies case.

For large ionisations, the ionisation maps of the constant QHOD and constant LF AGN models display larger H ii regions as compared with those of Galaxies model. This is necessary to compensate for the lack of numerous small bubbles, in order to achieve a similar neutral fraction. This is driven by the strong AGN clustering as suggested by the input G15 LF at its faint end.

The 50-50 model is, not surprisingly, intermediate between the Galaxies and AGN-only models. The small bubbles are less prominent owing to the lower $f_{\mathrm{esc,\star}}$. The maximum bubble sizes are comparable to but slightly larger than in the Galaxies model. Hence the star-forming galaxies tend to drive the topology even when substantial AGN are present. This comes from the facts that each halo mass bin (magnitude bin) includes fewer AGN than the number of possible hosting halos as implied by the AGN duty cycle estimates (e.g. Shankar et al. 2009), and hence star-forming galaxies would overcome the impact of those fewer AGN on the ionisation field topology. The difference between AGN (constant QHOD and constant LF) and star-forming galaxies (Galaxies and perhaps 50-50) dominated models increases as reionisation proceeds and becomes clear at late stage of reionisation (bottom row for ${x_{\rm HI}\sim 25\%}$ of Figure 6). We expect these trends to be reflected quantitatively in their 21cm power spectra.

3.5 The 21cm power spectrum

Using the ionisation fields of these models at fixed neutral fractions (Figure 6), we now compute our key EoR observable, namely the 21cm power spectrum. Assuming that the spin temperature is much higher than the CMB temperature, the 21cm brightness temperature takes the following form:

[TABLE]

where $dv/dr$ is the comoving gradient of the line of sight component of the peculiar velocity. Using this equation, it is straightforward to create the 21cm brightness temperature boxes from which we compute the 21cm power spectrum as follows: $\mathrm{\Delta^{2}_{21}\equiv k^{3}/(2\pi^{2}\,V)<|\delta T_{b}(k)|^{2}_{k}>}$.

In top panels of Figure 7, we show a comparison between our models’ predicted 21cm power spectrum, as well as with the AGN-dominated models of Kulkarni et al. (2017). Bottom panels show the ratio of each model 21cm power spectrum to the Galaxies model, in order to clearly display the models differences as a function of the scale $k$. Consistent with the ionisation maps (Figure 6), the Galaxies model has less power than the AGN-only models, which at a fixed $x_{\rm HI}$ are themselves rather similar. This shows that the 21cm power spectrum is somewhat insensitive to the method by which we populate AGN (constant QHOD versus constant LF), even though there is a clear difference in the reionisation histories due to differences in the ionising emissivity evolutions (see Figures 3 &  4). Relative to the Galaxies case, at early epochs the differences with the AGN-only models peaks at intermediate scales (k $\sim$ 0.3 $-$ 0.1), typical of the AGN bubble sizes. At later epochs, however, there is not much scale dependence to the variation, and the AGN models are simply about a factor of $1.5-2\times$ higher than the Galaxies model, owing to the stronger clustering of G15 AGN observations.

The 50-50 models yield similar 21cm power spectra both on small and large scales to the Galaxies case at the early ($x_{\rm HI}\sim 0.75$) and intermediate (${x_{\rm HI}\sim 0.5}$) stages of reionisation. This confirms our previous finding that star-forming galaxies are dominant in determining the ionised regions properties and hence their associated 21cm fluctuations during the early stages of reionisation. At later stages, there is an increasing difference between the 50-50 model and the Galaxies case, owing to the increased contribution of AGN at later epochs. Hence it may be optimal to look at the latter stages of reionisation in order to obtain more quantitative constraints on the fractional contribution of AGN.

Our AGN models agree with the Kulkarni et al. (2017) AGN model reasonably well for the large-scale 21cm power spectrum, while their small scale power is suppressed by a factor of $\sim\times 2-2.5$ relative to ours. As mentioned before in §3.4, our AGN model yields small-scale ionised regions owing to populating AGN randomly at all halo mass bins (magnitude bins) into their hosting halos using the actual number of AGN as suggested by the G15 LF observations, as opposed to a unity duty cycle in massive halos in the Kulkarni et al. (2017) AGN model. By constraining our LF to that observed, our AGN model predicts a non-unity duty cycle, which occasionally populates small halos with AGN and thus boosts the 21cm power spectrum on small scales.

3.6 Forecasting 21cm power spectra to constrain AGN models

Our work has shown that AGN driven-EoR models are photon-starved, in agreement with many others, and as such plausible AGN models are unlikely to fully drive reionisation. Nonetheless, a substantial AGN contribution such as in our 50-50 model is still allowable given current data, which begs the question, will future 21cm data provide more stringent constraints on the contribution of AGN to reionisation?

To answer this question, we focus our analysis on three different low-frequency radio interferometer designed to measure the 21cm EoR power spectrum: the Low Frequency Array (LOFAR)111http://www.lofar.org/, the Hydrogen Epoch of Reionisation Array (HERA)222http://reionization.org, and the Square Kilometer Array (SKA-Low)333https://www.skatelescope.org.

To forecast the power spectra for these facilities, we use the same recipe presented in Hassan et al. (2017), outlined as follows: We first select the redshift (observed frequency) at which we compute the 21cm power spectrum for each model. To establish a proper comparison, we operate these three array designs, with parameters summarized in Table 2, in a drift-scanning mode for 6 observing hours per day for 180 days at 8 MHz bandwidth. We then add the total uncertainty which includes the thermal noise and sample variance using the 21cmSense444https://github.com/jpober/21cmSense, a package for calculating the sensitivity of 21cm experiments to the EoR power spectrum. We refer to Parsons et al. (2012) for basics of the radio interferometer sensitivities, to Pober et al. (2013, 2014) for more details on observation strategies and foreground removal models, and to Hassan et al. (2017) for the experiments designs and configurations. We only change the foreground removal method to use the optimistic model developed in Pober et al. (2014) in which the foreground wedge extends to the Full-Width Half-Max (FWHM) of the experiments’ primary beam. This will extend our analysis to cover more large scales than using a moderate model in which foreground wedge extends only to 0.1 $h$ Mpc*-1* beyond horizon limit.

From Figure 7, we notice that the large variations between the Galaxies and constant QHOD AGN models occurs at the late stages of reionisation (${x_{\rm HI}\sim 25\%}$) when the HII bubbles begin to overlap. We thus create our 21cm mock observations at these epochs. Since the reionisation occurs very late in the constant QHOD AGN model (see Figure 4), we re-tune the model (using $A=1.25\times 10^{7}$, which we note would substantially overproduce the ionising emissivity constraints) to match the optical depth as obtained by Galaxies model. We then perform our 21cm mock observations at $z=8$ and ${x_{\rm HI}}\sim 0.3$ for both models in order to conduct a comparison at the same redshift and neutral fraction.

Figure 8 shows a 21cm mock observation comparison at z=8 (${x_{\rm HI}}\sim 0.3$) between Galaxies and the re-tuned constant QHOD AGN models. The shaded area shows the total uncertainties (1-$\sigma$ level) expected from the experimental designs summarised in Table 2 using the 21cmSense package. The vertical lines mark the scale at which the 1-$\sigma$ errorbars from a specific experiment overlap, corresponding to the sensitivity limit for each experiment where the models can be distinguished.

Given that it can only probe relatively large scales, LOFAR will have some difficulty discriminating between the Galaxies and constant QHOD models, as they lie within $\sim$ 2-$\sigma$ of each other. However, HERA should be able to distinguish between these models rather well on scales above about 10 Mpc, and the larger baselines of SKA-low will enable discrimination to significantly smaller scales. From Figure 8, we see that LOFAR, HERA, and SKA can discriminate between these models during the latter stages of reionisation at scales of $k<0.21$ Mpc*-1* ($>30$ Mpc), $k<0.53$ Mpc*-1* ($>12$ Mpc) and $k<1.66$ Mpc*-1* ($>4$ Mpc), respectively.

The fact that HERA and SKA can easily discriminate AGN-only and Galaxies-only models suggests that it may be possible to constrain the fractional contribution of AGN using such data. We thus repeat the same steps above, except now for the 50-50 model, tuned to match Galaxies $\tau_{e}$.

Figure 9 shows a forecasting comparison between the Galaxies and 50-50 models. Since the 50-50 model yields 21cm power amplitude that is closer to Galaxies than constant QHOD model does, the scales at which experiments overlap are shifted towards large scales (compare vertical dashed lines in figure 8 versus 9). This shows that LOFAR is unlikely to discriminate between these models, unless a very optimistic foreground removal is applied to detect the signal on large scales ($k<0.12$ Mpc*-1* ($>53$ Mpc)), which are highly contaminated by foregrounds (Pober et al., 2014). Given a successful foreground removal, HERA and SKA can discriminate between the Galaxies and 50-50 models during the latter stages of reionisation at scales of $k<0.46$ Mpc*-1* ($>14$ Mpc) and $k<0.65$ Mpc*-1* ($>10$ Mpc), respectively as shown by vertical dashed lines in Figure 9.

Note that we have shown 1-$\sigma$ uncertainties, which is unlikely to be sufficient to robustly discriminate between Galaxies and AGN models. If one instead requires 3-$\sigma$ to distinguish between the models, then LOFAR fails to distinguish between our models, but HERA and SKA are still successful albeit on scales somewhat larger than those obtained with the 1-$\sigma$ limit.

Figure 8 & 9 both illustrate how future 21cm observations could potentially help constrain the nature of the source population. While current observations cannot rule out the 50-50 model, in principle HERA and SKA should be able to do so straightforwardly, assuming they can reach their target sensitivities. Until these facilities come online, ancillary observations will continue to improve. Hence comprehensive models that are able to make predictions for, and constrain to, a wide variety of EoR data are vital for optimizing the scientific information extracted from future observations including 21cm data.

4 Conclusion

We have presented predictions for the 21cm power spectrum arising from AGN-driven reionisation models, and contrasted them with predictions from galaxy-driven models and models with a mixture of sources. The AGN source population is placed into galaxy halos using a physically motivated prescription based on the G15 AGN observations, deriving a Quasar Halo Occupancy Distribution of AGN at $z=5.75$ from this and using it to evolve the number of AGN to higher redshifts. This framework is implemented into our Time-integrated version of SimFast21, which self-consistently accounts for recombinations and the evolution of structure. We have calibrated these AGN models to reproduce ionising emissivity constraints, and compared them with models in which ionising radiation is dominated by star-forming galaxies. Our key findings are as follows:

• •

When tuned to match the G15 ionising emissivity constraints, AGN-only models produce very late reionisation at $z\ll 6$ (Figure 4). If we assume a constant halo occupancy for AGN as is consistent with other observational constraints and models, then the predicted Thomson optical depth is only 0.036, well below Planck constraints (see Figure 5). This strongly disfavours AGN as providing the dominant source for reionising photons.

• •

We determine a quasar halo occupancy distribution (QHOD) that is near unity for very massive halos, but drops to sub-percent level for more typical halos. This is directly interpretable as a duty cycle for AGN. This also explains why reionisation is so late in this AGN-only model, because AGN populate massive halos more frequently and hence their emissivity contribution grows strongly only at late epochs (Figure 3), thereby not ionising enough volume to match the optical depth measurements.

• •

Our results are consistent with those from the AGN-only models by Qin et al. (2017) using the DRAGONS semi-analytic models, who also found that AGN could not dominate reionisation. Our Galaxies model is also in broad agreement with their full model, supporting their result that star-forming galaxies can provide sufficient photons as a function of redshift to reionise.

• •

A model where we assume a constant AGN luminosity function at $z\geq 5.75$ can barely match the Planck $\tau_{e}$, but still reionises very late, and moreover it is not physically well motivated and disagrees with the measured evolution of the AGN luminosity function at $z\sim 4-6$ (Giallongo et al., 2015). It does, however, result in a global emissivity evolution similar to that assumed in Madau & Haardt (2015), but even in this case we do not confirm their result that such a model is viable, likely because we include the effects of recombinations along with a more sophisticated accounting of the clustering of sources and sinks.

• •

Our AGN-only model produces larger HII bubbles as compared with our Galaxies-only model (see Figure 6), consistent with results from another semi-numerical model by Kulkarni et al. (2017) (Figure 6). This results in a larger 21cm power spectrum amplitude by $\sim\,\times\,1.5-2$ as compared with that from the galaxies-only model (see Figure 7).

• •

We examine a model which includes a 50% contribution from galaxies and AGN (assuming a constant QHOD). We find that this model can barely satisfy current $\tau_{e}$ and $x_{\rm HI}$ constraints. At early epochs, the galaxies contribution dominates the power spectrum, but during the latter stages of reionisation the quasars contribution is more significant, and the power spectrum deviates more substantially from the galaxies-only case.

• •

Future 21cm observations by LOFAR, HERA and SKA can discriminate between the constant QHOD AGN and Galaxies models during late reionisation at scales of $k<0.21$ Mpc*-1* ($>30$ Mpc), $k<0.53$ Mpc*-1* ($>12$ Mpc) and $k<1.66$ Mpc*-1* ($>4$ Mpc), respectively (see Figure 8). HERA and SKA will also be able to distinguish between our 50-50 and galaxies-only models, and thus potentially constrain the fractional contribution of AGN to reionisation.

• •

We have assumed an optimistic model for the AGN photon output rate by extrapolating to very low luminosities ($M_{1450}>-18$). There are also suggestions that G15 over-estimates the number of high-$z$ AGN and thus their total emissivity (e.g. Parsa et al., 2017). In either scenario, our claim that AGN cannot dominate reionisation is further strengthened.

It might still be possible that our AGN-only models could match all EoR observations simultaneously, if we relax some of these models’ assumptions. For instance, our AGN-only models depend on two parameters, namely, the photon escape fraction $f_{\mathrm{esc,AGN}}$ and the black hole formation efficiency $A$ (see Equation 4). These results are already obtained using 100% $f_{\mathrm{esc,AGN}}$, and it is clear that those models would reionise much earlier if we adopt $f_{\mathrm{esc,AGN}}\gg 100\%$, which is not physical. The $A$ parameter has been tuned to allow these models to reproduce the G15 data at z=5.75. If we choose not to calibrate these models to match the G15 data, we then can adopt larger $A$ which results in a smaller QHOD and earlier reionization. For example, the QHOD in the most faint magnitude bin at z=5.75 decreases from $\sim 2.7\%$ to $1.7\%$ as the $A$ increases from $5\times 10^{5}$ to $10^{6}$. This shows that allowing the AGN-only models to form more efficient black holes (large $A$) result in a fewer number (less QHOD) of them since the QHOD quantifies the fraction of active halos with AGN. From Figure 4, we see that the reionization in the constant LF AGN model (with $A=10^{6}$) starts earlier than in the constant QHOD AGN model (with $A=5\times 10^{5}$), but because there is no evolution in the source population (fixed LF), the model reionises much later than the constant QHOD. It is then clear that if we adopt larger $A$ to form more effeicent black holes, our AGN-only models will yield an early reionization. In this case, these models might produce consistent reionization histories and optical depths as compared with the observations. Given the large uncertainties in the ionising emissivity measurements by G15 and Becker & Bolton (2013), these models might still be consistent with these measurements’ 2-$\sigma$ level at $A>10^{6}$. We expect that the AGN clustering remains unaffected at a specific neutral fraction for larger $A$ values and hence the future 21cm observations can still discriminate between our AGN-only versus Galaxy-only models. The AGN-only models also assume that $M_{\rm bh}-v_{\rm cir}$ correlation (Equation 4) is valid even at high redshifts. This adds more uncertainties since such a correlation has only been measured in the local universe. More physically motivated and self-consistent AGN modelling at high redshift is clearly required to understand their formation and evolution as a function of redshift. Analogously to our stellar source model, one may allow the black hole emissivity to scale super-linearly with black hole mass (${R_{\rm ion,AGN}}\propto M^{C}_{\rm bh}$). This new parameter $C$ would indeed affect the AGN clustering and the corresponding ionised bubble sizes, as previously seen in our star-forming galaxies models (${R_{\rm ion,\star}}\propto M^{C}_{\rm h}$). We leave investigating this dependence for future works.

Improved high redshift AGN observations are clearly desirable in order to more robustly determine the AGN luminosity function, particularly at the low-luminosity end. The results of such observations can be incorporated straightforwardly into the framework we have presented here, and will provide better constraints on the contribution of AGN to reionisation. Our framework illustrates that the QHOD is an effective approach to evolve the number of AGN during the EoR. Our analysis indicates that despite there being potentially more faint AGN than previously believed, star-forming galaxies still dominate the neutral gas topology and ionising photon budget.

D’Aloisio et al. (2016), Mitra et al. (2016), and Oñorbe et al. (2017) have all independently demonstrated that AGN-only models overheat the IGM, inconsistent with Ly-$\alpha$ temperature measurements (Becker et al., 2011), due to the early onset of He II reionization. Finlator et al. (2016) further have shown that these models also over-ionise the metals when compared with observed metal absorption line measurements (D’Odorico et al., 2013). Our findings corroborate these results via a different approach.

AGN could still be important for reionisation because of their long-range heating effects owing to their harder emission spectrum, as well as for setting the shape of the metagalactic ionising flux that is important for interpreting metal-line absorption data (e.g. Finlator et al., 2016). Early AGN are in and of themselves interesting in order to understand the emergence of supermassive black holes particularly at early epochs. Our results here suggest that future 21cm experiments will have a key role to play in constraining the amount of AGN activity and its contribution to the metagalactic flux during the EoR.

acknowledgements

The authors acknowledge helpful discussions with Jonathan Pober and Jonathan Chardin. We thank Girish Kulkarni and Yuxiang Qin for making their models’ data available to us. SH is supported by the Deutscher Akademischer Austauschdienst (DAAD) Foundation. RD and SH are supported by the South African Research Chairs Initiative and the South African National Research Foundation. MGS is supported by the South African Square Kilometre Array Project and National Research Foundation. Part of this work was conducted at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293. RD acknowledges long-term visitor support provided by the Simons Foundation’s Center for Computational Astrophysics, as well as the Distinguished Visitor Program at Space Telescope Science Institute, where some of this work was conducted. Computations were performed at the cluster “Baltasar-Sete-Sois”, supported by the DyBHo-256667 ERC Starting Grant, and the University of the Western Cape’s “Pumbaa" cluster.