- Split View
-
Views
-
Cite
Cite
Yuxiang Qin, Simon J. Mutch, Gregory B. Poole, Chuanwu Liu, Paul W. Angel, Alan R. Duffy, Paul M. Geil, Andrei Mesinger, J. Stuart B. Wyithe, Dark-ages reionization and galaxy formation simulation – X. The small contribution of quasars to reionization, Monthly Notices of the Royal Astronomical Society, Volume 472, Issue 2, December 2017, Pages 2009–2027, https://doi.org/10.1093/mnras/stx1909
- Share Icon Share
Abstract
Motivated by recent measurements of the number density of faint AGN at high redshift, we investigate the contribution of quasars to reionization by tracking the growth of central supermassive black holes in an update of the Meraxes semi-analytic model. The model is calibrated against the observed stellar mass function at z ∼ 0.6–7, the black hole mass function at z ≲ 0.5, the global ionizing emissivity at z ∼ 2–5 and the Thomson scattering optical depth. The model reproduces a Magorrian relation in agreement with observations at z < 0.5 and predicts a decreasing black hole mass towards higher redshifts at fixed total stellar mass. With the implementation of an opening angle of 80 deg for quasar radiation, corresponding to an observable fraction of ∼23.4 per cent due to obscuration by dust, the model is able to reproduce the observed quasar luminosity function at z ∼ 0.6–6. The stellar light from galaxies hosting faint active galactic nucleus (AGN) contributes a significant or dominant fraction of the UV flux. At high redshift, the model is consistent with the bright end quasar luminosity function and suggests that the recent faint z ∼ 4 AGN sample compiled by Giallongo et al. (2015) includes a significant fraction of stellar light. Direct application of this luminosity function to the calculation of AGN ionizing emissivity consequently overestimates the number of ionizing photons produced by quasars by a factor of 3 at z ∼ 6. We conclude that quasars are unlikely to make a significant contribution to reionization.
1 INTRODUCTION
The epoch of reionization (EoR) is the phase of the Universe when neutral hydrogen in the intergalactic medium (IGM) was reionized. Star-forming galaxies at high redshift are believed to be one of the dominant sources of ionizing UV photons provided one assumes a high average escape fraction of Lyman continuum radiation (fesc,* ≳ 10 per cent) and extends the UV luminosity function to faint dwarf galaxies (Kuhlen & Faucher-Giguère 2012; Duffy et al. 2014; Feng et al. 2016; Mesinger, Greig & Sobacchi 2016). However, the value of fesc,* is very uncertain. Observations of star-forming galaxies at low redshift usually indicate a much lower escape fraction. For example, by measuring the ratio of Ly α to H β line emission, Ciardullo et al. (2014) derived an escape fraction of 4.4 per cent, while Matthee et al. (2016) measured a median escape fraction of 1.6 per cent using stacking of H α-selected galaxies. Moreover, some theoretical works also suggest a low fesc,* (Gnedin, Kravtsov & Chen 2008; Hassan et al. 2016; Sun & Furlanetto 2016).
In order to reconcile the difference between low-redshift observations and the photon budget at high redshift, some propose a rapid increase of fesc,* with redshift (Haardt & Madau 2012; Khaire et al. 2016; Price, Trac & Cen 2016; Sharma et al. 2016) and with decreasing mass (Paardekooper, Khochfar & Dalla 2013; Kimm & Cen 2014; Wise et al. 2014). This is supported by identifying local analogues of high-redshift galaxies and extrapolating the observed fesc,* using indicators such as the [O iii]/[O ii] ratio to high redshift (Faisst 2016 and references therein). On the other hand, additional contributors to reionization may also be present. For example, Ma et al. (2016) included a binary population into a set of radiative transfer cosmological simulations and found that they produced significantly more ionizing photons among the old stellar population compared to a model without binaries, reducing the requirement of high escape fraction. In addition, these photons produced at later times can escape from galaxies more easily since the local feedback efficiently clears out nearby gas, leading to a lower required escape fraction on average. With a high escape fraction (fesc,q ∼ 1, Barkana & Loeb 2001), quasars (active galactic nucleus, AGN) are potential contributors to reionization (Volonteri & Gnedin 2009; Fontanot et al. 2014; Madau & Haardt 2015; Mitra, Roy Choudhury & Ferrara 2015), despite their relatively low number. Recently, Giallongo et al. (2015) identified faint AGN candidates in the CANDELS GOODS-South field and suggested that there is a high number density of faint AGN at z = 4–6. These faint quasars provide a new source of reionization (Madau & Haardt 2015). However, it is still debated whether there are enough luminous quasars at high redshift to make a significant contribution (Bouwens et al. 2015; Weigel et al. 2015; Manti et al. 2017; Parsa, Dunlop & Mclure 2017) and whether the escape fraction of high-redshift, low-luminosity AGN is of the order of unity (Cristiani et al. 2016; Micheva, Iwata & Inoue 2017).
To explore the consequences for galaxy formation and reionization from theses faint quasars, this paper describes the addition of a population of evolving black holes to the Meraxes1 semi-analytic model of galaxy formation and reionization (Mutch et al. 2016). This new model enables a detailed exploration of the relative role of quasars during the EoR. The paper is organized as follows. We begin with a description of the semi-analytic model in Section 2, in which the black hole growth model is introduced in detail. We present the black hole properties in Section 3 and explore reionization from quasars in Sections 4 and 5. Conclusions are given in Section 6. In this work, we adopt cosmological parameters from the Planck 2015 results (Ωm, Ωb, ΩΛ, h, σ8, ns = 0.308, 0.0484, 0.692, 0.678, 0.815, 0.968; Planck Collaboration XIII 2016).
2 MODELLING BLACK HOLE GROWTH
Built on halo2 merger trees constructed from the Tiamat collisionless N-body simulation (Poole et al. 2016), the Meraxes semi-analytic model (Mutch et al. 2016) was specifically designed to study galaxy formation and reionization at high redshift. The model computes galaxy properties according to different astrophysical processes including gas infall, cooling, star formation, supernova feedback, metal enrichment, mergers and reionization. In order to properly track the evolution of galaxies during reionization, Tiamat provide 100 snapshots between z = 35 and 5 with a time interval of ∼11.1 Myr and 64 additional snapshots between z = 5 and 2 separated equally in units of Hubble time. The mass resolution of Tiamat is ∼2.64 × 106 h−1 M⊙ and the box size is 67.8 h−1Mpc. Additionally, in order to obtain information of more massive objects and lower redshifts, we also take advantages of the dark matter halo merger trees generated from the Tiamat-125-HR simulation (Poople et al. in prep). The Tiamat-125-HR simulation shares identical cosmology with Tiamat but has a lower mass resolution of 0.12 × 109 h−1 M⊙ in a bigger simulation volume, with side lengths equal to 125 h−1 Mpc. The Tiamat-125-HR trees are constructed down to z = 0.56 with the same snapshot separation strategy as the Tiamat trees.
In the following subsections, we briefly describe the galaxy formation in Meraxes and introduce the new implementation of black hole growth and feedback in detail. More details about the implemented galaxy formation physics can be found in Mutch et al. (2016).
2.1 Galaxy evolution
In addition, the metals produced by supernovae enrich the environment, which then enhances the cooling rate (see equation 2). Moreover, mergers drive strong turbulence and hence increase the possibility of star formation. Major mergers generally introduce more energetic bursts than minor mergers since they induce strong inflows and easily trigger bar-like instabilities in the cloud (Somerville, Primack & Faber 2001). Therefore, following mergers, Meraxes also includes a starburst mechanism. Reionization feedback will be introduced in Section 4.1. There are more details of the semi-analytic model in Mutch et al. (2016). This work extends the current model with a detailed black hole growth prescription based on Croton et al. (2016) and is described in the following subsection.
2.2 Black hole growth
In the updated model, every newly formed galaxy is seeded with a central black hole5 of mass 1000 h−1 M⊙. The two gas reservoirs in the galaxy (i.e. hot and cold) lead to two different black hole growth scenarios, termed radio and quasar modes (Croton et al. 2016). In the normal quiescent state, black holes only accrete mass from the hot gas reservoir, resulting in radio emission in the centre of galaxy. However, mergers trigger rapid accretion on to the black hole from the cold gas disc, causing them to radiate as quasars. In this work, we do not distinguish AGN with different types and refer to them all as quasars unless specified otherwise.
2.2.1 Accretion of hot gas
We note that at high redshift, accreted hot gas does not contribute significantly to black hole growth, with the accreted mass during the radio mode typically being only ∼0.1 per cent of the mass from the quasar mode (see the Appendix) that is introduced in the next section. Additionally, because the radio mode accretion rate is approximately 3 orders of magnitude smaller than the Eddington accretion rate, we ignore the energy due to the radio mode for the calculation of the quasar luminosity.
2.2.2 Accretion of cold gas
- Δmbh,max < mEdd: In this case, there is inadequate mass to feed the central black hole at the Eddington rate for the entire time step. The duration of accretion can be calculated throughWhen the quasar is observed at a random time tobs, the bolometric luminosity is(21)\begin{equation} t_\mathrm{acc} = \ln \left(\frac{\Delta {m}_\mathrm{bh,cold}}{m_\mathrm{bh}} +1\right)\times \frac{\eta t_\mathrm{Edd}}{\epsilon }. \end{equation}when tobs ≤ tacc, otherwise Lbol = 0.(22)\begin{equation} L_\mathrm{bol} \equiv \epsilon m_\mathrm{bh}|_{t=t_\mathrm{obs}} \frac{c^2}{t_\mathrm{Edd}} = \epsilon m_\mathrm{bh}\exp \left(\frac{\epsilon t_\mathrm{obs}}{\eta t_\mathrm{Edd}}\right)\times \frac{c^2}{t_\mathrm{Edd}}, \end{equation}
Δmbh,max ≥ mEdd: In this case, the merger event delivers sufficient cold gas into the accretion disc. In the case of Δmbh,max > mEdd, instead of consuming this instantaneously, causing a super-Eddington accretion event, some of the mass is accreted by the central black hole, limited by the Eddington rate, while the rest, |$\Delta {m}_\mathrm{bh,max}^\prime = \Delta {m}_\mathrm{bh,max} - m_\mathrm{Edd}$|, is stored in the accretion disc to be consumed in the next time step. Similarly, when this quasar is observed at tobs, the bolometric luminosity can be calculated using equation (22).
It is suggested that during mergers, black holes undergo rapid accretion for a certain time period, which is followed by a long quiescent phase (Hopkins et al. 2005a,b,c,d). The assumption that black holes are either accreting at the Eddington rate (ε = 1, see equation 16) or stay quiescent has been shown to provide a good description of black hole growth for the majority of black holes at high redshift (Bonoli et al. 2009).
The energy injected into galactic gas during the quasar mode is given by κqηΔmbh, coldc2, where κq represents the mass coupling factor in the quasar mode. Unlike the radio mode, this energy generates a wind, which heats the gas in the cold disc into the hot reservoir. Depending on the amount of energy provided by the quasar, the wind can further unbind and eject the hot gas in a manner similar to the stellar feedback prescription presented in Section 2.1.
2.3 Quasar luminosity
In order to compare the predicted black hole population in our model with observations, the intrinsic B band and UV 1450 Å band luminosities of quasars are calculated as follows:
- We calculate the bolometric magnitude through(23)\begin{equation} M_\mathrm{bol} = 4.74 - 2.5\log _{10}\left(\frac{L_\mathrm{bol}}{\mathrm{L_{\odot }}}\right). \end{equation}
- We calculate the B-band magnitude in the Vega magnitude system using the bolometric correction proposed by Hopkins, Richards & Hernquist (2007),where(24)\begin{equation} M_\mathrm{bol} - M_\mathrm{B} = -2.5\log _{10}k_\mathrm{B}, \end{equation}(25)\begin{eqnarray} k_\mathrm{B} \equiv \frac{L_\mathrm{bol}}{L_\mathrm{B}}= 6.25\left(\frac{L_\mathrm{bol}}{10^{10}\,\,\mathrm{L_{\odot }}}\right)^{-0.37} + 9.00\left(\frac{L_\mathrm{bol}}{10^{10}\,\,\mathrm{L_{\odot }}}\right)^{-0.012}.\nonumber \\ \end{eqnarray}
- We convert the B-band magnitude from the Vega system to the AB system following Glikman et al. (2010)(26)\begin{equation} M_\mathrm{AB,B} - M_\mathrm{B} = -0.09. \end{equation}
- We extrapolate the B-band magnitude of which the effective wavelength is 4344 Å (Blanton & Roweis 2007) to the 1450 Å magnitude, assuming the quasar continuum between 1450 and 4344 Å has a power-law slope of αq, optical = 0.44 relative to wavelength (Schirber & Bullock 2003). Thus,(27)\begin{eqnarray} M_{1450} &=& M_\mathrm{AB,B} - 2.5\log _{10}\left(\frac{1450 \mathrm{\,\,\mathring{\rm A} }}{4344\mathrm{\,\,\mathring{\rm A} }}\right)^{\alpha _{q,\mathrm{optical}}} = M_\mathrm{AB,B}+0.524.\nonumber \\ \end{eqnarray}
We will further discuss the quasar luminosity function in Section 3.3.
3 GALAXY AND BLACK HOLE PROPERTIES
We summarize the relevant model parameters in Table 1 compared to the original value adopted in M16. All other parameters remain the same as M16. In this work, we constrain9 our model against:
the observed evolution of the galaxy stellar mass function (Pozzetti et al. 2007; Drory et al. 2009; Marchesini et al. 2009; González et al. 2011; Mortlock et al. 2011; Santini et al. 2012; Ilbert et al. 2013; Muzzin et al. 2013; Duncan et al. 2014; Tomczak et al. 2014; Grazian et al. 2015; Song et al. 2016; Huertas-Company et al. 2016; Stefanon et al. 2017; Davidzon et al. 2017) between z ∼ 0.6 and 7. We note that the observed stellar mass functions, based on a diet Salpeter IMF, a Chabrier (2003) IMF or a Kroupa (2001) IMF, were all converted into a standard Salpeter IMF by adding −0.15, 0.22 or 0.18 dex, respectively, to the logarithm of the stellar masses;
the observed black hole mass function (Graham et al. 2007; Shankar, Weinberg & Miralda-Escudé 2009; Vika et al. 2009; Davis et al. 2014; Mutlu-Pakdil, Seigar & Davis 2016) and Magorrian relation (Thornton et al. 2008; Jiang, Greene & Ho 2011; Mathur et al. 2012; Jiang et al. 2013; Reines, Greene & Geha 2013; Scott, Graham & Schombert 2013; Busch et al. 2014; Sanghvi et al. 2014; Yuan et al. 2014) at low redshift (z ≲ 0.5);
the latest integrated free electron Thomson scattering optical depth measurement (Planck Collaboration XLVI 2016);
the predicted global ionizing emissivity from Ly α opacities (Becker & Bolton 2013).
Parameter . | Section . | Equation . | Description . | Fiducial . | M16BH . | M16 . |
---|---|---|---|---|---|---|
αsf | 2.1 | 6 | Star formation efficiency | 0.08 | 0.03 | 0.03 |
αenergy | 2.1 | 7 | Energy coupling efficiency normalization | 1.0 | 0.5 | 0.5 |
αmass | 2.1 | 9 | Mass loading normalization | 15.0 | 9.0 | 9.0 |
η | 2.2.1 | 15 | Black hole efficiency of converting mass to energy | 0.06 | 0.06 | - |
ε | 2.2.1 | 16 | Eddington ratio | 1.0 | 1.0 | - |
kh | 2.2.1 | 14 | Black hole growth efficiency for the radio mode | 0.3 | 0.1 | - |
κr | 2.2.2 | 16 | Black hole feedback efficiency for the radio mode | 1.0 | 1.0 | - |
kc | 2.2.2 | 19 | Black hole growth efficiency for the quasar mode | 0.05 | 0.05 | - |
κq | 2.2.2 | – | Black hole feedback efficiency for the quasar mode | 0.0005 | 0.0005 | - |
Nγ,* | 4.1 | 28 | Mean number of ionizing photons produced per stellar baryon | 4,000 | 4,000 | 4,000 |
fesc,q | 4.1 | 28 | Ionizing photon escape fraction for quasars | 1.0a | 0.0 | - |
fesc,* | 4.1 | 28 | Ionizing photon escape fraction for the stellar component | fesc, *, zb, c | fesc, *, M16d | fesc, *, M16 |
Parameter . | Section . | Equation . | Description . | Fiducial . | M16BH . | M16 . |
---|---|---|---|---|---|---|
αsf | 2.1 | 6 | Star formation efficiency | 0.08 | 0.03 | 0.03 |
αenergy | 2.1 | 7 | Energy coupling efficiency normalization | 1.0 | 0.5 | 0.5 |
αmass | 2.1 | 9 | Mass loading normalization | 15.0 | 9.0 | 9.0 |
η | 2.2.1 | 15 | Black hole efficiency of converting mass to energy | 0.06 | 0.06 | - |
ε | 2.2.1 | 16 | Eddington ratio | 1.0 | 1.0 | - |
kh | 2.2.1 | 14 | Black hole growth efficiency for the radio mode | 0.3 | 0.1 | - |
κr | 2.2.2 | 16 | Black hole feedback efficiency for the radio mode | 1.0 | 1.0 | - |
kc | 2.2.2 | 19 | Black hole growth efficiency for the quasar mode | 0.05 | 0.05 | - |
κq | 2.2.2 | – | Black hole feedback efficiency for the quasar mode | 0.0005 | 0.0005 | - |
Nγ,* | 4.1 | 28 | Mean number of ionizing photons produced per stellar baryon | 4,000 | 4,000 | 4,000 |
fesc,q | 4.1 | 28 | Ionizing photon escape fraction for quasars | 1.0a | 0.0 | - |
fesc,* | 4.1 | 28 | Ionizing photon escape fraction for the stellar component | fesc, *, zb, c | fesc, *, M16d | fesc, *, M16 |
afesc,q = 0 is adopted for the StellarReion model in which stars are the only reionization source.
bfesc,* = 0 is adopted for the QuasarReion model in which quasars are the only reionization source.
|$^c{}f_\mathrm{esc,*,z} =\min [0.06\times (\frac{1+z}{6})^{0.5},1.0]$| for the fiducial and StellarReion models.
d|$f_\mathrm{esc,*,M16} = \min [0.04\times (\frac{1+z}{6})^{2.5},1.0]$|.
Parameter . | Section . | Equation . | Description . | Fiducial . | M16BH . | M16 . |
---|---|---|---|---|---|---|
αsf | 2.1 | 6 | Star formation efficiency | 0.08 | 0.03 | 0.03 |
αenergy | 2.1 | 7 | Energy coupling efficiency normalization | 1.0 | 0.5 | 0.5 |
αmass | 2.1 | 9 | Mass loading normalization | 15.0 | 9.0 | 9.0 |
η | 2.2.1 | 15 | Black hole efficiency of converting mass to energy | 0.06 | 0.06 | - |
ε | 2.2.1 | 16 | Eddington ratio | 1.0 | 1.0 | - |
kh | 2.2.1 | 14 | Black hole growth efficiency for the radio mode | 0.3 | 0.1 | - |
κr | 2.2.2 | 16 | Black hole feedback efficiency for the radio mode | 1.0 | 1.0 | - |
kc | 2.2.2 | 19 | Black hole growth efficiency for the quasar mode | 0.05 | 0.05 | - |
κq | 2.2.2 | – | Black hole feedback efficiency for the quasar mode | 0.0005 | 0.0005 | - |
Nγ,* | 4.1 | 28 | Mean number of ionizing photons produced per stellar baryon | 4,000 | 4,000 | 4,000 |
fesc,q | 4.1 | 28 | Ionizing photon escape fraction for quasars | 1.0a | 0.0 | - |
fesc,* | 4.1 | 28 | Ionizing photon escape fraction for the stellar component | fesc, *, zb, c | fesc, *, M16d | fesc, *, M16 |
Parameter . | Section . | Equation . | Description . | Fiducial . | M16BH . | M16 . |
---|---|---|---|---|---|---|
αsf | 2.1 | 6 | Star formation efficiency | 0.08 | 0.03 | 0.03 |
αenergy | 2.1 | 7 | Energy coupling efficiency normalization | 1.0 | 0.5 | 0.5 |
αmass | 2.1 | 9 | Mass loading normalization | 15.0 | 9.0 | 9.0 |
η | 2.2.1 | 15 | Black hole efficiency of converting mass to energy | 0.06 | 0.06 | - |
ε | 2.2.1 | 16 | Eddington ratio | 1.0 | 1.0 | - |
kh | 2.2.1 | 14 | Black hole growth efficiency for the radio mode | 0.3 | 0.1 | - |
κr | 2.2.2 | 16 | Black hole feedback efficiency for the radio mode | 1.0 | 1.0 | - |
kc | 2.2.2 | 19 | Black hole growth efficiency for the quasar mode | 0.05 | 0.05 | - |
κq | 2.2.2 | – | Black hole feedback efficiency for the quasar mode | 0.0005 | 0.0005 | - |
Nγ,* | 4.1 | 28 | Mean number of ionizing photons produced per stellar baryon | 4,000 | 4,000 | 4,000 |
fesc,q | 4.1 | 28 | Ionizing photon escape fraction for quasars | 1.0a | 0.0 | - |
fesc,* | 4.1 | 28 | Ionizing photon escape fraction for the stellar component | fesc, *, zb, c | fesc, *, M16d | fesc, *, M16 |
afesc,q = 0 is adopted for the StellarReion model in which stars are the only reionization source.
bfesc,* = 0 is adopted for the QuasarReion model in which quasars are the only reionization source.
|$^c{}f_\mathrm{esc,*,z} =\min [0.06\times (\frac{1+z}{6})^{0.5},1.0]$| for the fiducial and StellarReion models.
d|$f_\mathrm{esc,*,M16} = \min [0.04\times (\frac{1+z}{6})^{2.5},1.0]$|.
3.1 Black hole properties
The black hole mass functions at z ∼ 8.0–0.6 are shown with different colours in the left-hand panel of Fig. 1. The results calculated using the Tiamat and Tiamat-125-HR trees are shown with thick and thin lines, respectively. The shaded regions represent the 1σ Poisson uncertainties.10 Estimates of the local black hole mass function are shown with points and grey shaded regions. We see that the discrepancy between various observations is substantial. This is a result of inconsistent correlations between black hole mass and observable quantities, such as the Sérsic indices, bulge velocity dispersion or luminosity and galaxy geometry, and from the intrinsic scatter of these adopted scaling relations. Extrapolations of these observed scaling relations have impacts on the black hole mass function at the high-mass end, while different treatments of the spiral galaxy bulge can significantly change the low-mass end (Shankar et al. 2009). In this work, the model is therefore calibrated against the black hole mass function between 107.5 M⊙ and 109 M⊙. In the bottom left-hand panel of Fig. 1, we see that the mass function converges at lower redshifts above a black hole mass of 106 M⊙ (shown as the vertical dotted line). The different mass resolutions of Tiamat and Tiamat-125-HR result in different merger rates, especially when approaching the resolution limit. At high redshift, because the growth of black hole is dominated by the merger triggered quasar mode, the number density of small black holes is relatively lower in the Tiamat-125-HR result (e.g. comparing the z = 8.0 thick and thin lines). At low redshift, z ∼ 0.6, the model is in agreement with the observational estimations.
The middle panel of Fig. 1 shows the relation between black hole mass and stellar mass (the Magorrian relation11, Magorrian et al. 1998). The 2D histogram indicates the distribution of galaxies12 in logarithm from the fiducial model using the Tiamat-125-HR halo merger trees at z ∼ 0.6, while the solid line represents the mean. The Magorrian relations at z = 2, 5 and 7 from Tiamat-125-HR are also shown with dash–dotted, dashed and dash–dot–dotted lines, respectively. The right bottom subplot shows the z = 2, 5 and 7 Magorrian relations of the Tiamat result with thick lines compared with the z ∼ 0.6 Tiamat-125-HR Magorrian relation. Observations from the local Universe are indicated with different symbols (Thornton et al. 2008; Jiang et al. 2011; Mathur et al. 2012; Jiang et al. 2013; Reines et al. 2013; Scott et al. 2013; Busch et al. 2014; Sanghvi et al. 2014; Yuan et al. 2014). We see that the model predicts a similar Magorrian relation at z ∼ 0.6 compared to the local observations and we find an increasing normalization towards lower redshifts in the mass range of 1010 M⊙ < M* < 1012 M⊙.
The evolution of the Magorrian relation in our model is due to the black hole and stellar mass evolving with the underlying dark matter halo mass differently. In M16, we have shown that the median relation between stellar mass and virial mass does not evolve in our model and can be described by |$M_*\propto M_\mathrm{vir}^{7/5}$| in the range of 108 M⊙ < M* < 1011 M⊙, which is supported by a simple analytic model of supernova energy conversation (Wyithe & Loeb 2003). On the other hand, haloes with a given virial mass host less massive black holes at earlier times in our model (see the right-hand panel of Fig. 1). These result in an increasing normalization of the Magorrian relation towards lower redshifts.
From the right-hand panel of Fig. 1, we see that the Mbh–Mvir scaling relation does not get suppressed in massive haloes (at least to Mvir ∼ 1014 M⊙ or Mbh ∼ 109 M⊙). Both black holes and stars grow from the cold gas disc.13 However, AGN feedback significantly suppresses the cooling flow in massive galaxies (Croton et al. 2006), preventing stellar mass from growing. On the other hand, black holes are able to continue accreting until there is enough energy in feedback to overcome the halo potential and unbind the gas (Booth & Schaye 2010). Because of these, the slope becomes steeper in the Magorrian relation at Mbh > 108 M⊙.
3.2 Galaxy properties
Fig. 2 presents the galaxy stellar mass functions from our fiducial model for comparison with the available observational data (Pozzetti et al. 2007; Drory et al. 2009; Marchesini et al. 2009; González et al. 2011; Mortlock et al. 2011; Santini et al. 2012; Ilbert et al. 2013; Muzzin et al. 2013; Duncan et al. 2014; Tomczak et al. 2014; Grazian et al. 2015; Song et al. 2016; Huertas-Company et al. 2016; Stefanon et al. 2017; Davidzon et al. 2017) at redshifts 7 to ∼0.6. The results calculated using the Tiamat and Tiamat-125-HR trees are shown with thick and thin lines, respectively. The shaded regions represent the 1σ Poisson uncertainties. We see that the fiducial model is able to reproduce the observed galaxy stellar mass function across the redshift range of z = 7–0.6.
We also present two additional models in Fig. 2, M16 and M16BH. M16 adopts identical parameters as the redshift varying fesc,* model (|$f_\mathrm{esc,*,M16} = \min \left[0.04\times \left(\frac{1+z}{6}\right)^{2.5},1.0\right]$|, see Table 1) presented in M16. This model is able to reproduce the evolution of the stellar mass function at high redshift (z > 5). The redshift dependence of fesc, *, M16 was chosen to simultaneously reproduce the normalization and flat slope of the McQuinn, Peng Oh & Faucher-Giguère (2011) emissivity measurement at z ∼ 5 and the Planck 2015 optical depth measurement (Planck Collaboration XIII 2016, hereafter Planck15). Details of the M16 model at high redshift (z > 5) can be found in M16. In this work, we extend the model14 to lower redshifts and find that, without AGN feedback regulating galaxy formation, the model fails to reproduce the observed stellar mass functions at z < 2, especially in larger mass ranges (M* > 1011 M⊙). However, when AGN feedback is implemented, shown as the M16BH model, the model shows better agreement with observations (Croton et al. 2006). These two models also suggest that radio-mode feedback does not play a significant role in galaxy formation during the EoR, and because reionization is dominated by low-mass galaxies (Liu et al. 2016, see also Section 5), AGN feedback is expected to have no significant impact on reionization.15
With respect to M16, the fiducial model presented in this work employs a stronger star formation efficiency (αsf) with maximized supernova feedback (αenergy and αmass) and more intense radio mode feedback (kh), in order to gain better agreement with the observed stellar mass function in the intermediate mass range (109 M⊙ < M* < 1011 M⊙) at 1 < z < 2.
3.3 Quasar luminosity function
In the previous two subsections, we presented the predicted properties of the central massive black holes and host galaxies, together with their correlations. These show that our model is able to produce the evolution of the galaxy stellar mass function over a large time-scale (z ∼ 7–0.6), as well as the observed black hole properties at low redshift. Before we start exploring the contribution of quasars to reionization, we discuss the quasar luminosity function in this section.
Semi-analytic models can predict quasar bolometric luminosity through the modelled central black hole mass. In order to compare with observations, three corrections are usually required: bolometric corrections, duty cycles and obscured fractions. In our model, we calculate the duty cycle self-consistently by assuming Eddington accretion of available gas and a random observation time between snapshots (see Section 2.2.2). Black holes that are inactive when observed are not considered in the census. We convert the total luminosity into a particular band using the bolometric correction of Hopkins et al. (2007) (see Section 2.3). Finally, we need to take the obscuration due to the presence of a dusty torus surrounding the AGN into account. In practice, each quasar is weighted by |$1-\cos \frac{\theta }{2}$|, where θ represents the opening angle of AGN radiation, during the calculation of the quasar luminosity function. We note that the dependence of θ or the obscured fraction (the ratio of obscured to unobscured AGN) is very complicated. It has been suggested that θ may depend on wavelength, redshift and luminosity (Elvis et al. 1994; Hopkins et al. 2007). For simplicity, we only consider constant θ. The opening angle is a free parameter in our model, chosen to reproduce the quasar luminosity function amplitude. In this work, θ = 80 deg is adopted. This corresponds16 to a fraction of visible objects,17fobs, to be ∼23.4 per cent, in agreement with Hopkins et al. (2007), who suggests a luminosity-dependent observable fraction, fobs = 0.26(Lbol/1012.4 L⊙)0.082 for the B band.
We present the UV 1450 Å luminosity functions of quasars at redshift 6, 5, 4 and 3, and the B-band luminosity functions at z ∼ 2, 1.5, 1.3 and 0.6 in Fig. 3, compared to the observational data from Wolf et al. (2003), Hunt et al. (2004), Richards et al. (2005), Dijkstra et al. (2006), Bongiorno et al. (2007), Croom et al. (2009), Willott et al. (2010), Glikman et al. (2011), Shen & Kelly (2012), Masters et al. (2012), McGreer et al. (2013), Palanque-Delabrouille et al. (2013), Giallongo et al. (2015) and Jiang et al. (2016). The results calculated using the Tiamat and Tiamat-125-HR trees are shown with thick and thin lines, respectively. The shaded regions represent the 95 per cent confidence intervals around the mean using 100 000 bootstrap re-samples (due to the random number tobs in equation 22). Similarly to the black hole mass functions, the Tiamat-125-HR results show a lower number density at high redshift compared to the Tiamat results, and they converge at lower redshifts (z ≲ 4). We see that the model shows good agreement with observations across a large redshift range (z ∼ 6–0.6). At high redshift, the model is consistent with the samples of bright quasars (Willott et al. 2010; Shen & Kelly 2012; McGreer et al. 2013; Jiang et al. 2016), while it predicts a lower number density of faint quasars compared to the Giallongo et al. (2015, hereafter G15) sample. On the other hand, the model produces a significant number of faint quasars. The turnover (not shown here) of the predicted UV luminosity function is around M1450 ∼ −11, which is much fainter than the observed faintest quasars, M1450 ∼ −18 (G15). In addition, the observed population of bright quasars (M1450 < −23) at high redshift is not present in our model due to our limited simulation volume. All of these factors may have an impact on the contribution of AGN to the reionization history, which is discussed in Section 4.
3.4 The luminosity function of galaxies with AGN
The large number of faint AGN identified by G15 has prompted renewed discussion of the contribution of quasars to reionization (Madau & Haardt 2015; Mitra et al. 2015; Kulkarni et al. 2017). We therefore further discuss the low number density in the faint end predicted by the model compared with the G15 data.
When constructing an observed AGN luminosity function, the sample is typically identified spectroscopically or via colour–colour selections if spectra are not available, following which contamination from host galaxies is removed. For example, some observers model the surface brightness distributions of host galaxies and fit galaxies (e.g. Sérsic profiles) and point sources (using point-spread functions, Dunlop et al. 2003; Peng et al. 2006; Du et al. 2014; Martínez-Paredes et al. 2017) to the images.18 Others examined the SED using a combination of AGN and galaxy emission with possible extinction of the AGN flux when spectroscopic data are available (Bongiorno et al. 2007; Croom et al. 2009; Masters et al. 2012; Mechtley et al. 2012; Lyu, Rieke & Alberts 2016). In most cases, ignoring the contribution from host galaxies does not make a significant difference to bright quasars (Hopkins et al. 2007). Therefore, some AGN samples do not exclude stellar light (Wolf et al. 2003, Hunt et al. 2004, Richards et al. 2005, Willott et al. 2010, Shen & Kelly 2012, McGreer et al. 2013, Palanque-Delabrouille et al. 2013; G15). However, this may not be the case for faint AGN. In particular, the G15 AGN sample is selected using X-ray activity and no AGN-galaxy separation is possible. Thus, the total UV luminosity may have a large fraction of stellar light, suggesting that the G15 sample may be potentially impacted by stellar light contamination. Despite the recent claim that only 12 of the 22 reported X-ray detections in G15 are high-redshift AGN (Parsa et al. 2017), this conjecture is supported by Ricci et al. (2017), who used X-ray observations as a proxy and derived the quasar UV luminosity function down to much fainter ranges. They showed that the luminosity function is in agreement with UV/optical observations (e.g. Croom et al. 2009; Glikman et al. 2010; Masters et al. 2012; Palanque-Delabrouille et al. 2013), and have much lower amplitudes than the G15 results, and that the high number density of faint AGN in G15 can be explained by the contribution from the luminosity of host galaxy with M1450 ∼ −20.
In order to account for this, we calculate the galaxy UV luminosity by integrating model SEDs based on the modelled star formation history (Liu et al. 2016) and add the stellar light to mimic the observed total UV luminosity19 (AGN+galaxy). We present the resulting luminosity function of faint objects (M1450 ∼ −18 to −22) in Fig. 4 with dashed lines. The result is calculated using the Tiamat halo merger trees and is shown with shaded regions representing the 95 per cent confidence intervals around the mean using 100 000 bootstrap re-samples. The luminosity function calculated using only the AGN light (AGN only, as shown in Fig. 3) is shown as solid lines for comparison. We see that including stellar light can significantly increase the number density of faint AGN inferred at high redshift, by up to ∼1 order of magnitude. The fitting functions provided by G15 as shown with thin black dashed lines in Fig. 4 are more consistent with the AGN+galaxy luminosity function. If this is the case, the estimated emissivity at high redshift based on G15 is likely overestimated.
4 REIONIZATION FROM QUASARS
Motivated by the G15 sample, which suggests a numerous population of faint AGN at z = 4–6 (G15), Madau & Haardt (2015) extrapolated the emissivity calculated from G15 to higher redshifts and assessed a model of reionization, in which quasars are the dominant ionizing sources. They found that due to the high escape fraction of ionizing photons produced by those luminous objects, quasars are able to ionize the neutral hydrogen by z ∼ 5.7 if the high emissivity of quasars derived at z ∼ 5 continues to higher redshifts. Later, Mitra et al. (2015) revisited the model with a revised extrapolation and also found that quasars have a significant role during the EoR. However, their analysis still prefers models with a non-zero escape fraction of ∼12 per cent from galaxies. It is worth noting that with different formalism, Manti et al. (2017) fit the observed quasar UV luminosity function at z = 0.5–6.5, including the G15 sample. They recalculated the emissivity by integrating the luminosity function and confirmed a large number of ionizing photons from quasars at high redshift using the Schechter luminosity function. However, the fitting result using a double power law presents a rapidly decreasing emissivity at z > 6, suggesting that the extrapolated high-redshift quasar emissivity is strongly dependent on the assumed shape of the quasar luminosity function. Taking advantage of the Meraxes semi-analytic model with 21cmFAST (Mesinger, Furlanetto & Cen 2011), we investigate the contribution of quasars to reionization, within a frame work that accounts for black hole growth and feedback on star formation.20
4.1 Reionization model
4.2 Reionization history
There are two additional models in the top panel of Fig. 5 as well as in Fig. 6: StellarReion and QuasarReion. The ionizing source in the StellarReion model is only stars, with |$f_\mathrm{esc,*} = \min \left[0.06\times \left(\frac{1+z}{6}\right)^{0.5},1.0\right]$| and fesc,q = 0, while quasars are the only reionization contributor in the QuasarReion model, with fesc,q = 1 and fesc,* = 0 (see Table 1). Note that only changing the feedback from reionization has little impact on the stellar mass function (see M16), the quasar luminosity function or the Magorrian relation. By preventing gas infall, reionization only affects less massive objects in our model.
In Fig. 5, we see that the emissivity of quasars grows rapidly in the QuasarReion model (also in the fiducial model) by a factor of 10 from z ∼ 7–5. However, if quasars are the only reionization contributor, even with fesc,q = 1 the number of ionizing photons cannot reach the BB13 data. Moreover, due to the deficiency in the photon budget at high redshift, quasars can only start reionization at z ∼ 6 resulting in an end at z ∼ 3. Together with the predicted optical depth, our model rules out the quasar-only reionization scenario. We note that one may recalibrate the model with a more efficient black hole growth rate at high redshift (e.g. by incorporating a redshift dependence in the equation 19, see Bonoli et al. 2009), in order to match the G15 luminosity function and the estimated emissivity. Mitra et al. (2015) also suggest that if G15 emissivity is correct, quasar-only reionization is possible and it results in a small value of τe due to the rapid evolution of the Lyman-limit systems. However, simultaneously matching the model with the G15 faint AGN luminosity function and the other observations of bright systems at high redshift is difficult. For example, comparing to observations at z ∼ 6, our models produce a flatter luminosity function, which is more consistent with the bright quasar sample. This suggests that a mass-dependent black hole growth efficiency (e.g. kc, see Section 2.2.2) would be required, in order to steepen the luminosity function and produce more small quasars. In the following sections, we explore the relative contribution of quasars to reionization based only on the presented black hole growth model.
Comparing the neutral hydrogen fraction and the optical depth between the fiducial and StellarReion models also suggests that quasars do not have a significant role in reionization in this model. Their contribution helps finish reionization earlier by Δz ≲ 0.1 and decreases the optical depth by less than 10 per cent.
5 DISCUSSION
Our models are calibrated against the black hole – galaxy scaling relation and quasar luminosity function, in order to reproduce a realistic AGN catalogue for the study of the contribution of quasars to hydrogen reionization. However, it has recently been suggested that the black hole sample used to derive scaling relations is likely different from the entire population, leading to a selection bias25 (Bernardi et al. 2007). For instance, using Monte Carlo simulations Shankar et al. (2016) recovered the intrinsic scaling relation assuming the selection bias comes from unresolved black holes (e.g. for galaxies with a given velocity dispersion, σ*, it is more difficult to resolve smaller black holes; Batcheldor et al. 2010) and showed that such a bias can lead to factors of ≥3 and ∼50–100 higher normalizations of the Mbh–σ* and Mbh–M* relations, respectively. We note that a Magorrian relation with a smaller normalization can be achieved in our model using a smaller black hole growth efficiency (i.e. kc in equation 19). In this case, the reconstructed black hole population becomes less massive, more efficient AGN feedback (e.g. kh in equation 14) and a larger fraction of observable AGN in the UV band are required to simultaneously reproduce the observed stellar mass function and quasar luminosity function. With the model calibrated against the quasar luminosity function, the black hole–galaxy scaling relation is coupled with the observable AGN fraction – a lower Magorrian relation requires a larger fraction of observable AGN. We note that the total emissivity of quasars is integrated from the luminosity function. Therefore, scaling relations do not have a significant impact on reionization in this work. Noting the difficulty of observationally determining the fraction of obscured AGN and the large uncertainties in the black hole–galaxy scaling relation, in this section we use a range of models which predict similar Magorrian relations as shown in Fig. 1 and explore the contribution of quasars to reionization.
5.1 A larger opening angle
We have shown that with an opening angle of 80 deg, the model is able to reproduce the observed quasar luminosity function from z ∼ 6–0.6 (see Fig. 3). However, at high redshift, the model predicts significant stellar contribution to UV flux in the G15 sample (see Fig. 4), and consequently less ionizing photons from AGN. Based on this, we find that quasars do not have a significant role during EoR.
In the top panel of Fig. 7, we show the estimated total emissivity from G15 with the model proposed by Haardt & Madau (2012); Madau & Haardt (2015) and Mitra et al. (2015). The modelled emissivity from QuasarReion is shown for comparison. Our quasar reionization-only model (QuasarReion) predicts lower emissivities compared to these two estimations, with only a third of the G15 value at z ∼ 6. Consequently, in disagreement with Madau & Haardt (2015) and Mitra et al. (2015), we conclude that quasars cannot be the dominate sources during reionization. We could increase the emissivity by excluding the obscuration from dust (setting θ = 180, shown as QuasarReion_nodust in Fig. 7), which gives a closer emissivity compared to the G15 estimation and in agreement with the model proposed by Mitra et al. (2015). In this model, quasars have a more significant role during the EoR and can reionize the IGM alone by z ∼ 4.5. However, this model overestimates the number density of bright- and low-redshift AGN, leading to an incorrect evolution of the quasar luminosity function. A lower fraction of observable AGN, fobs, towards brighter luminosities and lower redshifts is required to solve this conflict. However, observations suggest the opposite trend in optical, infrared and X-ray bands (Hopkins et al. 2007) and more constraints are required to clearly establish a fobs–z relation.
5.2 Escape fraction of ionizing photons from quasars
In this section, we explore possible combinations of stars and quasars that could result in an overall photon budget at z > 5 consistent with the observed optical depth and ionizing flux at z ∼ 2–5. In the Section 4, as well as in M16 where only galaxies are considered, we demonstrated the requirement of an evolving escape fraction for stars to explain the observed emissivity at z ∼ 5. Noting this requirement, in this section, we assume constant escape fractions both for simplicity and to ease interpretation.
Motivated by the recent claim that the escape fraction of low-luminosity AGN is possibly less than unity at high redshift (Micheva et al. 2017), we run Meraxes with different combinations of fesc,* and fesc,q, without any changes to the other parameters. In the top panel of Fig. 8, the left panel shows the Thomson scattering optical depth, τe. For comparison, shaded regions are shown corresponding to the best fit and 1σ range of the Planck16 measurements. Based on the BB13 data at z ∼ 2–5 and the Planck16 measurement, the top right-hand two panels show the 68, 90 and 99 per cent confidence contours on each parameter of the best fit via the standard minimum-χ2 technique. The 2D histogram shows the distribution of the ratios of quasar emissivity to stellar emissivity at z ∼ 6. We see that a lower escape fraction of ionizing photons from stars, fesc,*, requires a higher contribution from quasars, in order to reach the observational constraint. This also returns a higher ratio of quasar emissivity to stellar emissivity. However, because there is not a significant number of quasars at high redshift, changing fesc,q has little impact to the optical depth. Through the best-fitting contours, we see that if the escape fraction of ionizing photons from stars is only a few per cent (<5 per cent; Ciardullo et al. 2014; Matthee et al. 2016b), the model requires fesc,q ∼ 1.0.
5.3 An evolving escape fraction
The bottom panels of Fig. 8 show the optical depth and the best-fitting confidence limits as functions of the normalization of fesc,*(z), fesc,*|z = 5 and the scaling, β, when fesc,q = 1. We see that because a larger scaling suppresses the escape fraction at lower redshifts, which results in a lower emissivity at z ≲ 5, a larger normalization is required. In addition, a larger β gives less ionizing photons at z ∼ 8, which slows the process of reionization and consequently increases the optical depth. When fesc,*|z = 5 reaches 0, the model becomes quasar dominated, returning a low τe of around a half of the Planck16 measurement (see Fig. 6). We see that when including the contribution from quasars, the model prefers a combination of fesc,*|z = 5 ∼ 6 per cent and β ∼ 0.5. This corresponds to fesc,* ∼ 6.5 per cent at z = 6 with a ratio between the emissivities of quasars and stars of ∼0.12.
5.4 Faintest and brightest quasar contributors
In addition to the possibility that quasars do not have a very high escape fraction (Barkana & Loeb 2001), we note that the very faintest quasars predicted in the model are not observed. The recent detection by G15 only reaches to M1450 ∼ −18, while our fiducial models predict a significant population of faint quasars down to M1450 ∼ −11. Whether those undetected quasars are able to contribute a significant amount of ionizing photons is still unknown. For instance, they might be buried in the dust with a large obscuration fraction (Hopkins et al. 2007). The critical mass above which quasars can contribute ionizing photons, coupled with the previously discussed escape fraction, represents limiting cases of a mass-dependent escape fraction for quasars. In the top panel of Fig. 9, we present the cumulative fraction of ionizing photons as a function of black hole mass (or the corresponding UV magnitude M1450 during the Eddington state as shown in the top axis) assuming fesc,q = 1 from the fiducial model using the Tiamat simulation. We see that quasars fainter than M1450 = −18 contribute approximately 80 per cent of the total emissivity at z ∼ 7 with a decreasing contribution towards lower redshifts (10 per cent at z ∼ 3, the end of the EoR in the QuasarReion model). This suggests that the number of fainter quasars becomes relatively smaller at later times, which can also be observed from the slope of the predicted quasar luminosity function becoming flatter from z ∼ 7–3 (see Fig. 3). However, at redshifts higher than z ∼ 6, the total emissivity from quasars is low. For example, the total emissivity at z ∼ 7 is five times lower than z ∼ 5, suggesting that faint quasars below current observational limits provide only a small contribution to reionization.
In addition, the AGN light curve adopted in this work, which assumes that black holes are either accreting with ε = 1 or stay quiescent (ε = 0) depending on the amount of accretion mass (see Section 2.2.2), has been shown to underestimate the number density of faint AGN at low redshift (Bonoli et al. 2009). For instance, allowing ε to decrease progressively when the accretion disc has been mostly consumed predicts more faint AGN with MB ∼ −16 by a factor of 2 at z ∼ 0.1. However, the impact becomes insignificant at brighter ranges and higher redshifts. Due to the small contribution of ionizing photons from faint quasars, AGN light curves are therefore not expected to have a significant impact on our conclusions regarding reionization.
On the other hand, due to the limited simulation volume the brightest quasars at high redshift (z ≥ 4) in our model only reach M1450 ∼ −23, above which the contribution of ionizing photons is not considered. However, the G15 emissivity accounts for bright quasars up to M1450 = −28, 100 times brighter than the brightest quasar in our model. In order to estimate the emissivity of the missing bright quasars, we integrate the fitting functions provided by G15 with a magnitude interval of −28 < M1450 < −23. We find that the total emissivity at high redshift increases by less than 1 per cent with the inclusion of the ionizing photons from these quasars. Therefore, the conclusion that quasars do not have a significant role during the EoR is not affected by the volume size. However, with the flattening luminosity function at lower redshifts (z < 4, indicated with the vertical dash–dotted line in Fig. 7), bright quasars become more important and their contribution to reionization is not ignorable.27 This results in a lower emissivity of quasars in our model compared to the Haardt & Madau (2012) model at low redshift (see Fig. 7). We note that including these objects will bring forward reionization in the QuasarReion model, but have no impact to the fiducial model.
5.5 Black hole seed mass
Our choice of black hole seed mass, 1000 h−1 M⊙, lies between the light seed (∼102–103 M⊙) from a remnant Pop III star and the massive seed (∼103–105 M⊙) from the direct collapse of a gas cloud at early times (Greene 2012). The massive seed mass is frequently used to initialize massive haloes (≳106–1012 M⊙) in hydrodynamic simulations (Springel, Di Matteo & Hernquist 2005; Vogelsberger et al. 2014; Schaye et al. 2014; Feng et al. 2016), while the ∼103 M⊙ seeds are also often adopted in semi-analytic models (Somerville et al. 2008; Bonoli et al. 2009). We note that this seed mass assumption only affects the black hole mass at early times and in the least massive galaxies. The main conclusions of this work are not significantly affected by this assumption. For example, with exactly the same adopted parameters (see Table 1) but 10 times larger seed mass, the properties such as the black hole mass function, the Magorrian relation and the UV luminosity function are changed by less than 5 per cent in massive galaxies (M* > 109 M⊙). On the other hand, the model predicts a significant number of less massive black holes with masses ∼105 M⊙, which is more than ∼1 order of magnitude larger than the Shankar et al. (2009) sample. However, this has negligible impact on the total instantaneous emissivity, and consequently, the reionization history does not change significantly. Fig. 10 presents the evolution of emissivity, neutral hydrogen fraction and optical depth for the fiducial and QuasarReion models with larger black hole seed masses of 104 h−1 M⊙. Compared to the original models, we see that the quasar emissivity increases with a larger seed mass, while the stellar emissivity decreases due to the stronger feedback from black holes. However, the changes are negligible, resulting in a small perturbation to the reionization history and optical depth.
6 CONCLUSIONS
We have updated the Meraxes semi-analytic model (M16) with a detailed prescription of black hole evolution as part of the Dark-ages Reionization And Galaxy formation Observables from Numerical Simulations (DRAGONS) project to study the role of AGN in reionization and galaxy formation at high redshift. The model is calibrated against the observed stellar mass function (z ∼ 7–0.6), black hole mass function (z ≲ 0.5), quasar luminosity function (z ∼ 6–0.6), ionizing emissivity (z ∼ 5–2) and the Thomson scattering optical depth. The model is in agreement with the observed Magorrian relation at low redshift (z < 0.5) and predicts a decreasing black hole mass towards higher redshifts at a fixed stellar mass. An opening angle of 80 deg, which corresponds to an un-obscured fraction of ∼23.4 per cent, allows the model to reproduce the observed quasar luminosity function across a large redshift range (z ∼ 6–0.6).
Our model suggests that the radiation observed from recently discovered faint AGN at high redshift G15 may include a significant fraction of UV flux from stars. Previous direct estimates of quasar contributions to reionization based on these observations (Madau & Haardt 2015; Mitra et al. 2015) therefore result in an overestimate of the emissivity of quasars by a factor of 3 at z ∼ 6.
When we include the contribution of AGN to reionization, we find that quasars do not dominate the ionizing photon budget at z > 6. In a quasar-only reionization model, where the escape fractions of ionizing photons are 1 and 0 for quasars and stars, respectively, we find that reionization happens very late, z ∼ 3, with a Thomson scattering optical depth of only half of the Planck16 measurement (Planck16). However, at low redshift, quasars are able to provide a large number of ionizing photons. With quasars contributing all of their ionizing photons (fesc,q = 1), our model prefers a redshift-dependent escape fraction for stars, having the form of |$f_\mathrm{esc,*}\left(z\right)=\min \left[0.06\times \left(\frac{1+z}{6}\right)^{0.5}, 1\right]$|. This corresponds to quasars contributing 10 per cent of the total ionizing photons at z ∼ 6.
Acknowledgements
We would like to thank the anonymous referees for providing helpful suggestions that improves the paper substantially. This research was supported by the Victorian Life Sciences Computation Initiative (VLSCI), grant no. UOM0005, on its Peak Computing Facility hosted at the University of Melbourne, an initiative of the Victorian Government, Australia. Part of this work was performed on the gSTAR national facility at Swinburne University of Technology. gSTAR is funded by Swinburne and the Australian Governments Education Investment Fund. This research programme is funded by the Australian Research Council through the ARC Laureate Fellowship FL110100072 awarded to JSBW. This work was supported by the Flagship Allocation Scheme of the NCI National Facility at the ANU, generous allocations of time through the iVEC Partner Share and Australian Supercomputer Time Allocation Committee. AM acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant no. 638809 – AIDA).
Note that in this work, haloes are defined as the substructures of fof groups. Central haloes (the most massive halo in a fof group) and satellites co-evolve in Meraxes and both contribute to reionization.
A central halo consist of the majority of its fof particles and dominates the gravitational potential of the entire system. In the code, the infall gas is added into the central halo based on the baryon fraction of its fof group while the satellites do not get any fresh gas (see more details in Mutch et al. 2016).
Meraxes assumes the Salpeter IMF with a mass range of 0.1–120 M⊙. We use the lifetime-mass relation of stars (Portinari, Chiosi & Bressan 1997) to calculate the fraction (ηSNII) of stellar mass in a single stellar population that have reached the supernova stage after a certain period of time. For example, assuming stars more massive than 8 M⊙ will reach the type-II supernova stage after ∼40 Myr, ηSNII = 7.432 × 10−3.
Dual or multiple AGN are not considered.
It is generally true for the radio mode. However, the quasar mode introduced in Section 2.2.2 can induce significant increases in black hole mass during a single time step, leading to a non-negligible change of the accretion rate.
We do not limit the heating radius (see Croton et al. 2016) to only move outwards.
Ignoring the loss due to gravitational wave emission.
Calibration using the MCMC technique is ongoing. In this work, the calibration is performed by hand.
In order to avoid crowded presentations, only the Tiamat-125-HR uncertainty is shown when the Tiamat and Tiamat-125-HR results are both present in one plot.
A black hole mass – galaxy property scaling relation usually accounts for the bulge property. However, the majority of systems with black holes are expected to be bulge dominated. Therefore, using the total stellar mass as a proxy does not lead to a significant bias in Fig. 1. We also note that recent studies suggest that while the Mbh–M* relation is significantly biased (see Section 5), it is the bulge velocity dispersion that connects supermassive black holes and host galaxies. However, interpreting the scaling relation is beyond the scope of this work. We leave a detailed analysis of the black hole–galaxy scaling relation to future work when bulge properties are included (e.g. Tonini et al. 2016).
We exclude the recently identified haloes in order to minimize the effect of black hole seeding. However, this does not have a significant impact in Fig 1.
Radio mode accretion and merger-driven starburst are relatively less important to the growth.
We note that the Tiamat halo merger trees have been improved (Poole et al. in preparation) since M16. However, the impact on galaxy formation is trivial.
Ignoring the impact on the ionizing photon escape fraction from AGN feedback.
An opening angle of θ = 80 deg corresponds to a solid angle of |$\Omega = 2\pi \,\,(1-\cos \frac{\theta }{2})\sim 1.47$|. Considering a symmetric radiation from both sides of the accretion disc, the un-obscured fraction is |$f_\mathrm{obs}\equiv \frac{2\Omega }{4\pi }\sim 0.234$|.
The opening angle is used to illustrate the observable fraction, assuming an orientation model. We note that the observable fraction can also be interpreted by the line-of-sight absorption column density when the evolutionary model is assumed.
Also with a constant to model the sky background.
Dust attenuation is not considered because these relatively faint galaxies have little dust in our models (Liu et al. 2016).
The Tiamat-125-HR halo merger trees cannot resolve small objects; therefore, we cannot consider ionizing photons from faint galaxies or quasars. We only show the result calculated using the Tiamat trees in the following sections.
The mass of black hole seed is subtracted because they do not produce any ionizing photons. Reionization from the progenitor of black hole seeds will be considered in the future when Population III stars are implemented.
One may also define fesc,qfobs as the quasars escape fraction.
The stellar mass recycled to the ISM through supernova also contributes ionizing photons in the H ii bubble.
The systematic error due to recombination radiation is ignored.
If it is not a selection bias, it is possible that not every galaxy hosts a central massive black hole.
We fit the predicted quasar luminosity function at z = 2 using a single power law at M1450 > −24 and estimate the number of missing photons with a magnitude interval of −24 < M1450 < −22.5. We find the emissivity of quasars can be increased by a factor of 2 with the inclusion of the missing bright quasars.
REFERENCES
APPENDIX A: CALCULATING THE MEAN NUMBER OF IONIZING PHOTONS PRODUCED PER BLACK HOLE
During one time step, for a black hole with a given initial mass of MBH, its bolometric luminosity at Eddington rate can be calculated through the right hand of equation (16). Since the bolometric correction (Hopkins et al. 2007) adopted in this work is dependent on the bolometric luminosity (see equation 25), the UV magnitude, M1450 of the quasar changes during its accretion, so does the emissivity. In our model, because the accretion mass is always smaller than the black hole mass (ΔMBH < MBH, see Fig. A1), for the sake of calculation speed, we estimate the mean number of ionizing photons produced per black hole, Nγ, q with the bolometric luminosity at the beginning of accretion. We calculate Nγ, q as follows:
We calculate the UV magnitude, M1450 using equa-tions (22)–(27).
- We calculate the UV flux with M1450 in units of erg s−1Hz−1 through(A1)\begin{equation} F_{1450} = 10^{\left(M_{1450}-48.6\right)/-2.5}\times 4\pi \left(\frac{10\,\,\mathrm{pc}}{1\mathrm{\,\,cm}}\right)^2. \end{equation}
- We calculate the flux at Lyman limit following G15where αq,optical = 0.44 and αq = 1.57 correspond to a double power-law AGN SED.(A2)\begin{equation} F_{912} = F_{1450}\left(\frac{1200}{1450}\right)^{\alpha _{q,\mathrm{optical}}}\left(\frac{912}{1200}\right)^{\alpha _\mathrm{q}}, \end{equation}
- We calculate the instantaneous emissivity by(A3)\begin{equation} \dot{N}_\mathrm{ion} \equiv \int _{\nu _{912}}^{\infty } F_{912}\left(\frac{\nu }{\nu _{912}}\right)^{-\alpha _\mathrm{q}} \frac{\text{d}\nu }{h\nu } = \frac{F_{912}}{h\alpha _\mathrm{q}}. \end{equation}
- The duration of accreting mass, ΔMBH can be calculated throughTherefore, the total number of ionizing photons emitted is |$\dot{N}_\mathrm{ion}t_\mathrm{acc}$| and the mean number of ionizing photons produced per black hole is(A4)\begin{equation} t_\mathrm{acc} = \ln \left(\frac{\Delta {M}_\mathrm{BH}}{M_\mathrm{BH}} +1\right)\times \frac{\eta t_\mathrm{Edd}}{\epsilon }. \end{equation}(A5)\begin{equation} N_{\gamma ,q} = \frac{\int _{0}^{t_\mathrm{acc}}\dot{N}_\mathrm{ion}\left(t\right)\mathrm{d}t}{\left(1-\eta \right)\Delta {M}_\mathrm{BH}/m_\mathrm{p}} \approx \frac{\dot{N}_\mathrm{ion}|_{t=t_\mathrm{acc}/2}t_\mathrm{acc}}{\left(1-\eta \right)\Delta {M}_\mathrm{BH}/m_\mathrm{p}}, \end{equation}
where the last step adopts the instantaneous emissivity at the middle of accretion for the sake of computational speed.
We note that during the accretion, with an exponential increase of black hole mass, the AGN bolometric luminosity, Lbol increases exponentially (equation 22). Since the B-band bolometric correction, kB, decreases with increasing luminosity following a double power law (equation 25), |$\dot{N}_\mathrm{ion}$| is a convex function of time. Therefore, the approximation in equation (A5) underestimates the number of ionizing photons produced by black hole. In order to test whether this has a significant impact on our conclusion, we rerun the QuasarReion model assuming a constant bolometric correction with |$k_\mathrm{B}\left(t\right) \approx k_\mathrm{B}|_{t=t_\mathrm{acc}}$|. Eliminating the complex dependence of time from kB, Nγ, q can be analytically calculated by integrating the AGN light curve. However, we note that since |$k_\mathrm{B}\left(t\right) \le k_\mathrm{B}|_{t=t_\mathrm{acc}}$|, this approximation overestimates Nγ, q.
Fig. A2 presents the evolution of emissivity, neutral hydrogen fraction and optical depth for different QuasarReion models assuming constant |$\dot{N}_\mathrm{ion}$| (QuasarReion) and kB (QuasarReion_kB), respectively. Since the time interval between two snapshots is much smaller than the Eddington accretion time-scale (tEdd ∼ 450 Myr), the black hole mass increment is still within the linear regime and therefore we see that the impact from the calculation of Nγ, q is not significant.