Dark-ages Reionization and Galaxy Formation Simulation - X. The small contribution of quasars to reionization

TL;DR

This study uses an updated semi-analytic model to assess the role of quasars in cosmic reionization, finding their contribution to be minimal compared to galaxies, especially at high redshift.

Contribution

The paper introduces a refined model that incorporates recent faint AGN data and dust obscuration effects to evaluate quasar impact on reionization.

Findings

Quasars contribute minimally to reionization at high redshift.

Faint AGN luminosity functions overestimate ionizing photons from quasars.

Stellar light from faint AGN hosts significantly affects UV flux measurements.

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\sim0.6-7$, the black hole mass function at $z\lesssim0.5$, the global ionizing emissivity at $z\sim2-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 ${\sim}23.4$ per cent due to obscuration by dust, the model is able to reproduce the observed quasar luminosity function at $z\sim0.6-6$. The stellar light from galaxies hosting faint 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\sim4$ 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\sim6$. We conclude that quasars are unlikely to make a significant contribution to reionization.

Figures

Tables

Table 1. Table 1: A list of relevant parameters in the model with description and adopted value. The redshift varying f esc , ∗ subscript 𝑓 esc f_{\mathrm{esc,*}} model (Mutch et al . , 2016 , hereafter M16 ) is shown for comparison. The M16BH model includes black hole feedback compared to the M16 model but without reionization from quasars. We refer the interested reader to Mutch et al . ( 2016 ) .

Parameter	Section	Equation	Description	Fiducial	M16BH	M16
${\alpha }_{\mathrm{sf}}$	2.1	6	Star formation efficiency	0.08	0.03	0.03
${\alpha }_{\mathrm{energy}}$	2.1	7	Energy coupling efficiency normalization	1.0	0.5	0.5
${\alpha }_{\mathrm{mass}}$	2.1	9	Mass loading normalization	15.0	9.0	9.0
\hdashline$\eta$	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	-
$k_h$	2.2.1	14	Black hole growth efficiency for the radio mode	0.3	0.1	-
${\kappa }_r$	2.2.2	16	Black hole feedback efficiency for the radio mode	1.0	1.0	-
$k_c$	2.2.2	19	Black hole growth efficiency for the quasar mode	0.05	0.05	-
${\kappa }_q$	2.2.2	-	Black hole feedback efficiency for the quasar mode	0.0005	0.0005	-
$N_{\gamma ,\ast }$	4.1	28	Mean number of ionizing photons produced per stellar baryon	4,000	4,000	4,000
$f_{\mathrm{esc},q}$	4.1	28	Ionizing photon escape fraction for quasars	1.0a	0.0	-
$f_{\mathrm{esc},\ast }$	4.1	28	Ionizing photon escape fraction for the stellar component	$f_{\mathrm{esc},\ast ,z}$b,c	$f_{\mathrm{esc},\ast ,\mathrm{M16}}$d	$f_{\mathrm{esc},\ast ,\mathrm{M16}}$

Equations

$$\Delta m_{\mathrm{hot}}=\max \left[0,{\chi }_{\mathrm{r}}{f}_{\mathrm{b}}{M}_{\mathrm{vir}}-\left(m_{\star }+m_{\mathrm{cold}}+m_{\mathrm{hot}}+m_{\mathrm{eject}}\right)\right],$$

$$t_{\mathrm{cool}}\left(r\right)=\frac{3\bar{\mu }{m}_{\mathrm{p}}k{T}_{\mathrm{hot}}}{2{\rho }_{\mathrm{hot}}\left(r\right)\Lambda \left(T_{\mathrm{hot}},Z_{\mathrm{hot}}\right)},$$

$$t_{\mathrm{cool}}\left(r\right)=\frac{6\pi \bar{\mu }{m}_{\mathrm{p}}k{T}_{\mathrm{vir}}{R}_{\mathrm{vir}}{r}^2}{m_{\mathrm{hot}}\Lambda \left(T_{\mathrm{vir}},Z_{\mathrm{hot}}\right)},$$

$$r_{\mathrm{cool}}=\sqrt{\frac{m_{\mathrm{hot}}\Lambda \left(T_{\mathrm{vir}},Z_{\mathrm{hot}}\right)}{6\pi \bar{\mu }{m}_{\mathrm{p}}k{T}_{\mathrm{vir}}{V}_{\mathrm{vir}}}},$$

$$m_{\mathrm{cool}}\equiv m_{\mathrm{hot}}\times \min \left[1,\min \left(1,\frac{r_{\mathrm{cool}}}{R_{\mathrm{vir}}}\right)\times \frac{\Delta t}{t_{\mathrm{cool}}}\right],$$

$$\Delta m_{\star }=\min \left[\frac{{\alpha }_{\mathrm{sf}}\max \left(0,m_{\mathrm{cold}}-m_{\mathrm{crit}}\right)\Delta t}{t_{\mathrm{dyn},\mathrm{disc}}},m_{\mathrm{cold}}\right],$$

$$\begin{array}{rl}{e}_{\mathrm{total}}=& \min \left\{{\alpha }_{\mathrm{energy}}\left[0.5+{\left(\frac{V_{\max }}{70\mathrm{km}{\mathrm{s}}^{-1}}\right)}^{-2.0}\right],1\right\}\\ & \times 10^{51}\mathrm{erg}\times {\eta }_{\mathrm{SN}\mathrm{II}}\Delta m_{\star },\end{array}$$

$$E_{\mathrm{total}}={\sum }_{i=j-4}^{i=j}{e}_{\mathrm{total},i}\left(\Delta m_{\star ,i},V_{\max ,\mathrm{i}},{\eta }_{\mathrm{SNII},\mathrm{i}}\right).$$

$$\begin{array}{r}{E}_{\mathrm{reheat}}=\frac{1}{2}{\alpha }_{\mathrm{mass}}\Delta m_{\star }{V}_{\mathrm{vir}}^2,\end{array}$$

$$m_{\mathrm{reheat}}=\min \left[m_{\mathrm{cold}},\frac{\min \left(E_{\mathrm{total}},E_{\mathrm{reheat}}\right)}{0.5V_{\mathrm{vir}}^2}\right].$$

$$\Delta m_{\mathrm{eject}}=\min \left[m_{\mathrm{hot}},\frac{\max \left(0,E_{\mathrm{total}}-0.5m_{\mathrm{reheat}}{V}_{\mathrm{vir}}^2\right)}{0.5V_{\mathrm{vir}}^2}\right],$$

$${\dot{m}}_{\mathrm{Bondi}}=\frac{2.5\pi G^2{m}_{\mathrm{bh}}^2{\rho }_{\mathrm{hot}}}{c_{\mathrm{s}}^3},$$

$$\frac{{\rho }_{\mathrm{hot}}}{c_{\mathrm{s}}^3}=\frac{3\mu m_{\mathrm{p}}kT}{4G\Lambda m_{\mathrm{bh}}}.$$

$$\Delta m_{\mathrm{bh},\mathrm{hot}}=\min \left(m_{\mathrm{hot}},m_{\mathrm{Edd}},k_{\mathrm{h}}{\dot{m}}_{\mathrm{Bondi}}\Delta t\right),$$

$$\Delta m_{\mathrm{bh}}=(1-\eta )\Delta m_{\mathrm{bh},\mathrm{hot}}.$$

$$\eta {\dot{m}}_{\mathrm{bh}}(t)c^2=ϵ\frac{4\pi G{m}_{\mathrm{bh}}(t)m_{\mathrm{p}}c}{{\sigma }_{\mathrm{T}}},$$

$$m_{\mathrm{Edd}}=m_{\mathrm{bh}}\left[\exp \left(\frac{ϵ\Delta t}{\eta t_{\mathrm{Edd}}}\right)-1\right],$$

$$m_{\mathrm{heat}}=\frac{{\kappa }_{\mathrm{r}}\eta \Delta m_{\mathrm{bh},\mathrm{hot}}{c}^2}{0.5V_{\mathrm{vir}}^2},$$

$$\Delta m_{\mathrm{bh},\max }=\min \left[m_{\mathrm{cold}},\frac{k_{\mathrm{c}}\gamma m_{\mathrm{cold}}}{{\left(1+\frac{280\mathrm{km}{\mathrm{s}}^{-1}}{V_{\mathrm{vir}}}\right)}^2}\right]+\Delta m_{\mathrm{bh},\max }^{\prime },$$

$$\Delta m_{\mathrm{bh},\mathrm{cold}}=\min \left(m_{\mathrm{Edd}},\Delta m_{\mathrm{bh},\max }\right),$$

$$t_{\mathrm{acc}}=\ln \left(\frac{\Delta m_{\mathrm{bh},\mathrm{cold}}}{m_{\mathrm{bh}}}+1\right)\times \frac{\eta t_{\mathrm{Edd}}}{ϵ}.$$

$$L_{\mathrm{bol}}\equiv ϵm_{\mathrm{bh}}{|}_{t=t_{\mathrm{obs}}}\frac{c^2}{t_{\mathrm{Edd}}}=ϵm_{\mathrm{bh}}\exp \left(\frac{ϵt_{\mathrm{obs}}}{\eta t_{\mathrm{Edd}}}\right)\times \frac{c^2}{t_{\mathrm{Edd}}},$$

$$M_{\mathrm{bol}}=4.74-2.5{\log }_{10}\left(\frac{L_{\mathrm{bol}}}{{\mathrm{L}}_{\odot }}\right).$$

$$M_{\mathrm{bol}}-M_{\mathrm{B}}=-2.5{\log }_{10}{k}_{\mathrm{B}},$$

$$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}.$$

$$M_{\mathrm{AB},\mathrm{B}}-M_{\mathrm{B}}=-0.09$$

$$M_{1450}=M_{\mathrm{AB},\mathrm{B}}-2.5{\log }_{10}{\left(\frac{1450\mathring{\mathrm{A}}}{4344\mathring{\mathrm{A}}}\right)}^{{\alpha }_{q,\mathrm{optical}}}=M_{\mathrm{AB},\mathrm{B}}+\mathrm{0.524.}$$

$$N_{\ast }{N}_{\gamma ,\ast }{f}_{\mathrm{esc},\ast }+N_{\mathrm{q}}{N}_{\gamma ,q}{f}_{\mathrm{obs}}{f}_{\mathrm{esc},\mathrm{q}}\ge \left(1+{\bar{N}}_{\mathrm{rec}}\right)N_{{\mathrm{H}}_{\mathrm{I}}},$$

$$\frac{{\xi }_{\ast }{m}_{\ast }+{\xi }_{\mathrm{q}}{m}_{\mathrm{bh}}}{\frac{4}{3}\pi R^3{\Omega }_{\mathrm{m}}{\rho }_{\mathrm{c}}\left(z\right)}\ge 1,$$

$${\xi }_{\ast }=\frac{N_{\gamma ,\ast }{f}_{\mathrm{esc},\ast }}{f_{\mathrm{b}}(1-0.75Y_{\mathrm{He}})}\mathrm{and}{\xi }_{\mathrm{q}}=\frac{N_{\gamma ,q}{f}_{\mathrm{obs}}{f}_{\mathrm{esc},\mathrm{q}}}{f_{\mathrm{b}}(1-0.75Y_{\mathrm{He}})},$$

Full text

††pagerange:

Dark-ages Reionization and Galaxy Formation Simulation - X.

The small contribution of quasars to reionization–12††pubyear: 2017

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 ($f_{\mathrm{esc,*}}{\gtrsim}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 et al., 2016). However, the value of $f_{\mathrm{esc,*}}$ 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$\alpha$ to H$\beta$ 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$\alpha$-selected galaxies. Moreover, some theoretical works also suggest a low $f_{\mathrm{esc,*}}$ (Gnedin et al., 2007; 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 $f_{\mathrm{esc,*}}$ with redshift (Haardt & Madau, 2012; Khaire et al., 2016; Price et al., 2016; Sharma et al., 2016) and with decreasing mass (Paardekooper et al., 2013; Kimm & Cen, 2014; Wise et al., 2014). This is supported by identifying local analogues of high-redshift galaxies and extrapolating the observed $f_{\mathrm{esc,*}}$ using indicators such as the [Oiii]/[Oii] 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 ($f_{\mathrm{esc,q}}\sim 1$, Barkana & Loeb 2001), quasars (AGN) are potential contributors to reionization (Volonteri & Gnedin, 2009; Fontanot et al., 2014; Madau & Haardt, 2015; Mitra et al., 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 et al., 2017) and whether the escape fraction of high-redshift low-luminosity AGN is of order of unity (Cristiani et al., 2016; Micheva et al., 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 Meraxes111http://dragons.ph.unimelb.edu.au 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 ($\Omega_{\mathrm{m}},\Omega_{\mathrm{b}},\Omega_{\mathrm{\Lambda}},h,\sigma_{8},n_{\mathrm{s}}$ = 0.308, 0.0484, 0.692, 0.678, 0.815, 0.968; Planck Collaboration et al. 2016a).

2 Modelling Black hole growth

Built on halo222Note 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. 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 ${\sim}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 ${\sim}2.64\times 10^{6}h^{-1}\mathrm{M}_{\odot}$ and the box size is 67.8$h^{-1}$Mpc. 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\times 10^{9}h^{-1}\mathrm{M}_{\odot}$ in a bigger simulation volume, with side lengths equal to $125h^{-1}\mathrm{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 the Meraxes semi-analytic model, cooling and star formation are assumed to be negligible in haloes below the atomic cooling mass threshold, ${\sim}10^{8}\mathrm{M}_{\odot}$. Thus, once a halo grows larger than the atomic cooling limit it is designated as a galaxy, which has three baryonic components: gas, stars and a central black hole. During one time step, $\Delta t$, additional gas falls into the hot gas component of a galaxy from the IGM when the mass fraction of baryons in the halo is lower than the cosmic mean value, $f_{\mathrm{b}}{=}\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}$:

[TABLE]

where $m_{\mathrm{\star}},\ m_{\mathrm{cold}},\ m_{\mathrm{hot}}$ and $m_{\mathrm{eject}}$ are the masses of the stellar component, cold gas, hot gas and ejected gas, respectively, and $M_{\mathrm{vir}}$ is the virial mass of the host halo in which the galaxy forms.333A 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). $\chi_{\mathrm{r}}$ is a baryon fraction modifier to take account of reionization feedback and will be introduced later in Section 4.1.

Some of the hot gas, $m_{\mathrm{cool}}$, cools and collapses on to the cold disc. Assuming the cooling process is in quasi-static thermal equilibrium, one can calculate the cooling time by

[TABLE]

where $\mu=0.59$, $m_{\mathrm{p}}$ and $k$ are the mean molecular weight for fully ionized gas, the mass of a proton and the Boltzmann constant, respectively. $T_{\mathrm{hot}}$ and $\rho_{\mathrm{hot}}$ are the temperature and density of the hot gas. $\Lambda$ is the cooling function (Sutherland & Dopita, 1993) determined by the temperature and metallicity, $Z_{\mathrm{hot}}$. Assuming the hot gas shares the same temperature as the host halo ($T_{\mathrm{hot}}=T_{\mathrm{vir}}$) due to shock heating, and follows a singular isothermal sphere density profile, equation (2) becomes

[TABLE]

where $R_{\mathrm{vir}}$ is the virial radius. Following Croton et al. (2006), we calculate the cooling radius, $r_{\mathrm{cool}}$, at which the cooling time is equal to the halo dynamical time, through

[TABLE]

where $V_{\mathrm{vir}}$ is the virial velocity. Cooling is sufficient within $r_{\mathrm{cool}}$ and we estimate the cooling mass by

[TABLE]

which is removed from the hot gas reservoir and redistributed into the cold gas disc, $\Delta m_{\mathrm{hot}}=-\Delta m_{\mathrm{cold}}{=}-m_{\mathrm{cool}}$. From equation (5), we see that depending on the ratio of $r_{\mathrm{cool}}$ to $R_{\mathrm{vir}}$, cooling is separated into two regimes: static hot halo ($r_{\mathrm{cool}}>R_{\mathrm{vir}}$) and rapid cooling ($r_{\mathrm{cool}}<R_{\mathrm{vir}}$, see more in Croton et al. 2006). In the case of a static hot halo, cooling is in thermal equilibrium and the cooling rate is determined by the continuity equation. On the other hand, if the cooling time is shorter than the dynamical time this leads to a rapid cooling process, with a free-fall of all the hot gas on to the cold gas disc. We note that $m_{\mathrm{cool}}$ is set to be zero when the mass of the host halo drops below the atomic cooling due to stripping.

When the galaxy collects enough cold gas, $m_{\mathrm{cold}}{>}m_{\mathrm{crit}}$, based on Kennicutt Jr (1998) and Kauffmann (1996), it forms new stars, $\Delta m_{\star}$, through a burst

[TABLE]

where $t_{\mathrm{dyn,disc}}=1.5\sqrt{2}\lambda R_{\mathrm{vir}}/V_{\mathrm{vir}}$ is the dynamical time of the cold gas disc ($\lambda$ is the spin parameter of the host halo) and $\alpha_{\mathrm{sf}}$ is a free parameter corresponding to the star formation efficiency. This mass is removed from the cold gas reservoir, $\Delta m_{\mathrm{cold}}=-\Delta m_{\star}$.

The new stellar mass is assumed to form following a Salpeter (1955) initial mass function (IMF). Ultimately some of these stars, $\eta_{\mathrm{SN{II}}}\Delta m_{\star}$, recycle their mass back to the interstellar medium (ISM) through type-II supernovae explosion.444Meraxes assumes the Salpeter IMF with a mass range of $0.1-120\mathrm{M}_{\odot}$. We use the lifetime-mass relation of stars (Portinari et al., 1997) to calculate the fraction ($\eta_{\mathrm{SN{II}}}$) 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\mathrm{M}_{\odot}$ will reach the type-II supernova stage after ${\sim}40$ Myr, $\eta_{\mathrm{SN{II}}}=7.432\times 10^{-3}$. The energy produced by these supernovae heats the ISM and converts some of the cold gas to the hot phase or, if there is sufficient energy, even ejects a fraction of hot gas from the galaxy. Assuming that the efficiency for supernovae energy coupling with the ISM scales with mass and is in proportion to a free parameter (Guo et al., 2011), $\alpha_{\mathrm{energy}}$, then the total energy released by supernovae that couples to the ISM is

[TABLE]

where $V_{\mathrm{max}}$ is the maximum circular velocity of the host halo. Since a massive star (${\gtrsim}8\mathrm{M}_{\odot}$) takes ${\sim}40$ million years, or 4 snapshots, before reaching its type-II supernova stage, Meraxes accounts for supernovae not only from the current snapshot, $j$, but also from the stars formed in the previous 4 snapshots. Therefore, the total energy released during one snapshot is

[TABLE]

Since the hot gas shares the same virial temperature as the host halo, assuming the mass loading factor for reheating cold gas is a free parameter, $\alpha_{\mathrm{mass}}$, the energy utilized in gas heating can be calculated by

[TABLE]

where $V_{\mathrm{vir}}$ is the virial velocity. Depending on the available energy, $E_{\mathrm{total}}$ and the required energy for re-heating, $E_{\mathrm{reheat}}$, the reheated mass, $m_{\mathrm{reheat}}$, is

[TABLE]

This mass is removed from the cold gas reservoir and redistributed in the hot gas component, $\Delta m_{\mathrm{cold}}=-\Delta m_{\mathrm{hot}}=-m_{\mathrm{reheat}}$. If there is still some energy left after reheating, the supernovae feedback will further remove hot gas from the galaxy, adding it to the ejected component

[TABLE]

which is removed from the hot gas reservoir, $\Delta m_{\mathrm{hot}}=-\Delta m_{\mathrm{eject}}$.

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 et al., 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 hole555Dual or multiple AGN are not considered. of mass $1000\ h^{-1}\mathrm{M}_{\odot}$. 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

Whenever there is a static hot gas reservoir, $m_{\mathrm{hot}}$, around the galaxy, some of it will cool, $m_{\mathrm{cool}}$, and form a cold gas disc, $m_{\mathrm{cold}}$ (see the previous section), while some will be directly accreted by the central black hole. We adopt the Bondi-Hoyle accretion model proposed in Croton et al. (2016) to describe this hot gas accretion. The Bondi-Hoyle accretion rate is given by

[TABLE]

where $G,m_{\mathrm{bh}},c_{\mathrm{s}}$ and $\rho_{\mathrm{hot}}$ are the gravitational constant, the black hole mass, the speed of sound and the density of the hot gas reservoir, respectively. We define the Bondi radius as $r_{\mathrm{Bondi}}=\dfrac{2Gm_{\mathrm{bh}}}{c_{\mathrm{s}}^{2}}$, and by equating the sound crossing travel time from the centre to the Bondi radius with the local cooling time (see equation 2), we obtain

[TABLE]

Then, assuming the accretion rate does not change during one time step,666It 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. $\Delta t$, the accretion mass is

[TABLE]

where $k_{\mathrm{h}}$ is a free parameter, used to adjust the efficiency of black hole growth in the radio mode since black holes may not be accreting at the full Bondi-Hoyle accretion rate. $m_{\mathrm{Edd}}$ is the Eddington limit, which will be introduced later. This mass is removed from the hot gas reservoir, $\Delta m_{\mathrm{hot}}=-\Delta{m}_{\mathrm{bh,hot}}$. Assuming a fraction, $\eta$, of the accreted mass is radiated, the black hole only grows by

[TABLE]

Moreover, the radiation acting outward is limited by the Eddington luminosity

[TABLE]

where $\epsilon$ is the Eddington ratio ($\epsilon{=}1$ in this work), $m_{\mathrm{bh}}(t)$ and $\dot{m}_{\mathrm{bh}}(t)$ are the black hole mass and growth rate at time $t$, respectively, $\sigma_{\mathrm{T}}$ is the Thomson cross-section area and $c$ is the speed of light. By integrating equation 16 through one time step, this provides the second limitation in equation (14)

[TABLE]

where $t_{\mathrm{Edd}}\equiv\dfrac{\sigma_{\mathrm{T}}c}{4\pi Gm_{\mathrm{p}}}\approx 450\mathrm{Myr}$ is the Eddington accretion time-scale and $m_{\mathrm{bh}}$ is the black hole mass at the beginning of this time step.

Assuming adiabatic heating and that a fraction of the radiated energy, $\kappa_{\mathrm{r}}$, is coupled to the surrounding gas and therefore can contribute to feedback, then the heated mass due to black hole feedback can be calculated through

[TABLE]

which is subtracted from the cooling flow, $\Delta m_{\mathrm{cool}}{=}{-}m_{\mathrm{heat}}$. However, if the heating from black hole feedback is too strong, $m_{\mathrm{heat}}{>}m_{\mathrm{cool}}$ (see equation 5), it will significantly suppress the cooling flow, which will consequently restrain the black hole accretion of hot gas. This suppression is referred to as radio mode feedback. In this case, following777We do not limit the heating radius (see Croton et al. 2016) to only move outwards. Croton et al. (2016), $\Delta m_{\mathrm{bh,hot}}$ is rescaled to be the amount of mass within the cooling radius, $\Delta m_{\mathrm{bh,hot}}=\Delta m_{\mathrm{bh,hot}}\times\dfrac{m_{\mathrm{cool}}}{m_{\mathrm{heat}}}$ and the heated mass consequently shrinks to be the cooling mass, $m_{\mathrm{heat}}{=}m_{\mathrm{cool}}$, resulting in a complete quenching of cooling, $\Delta m_{\mathrm{cool}}=-m_{\mathrm{heat}}\rightarrow m_{\mathrm{cool}}=0$. Radio mode feedback limits gas condensation in cluster cooling flows and regulates star formation in massive galaxies (Croton et al., 2006). We note that this does not have a significant impact on the results from Tiamat at $z\geq 2$ due to the limited number of massive objects in the box. However, black hole feedback is important in more massive objects at lower redshifts, which can be observed using the Tiamat-125-HR halo merger trees.

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 ${\sim}0.1$ per cent of the mass from the quasar mode (see Appendix A) 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

Many binary AGN have been detected in merging galaxies (Shields et al., 2012; Comerford et al., 2013; Comerford & Greene, 2014; Comerford et al., 2015; Müller-Sánchez et al., 2015). This suggests that galaxy mergers might trigger AGN activity, which is also supported by hydrodynamic simulations of galaxy mergers (Capelo et al., 2015; Volonteri et al., 2015a, b; Steinborn et al., 2016). Following Croton et al. (2016), when a merger occurs between two galaxies, the central black holes coalesce and their masses combine.888Ignoring the loss due to gravitational wave emission. As mergers drive strong gas inflows towards the central region (Capelo et al., 2015), cold gas is funnelled to the central black hole of the resulting merged galaxy, significantly increasing its mass. The amount of accretion can be estimated by (Bonoli et al., 2009; Croton et al., 2016)

[TABLE]

where $m_{\mathrm{cold}}$ is the amount of mass available in the cold gas disc, $k_{\mathrm{c}}$ is a free parameter used to modulate the strength of black hole accretion and $\gamma{\leq}1$ is the mass ratio between the two merging galaxies, respectively. The term $\Delta{m}_{\mathrm{bh,max}}^{\prime}$ corresponds to the accretion mass left from the quasar mode in the previous snapshot, which will be introduced later.

Unlike the radio mode, black holes grow dramatically during the quasar mode. The AGN activity lifetime (a few $10^{7}$ yr up to a Gyr Fiore et al. 2012) is much longer than the 11Myr time step at $z>5$. In some cases, this leaves the central black hole insufficient time to consume all of the newly accreted gas, $\Delta{m}_{\mathrm{bh,max}}$, at the Eddington limit. Therefore, the mass actually accreted by the black hole is

[TABLE]

which is removed from the cold gas reservoir, $\Delta m_{\mathrm{cold}}=-\Delta{m}_{\mathrm{bh,cold}}$. In our model, during the quasar mode black holes are assumed to either accrete and radiate at the Eddington rate or stay quiescent if the accretion mass is not sufficient. Therefore, depending on the total available mass brought in, $\Delta{m}_{\mathrm{bh,max}}$, and the Eddington limit, $m_{\mathrm{Edd}}$, there are two possible scenarios when the central black hole is undergoing a merger:

$\Delta{m}_{\mathrm{bh,max}}<m_{\mathrm{Edd}}$. 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 through

[TABLE]

When the quasar is observed at a random time $t_{\mathrm{obs}}$, the bolometric luminosity is

[TABLE]

when $t_{\mathrm{obs}}\leqslant t_{\mathrm{acc}}$, otherwise $L_{\mathrm{bol}}=0$. 2. 2.

$\Delta{m}_{\mathrm{bh,max}}\geqslant m_{\mathrm{Edd}}$. In this case, the merger event delivers sufficient cold gas into the accretion disc. In the case of $\Delta{m}_{\mathrm{bh,max}}>m_{\mathrm{Edd}}$, 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 $t_{\mathrm{obs}}$, 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 ($\epsilon=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 $\kappa_{\mathrm{q}}\eta\Delta m_{\mathrm{bh,cold}}c^{2}$, where $\kappa_{\mathrm{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 $\mathrm{\AA}$ band luminosities of quasars are calculated as follows:

We calculate the bolometric magnitude through

[TABLE] 2. 2.

We calculate the B-band magnitude in the Vega magnitude system using the bolometric correction proposed by Hopkins et al. (2007),

[TABLE]

where

[TABLE] 3. 3.

We convert the B-band magnitude from the Vega system to the AB system following Glikman et al. (2010)

[TABLE] 4. 4.

We extrapolate the B-band magnitude of which the effective wavelength is 4344 $\mathrm{\AA}$ (Blanton & Roweis, 2007) to the 1450$\mathrm{\AA}$ magnitude, assuming the quasar continuum between 1450 and 4344$\mathrm{\AA}$ has a power-law slope of $\alpha_{q,\mathrm{optical}}=0.44$ relative to wavelength (Schirber & Bullock, 2003). Thus

[TABLE]

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 Mutch et al. (2016). All other parameters remain the same as Mutch et al. (2016). In this work, we constrain999Calibration using the MCMC technique is ongoing. In this work, the calibration is performed by hand. 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., 2015; Huertas-Company et al., 2016; Stefanon et al., 2016; Davidzon et al., 2017) between $z{\sim}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; 2. 2.

the observed black hole mass function (Graham et al., 2007; Shankar et al., 2009; Vika et al., 2009; Davis et al., 2014; Mutlu-Pakdil et al., 2016) and Magorrian relation (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) at low redshift ($z\lesssim 0.5$); 3. 3.

the latest integrated free electron Thomson scattering optical depth measurement (Planck Collaboration et al., 2016b); 4. 4.

the predicted global ionizing emissivity from Ly$\alpha$ opacities (Becker & Bolton, 2013).

3.1 Black hole properties

The black hole mass functions at $z\sim 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\sigma$ Poisson uncertainties.101010In 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. 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 $10^{7.5}\mathrm{M}_{\odot}$ and $10^{9}\mathrm{M}_{\odot}$. 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 $10^{6}\mathrm{M}_{\odot}$ (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\sim 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 relation111111A 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 $M_{\mathrm{bh}}-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)., Magorrian et al. 1998). The 2D histogram indicates the distribution of galaxies121212We 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. in logarithm from the fiducial model using the Tiamat-125-HR halo merger trees at $z\sim 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\sim 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\sim 0.6$ compared to the local observations and we find an increasing normalization towards lower redshifts in the mass range of $10^{10}\mathrm{M}_{\odot}<M_{*}<10^{12}\mathrm{M}_{\odot}$.

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 Mutch et al. (2016), 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 $10^{8}\mathrm{M}_{\odot}<M_{*}<10^{11}\mathrm{M}_{\odot}$, 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 $M_{\mathrm{bh}}-M_{\mathrm{vir}}$ scaling relation does not get suppressed in massive haloes (at least to $M_{\mathrm{vir}}\sim 10^{14}\mathrm{M}_{\odot}$ or $M_{\mathrm{bh}}\sim 10^{9}\mathrm{M}_{\odot}$). Both black holes and stars grow from the cold gas disc.131313Radio mode accretion and merger-driven starburst are relatively less important to the growth. 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 $M_{\mathrm{bh}}>10^{8}\mathrm{M}_{\odot}$.

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., 2015; Huertas-Company et al., 2016; Stefanon et al., 2016; Davidzon et al., 2017) at redshifts 7 to ${\sim}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\sigma$ 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 $f_{\mathrm{esc,*}}$ model ($f_{\mathrm{esc,*,M16}}=\min\left[0.04\times\left(\dfrac{1+z}{6}\right)^{2.5},1.0\right]$, see Table 1) presented in Mutch et al. (2016). This model is able to reproduce the evolution of the stellar mass function at high redshift ($z>5$). The redshift dependence of $f_{\mathrm{esc,*,M16}}$ was chosen to simultaneously reproduce the normalization and flat slope of the McQuinn et al. (2011) emissivity measurement at $z\sim 5$ and the Planck 2015 optical depth measurement (Planck Collaboration et al. 2016a, hereafter Planck15). Details of the M16 model at high redshift ($z>5$) can be found in Mutch et al. (2016). In this work, we extend the model141414We note that the Tiamat halo merger trees have been improved (Poole et al. in preparation) since Mutch et al. (2016). However, the impact on galaxy formation is trivial. 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_{*}>10^{11}\mathrm{M}_{\odot}$). 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.151515Ignoring the impact on the ionizing photon escape fraction from AGN feedback.

With respect to M16, the fiducial model presented in this work employs a stronger star formation efficiency ($\alpha_{\mathrm{sf}}$) with maximized supernova feedback ($\alpha_{\mathrm{energy}}$ and $\alpha_{\mathrm{mass}}$) and more intense radio mode feedback ($k_{\mathrm{h}}$), in order to gain better agreement with the observed stellar mass function in the intermediate mass range ($10^{9}\mathrm{M}_{\odot}<M_{*}<10^{11}\mathrm{M}_{\odot}$) 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\sim 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\dfrac{\theta}{2}$, where $\theta$ represents the opening angle of AGN radiation, during the calculation of the quasar luminosity function. We note that the dependence of $\theta$ or the obscured fraction (the ratio of obscured to unobscured AGN) is very complicated. It has been suggested that $\theta$ may depend on wavelength, redshift and luminosity (Elvis et al., 1994; Hopkins et al., 2007). For simplicity, we only consider constant $\theta$. The opening angle is a free parameter in our model, chosen to reproduce the quasar luminosity function amplitude. In this work, $\theta=80$ deg is adopted. This corresponds161616An opening angle of $\theta=80$ deg corresponds to a solid angle of $\Omega=2\pi(1-\cos\dfrac{\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\dfrac{2\Omega}{4\pi}\sim 0.234$. to a fraction of visible objects,171717The 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. $f_{\mathrm{obs}}$, to be $\sim 23.4$ per cent, in agreement with Hopkins et al. (2007), who suggests a luminosity-dependent observable fraction, $f_{\mathrm{obs}}=0.26\left(L_{\mathrm{bol}}/10^{12.4}\mathrm{L}_{\odot}\right)^{0.082}$ for the B band.

We present the UV 1450 $\mathrm{\AA}$ luminosity functions of quasars at redshift 6, 5, 4 and 3, and the B-band luminosity functions at $z\sim 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 100000 bootstrap re-samples (due to the random number $t_{\mathrm{obs}}$ 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\lesssim 4$). We see that the model shows good agreement with observations across a large redshift range ($z\sim 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 the 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 $M_{1450}{\sim}{-}11$, which is much fainter than the observed faintest quasars, $M_{1450}{\sim}{-}18$ (G15). In addition, the observed population of bright quasars ($M_{1450}{<}-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.181818Also with a constant to model the sky background. 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 et al., 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 $M_{1450}\sim-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 luminosity191919Dust attenuation is not considered because these relatively faint galaxies have little dust in our models (Liu et al., 2016). (AGN+galaxy). We present the resulting luminosity function of faint objects ($M_{1450}{\sim}{-}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 100000 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 ${\sim}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{\sim}5.7$ if the high emissivity of quasars derived at $z{\sim}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 ${\sim}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 et al., 2011), we investigate the contribution of quasars to reionization, within a frame work that accounts for black hole growth and feedback on star formation.202020The 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.

4.1 Reionization model

The semi-numerical reionization code 21cmFAST (Mesinger et al., 2011) uses an excursion set formalism to identify H ii bubbles in which the cumulative number of ionizing photons is more than the number of absorbing atoms:

[TABLE]

where $N_{*}$ and $N_{\mathrm{q}}$ are the numbers of baryons in stars and quasars,212121The 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. $N_{\gamma,*}\sim 4000$ (Loeb & Barkana, 2001) and $N_{\gamma,q}$ (see the calculation in Appendix A) are the mean numbers of ionizing photons produced per baryon incorporated into the stellar or quasar components. The parameters $f_{\mathrm{esc,*}}$ and $f_{\mathrm{esc,q}}$ (see Table 1) are the escape fractions of ionizing photons produced by stars and quasars. $f_{\mathrm{obs}}{\sim}0.234$ represents the observable fraction due to obscuration222222One may also define $f_{\mathrm{esc,q}}f_{\mathrm{obs}}$ as the quasars escape fraction. (see Section 3.3). $N_{\mathrm{H_{I}}}$ is the cumulative number of atoms being ionized and $\bar{N}_{\mathrm{rec}}$ is the mean number of recombinations per baryon. Inhomogeneous recombinations are ignored, which can have a large impact (Sobacchi & Mesinger, 2014). In this work, $\bar{N}_{\mathrm{rec}}$ is set to be 0 as suggested by the high-redshift Ly$\alpha$ forest in the IGM (Bolton & Haehnelt, 2007; McQuinn et al., 2011). Expanding equation (28) gives

[TABLE]

with

[TABLE]

where $m_{*}$ is cumulative stellar mass that excludes the loss due to supernova232323The stellar mass recycled to the ISM through supernova also contributes ionizing photons in the H ii bubble. and $R$ is the radius of the H ii bubble. $\xi_{*}$ and $\xi_{\mathrm{q}}$ are the H ii ionizing efficiencies for stars and quasars, $f_{\mathrm{b}}{=}\dfrac{\Omega_{\mathrm{b}}}{\Omega_{\mathrm{m}}}$ and $Y_{\mathrm{He}}{=}0.24$ are the fraction of baryons in the Universe and the fraction of helium, and $\rho_{\mathrm{c}}\left(z\right)$ is the comoving critical density of the universe, respectively.

When the local volume around a galaxy is ionized, the UV background provides an extra heating mechanism, which modifies the baryonic fraction of the host halo (see $\chi_{\mathrm{r}}$ in equation 1). Following Sobacchi & Mesinger (2013), when the virial mass, $M_{\mathrm{vir}}$ is smaller than a filtering mass, which can be calculated through

[TABLE]

the suppression of gas becomes significant ($\chi_{\mathrm{r}}\equiv 2^{-{M_{\mathrm{filt}}}/{M_{\mathrm{vir}}}}$). Here, $z_{\mathrm{ion}}$ is the redshift when the local volume is first ionized, which is determined by the criteria given by equation (28). $J_{21}$ represents the intensity of the local UV background. This can be calculated through

[TABLE]

where $h$ is the Planck constant. $\lambda_{\mathrm{mfp}}$ is the mean-free path of ionizing photons, which is approximated by the H ii bubble radius, $R$. The parameters $\alpha_{*}=5.0$ (Loeb & Barkana, 2001) and $\alpha_{\mathrm{q}}=1.57$ (G15) are the spectral indexes for a stellar-driven and a quasar-driven spectra in the UV band. $\dot{m}_{*}$ and $\dot{m}_{\mathrm{q}}$ are the growth rates of stellar mass and black hole mass, respectively.

4.2 Reionization history

Fig. 5 shows the instantaneous ionizing emissivity averaged over the entire Tiamat simulation volume from different models as a function of redshift

[TABLE]

and

[TABLE]

where $\dot{N}_{\mathrm{ion,*}}$, $\dot{N}_{\mathrm{ion,q}}$, $\dot{m}_{\mathrm{*,tot}}$ and $\dot{m}_{\mathrm{q,tot}}$ are the ionizing emissivities for the stellar component and quasars, and the total growth rates of stellar mass excluding the loss from supernovae and black holes in the model, respectively. From the Ly$\alpha$ opacity, several measurements of the total emissivity (AGN and stars) at high redshift have been estimated (Bolton & Haehnelt, 2007; McQuinn et al., 2011; Becker & Bolton, 2013). There are relatively large uncertainties in these measurements. In this work, we compare our models using the most recent data242424The systematic error due to recombination radiation is ignored. from Becker & Bolton (2013, hereafter BB13). The top panel of Fig. 5 shows that the fiducial model agrees with the BB13 data and the bottom shows the ratio between ionizing photons from black holes and stars in the fiducial model, suggesting that during the EoR, quasars are subdominant in our model. The evolutions of the mass-weighted global neutral hydrogen fraction and the integrated Thomson scattering optical depth $\tau_{e}$ (see M16 for the calculation of $\tau_{e}$) are shown in Fig. 6. It shows that the fiducial model has a reasonable reionization history, with the mean global neutral hydrogen fraction decreasing from 90 per cent at $z\sim 10$ to 0 by $z\sim 6$ and a Thomson scattering optical depth in agreement with the latest Planck limits (Planck Collaboration et al. 2016b, hereafter Planck16).

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(\dfrac{1+z}{6}\right)^{0.5},1.0\right]$ and $f_{\mathrm{esc,q}}=0$, while quasars are the only reionization contributor in the QuasarReion model, with $f_{\mathrm{esc,q}}=1$ and $f_{\mathrm{esc,*}}=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\sim 7-5$. However, if quasars are the only reionization contributor, even with $f_{\mathrm{esc,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\sim 6$ resulting in an end at $z\sim 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 $\tau_{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\sim 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. $k_{\mathrm{c}}$, 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 $\Delta z{\lesssim}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 bias252525If it is not a selection bias, it is possible that not every galaxy hosts a central massive black hole. (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, $\sigma_{*}$, it is more difficult to resolve smaller black holes; Batcheldor et al. 2010) and showed that such a bias can lead to factors of $\geq 3$ and ${\sim}50-100$ higher normalizations of the $M_{\mathrm{bh}}-\sigma_{*}$ and $M_{\mathrm{bh}}-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. $k_{\mathrm{c}}$ in equation 19). In this case, the reconstructed black hole population becomes less massive, more efficient AGN feedback (e.g. $k_{\mathrm{h}}$ 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{\sim}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{\sim}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 $\theta=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\sim 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, $f_{\mathrm{obs}}$, 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 $f_{\mathrm{obs}}-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\sim 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\sim 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 $f_{\mathrm{esc,*}}$ and $f_{\mathrm{esc,q}}$, without any changes to the other parameters. In the top panel of Fig. 8, the left-hand panel shows the Thomson scattering optical depth, $\tau_{e}$. For comparison, shaded regions are shown corresponding to the best fit and 1$\sigma$ range of the Planck16 measurements. Based on the BB13 data at $z\sim 2-5$ and the Planck16 measurement, the top right two panels show the 68, 90 and 99 per cent confidence contours on each parameter of the best fit via the standard minimum-$\chi^{2}$ technique. The 2D histogram shows the distribution of the ratios of quasar emissivity to stellar emissivity at $z\sim 6$. We see that a lower escape fraction of ionizing photons from stars, $f_{\mathrm{esc,*}}$ 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 $f_{\mathrm{esc,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 percent ($<5$ per cent; Ciardullo et al. 2014; Matthee et al. 2016b), the model requires $f_{\mathrm{esc,q}}\sim 1.0$.

5.3 An evolving escape fraction

Although the escape fraction of ionizing photons depends on the local environment, many theoretical and observational works suggest an evolving or mass-dependent escape fraction with a decreasing average value at lower redshifts or in more massive galaxies (Kuhlen & Faucher-Giguère, 2012; Haardt & Madau, 2012; Paardekooper et al., 2013; Kimm & Cen, 2014; Wise et al., 2014; Bauer et al., 2015; Price et al., 2016). Therefore, as discussed in M16, in order to simultaneously match the normalization and flat slope of the observed ionizing emissivity at $z\lesssim 6$, and the Planck $\tau_{e}$ measurements,262626M16 constrained the model using the observed emissivities from McQuinn et al. (2011), which are lower than the BB13 data, and the Planck15 result, which has a larger optical depth than the Planck16 data. a redshift-dependent escape fraction for galaxies was proposed:

[TABLE]

We have shown two models with evolving escape fractions (fiducial and StellarReion) in Section 4. In this section, we further explore the possible evolution of the escape fraction by running the semi-analytic model with different combinations of $f_{\mathrm{esc,*}}|_{z=5}$ and $\beta$. Note that all of the ionizing photons from quasars are included in this section ($f_{\mathrm{esc,q}}=1$) in order to investigate the evolution of the stellar escape fraction with the contribution of quasars to reionization.

The bottom panels of Fig. 8 show the optical depth and the best-fitting confidence limits as functions of the normalization of $f_{\mathrm{esc,*}}(z)$, $f_{\mathrm{esc,*}}|_{z=5}$ and the scaling, $\beta$, when $f_{\mathrm{esc,q}}=1$. We see that because a larger scaling suppresses the escape fraction at lower redshifts, which results in a lower emissivity at $z\lesssim 5$, a larger normalization is required. In addition, a larger $\beta$ gives less ionizing photons at $z\sim 8$, which slows the process of reionization and consequently increases the optical depth. When $f_{\mathrm{esc,*}}|_{z=5}$ reaches 0, the model becomes quasar-dominated, returning a low $\tau_{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 $f_{\mathrm{esc,*}}|_{z=5}\sim 6$ per cent and $\beta\sim 0.5$. This corresponds to $f_{\mathrm{esc,*}}\sim 6.5$ per cent at $z=6$ with a ratio between the emissivities of quasars and stars of $\sim 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 $M_{1450}{\sim}{-}18$, while our fiducial models predict a significant population of faint quasars down to $M_{\mathrm{1450}}{\sim}{-}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 $M_{1450}$ during the Eddington state as shown in the top axis) assuming $f_{\mathrm{esc,q}}=1$ from the fiducial model using the Tiamat simulation. We see that quasars fainter than $M_{1450}=-18$ contribute approximately 80 per cent of the total emissivity at $z{\sim}7$ with a decreasing contribution towards lower redshifts (10 per cent at $z{\sim}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{\sim}7-3$ (see Fig. 3). However, at redshifts higher than $z{\sim}6$, the total emissivity from quasars is low. For example, the total emissivity at $z{\sim}7$ is five times lower than $z{\sim}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 $\epsilon=1$ or stay quiescent ($\epsilon=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 $\epsilon$ to decrease progressively when the accretion disc has been mostly consumed predicts more faint AGN with $M_{\mathrm{B}}\sim-16$ by a factor of 2 at $z\sim 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\geq 4$) in our model only reach $M_{\mathrm{1450}}{\sim}{-}23$, above which the contribution of ionizing photons is not considered. However, the G15 emissivity accounts for bright quasars up to $M_{1450}{=}{-}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<M_{1450}<-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.272727We fit the predicted quasar luminosity function at $z=2$ using a single power law at $M_{1450}>-24$ and estimate the number of missing photons with a magnitude interval of $-24<M_{1450}<-22.5$. We find the emissivity of quasars can be increased by a factor of 2 with the inclusion of the missing bright quasars. 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}\mathrm{M}_{\odot}$ lies between the light seed (${\sim}10^{2}-10^{3}\mathrm{M}_{\odot}$) from a remnant Pop III star and the massive seed (${\sim}10^{3}-10^{5}\mathrm{M}_{\odot}$) from the direct collapse of a gas cloud at early times (Greene, 2012). The massive seed mass is frequently used to initialize massive haloes ($\gtrsim 10^{6}-10^{12}\mathrm{M}_{\odot}$) in hydrodynamic simulations (Springel et al., 2005; Vogelsberger et al., 2014; Schaye et al., 2014; Feng et al., 2016) while the ${\sim}10^{3}\mathrm{M}_{\odot}$ 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_{*}>10^{9}\mathrm{M}_{\odot}$). On the other hand, the model predicts a significant number of less massive black holes with masses ${\sim}10^{5}\mathrm{M}_{\odot}$, which is more than ${\sim}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 $10^{4}\ h^{{-}1}\mathrm{M}_{\odot}$. 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\sim 7-0.6$), black hole mass function ($z\lesssim 0.5$), quasar luminosity function ($z\sim 6-0.6$), ionizing emissivity ($z\sim 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 ${\sim}23.4$ per cent, allows the model to reproduce the observed quasar luminosity function across a large redshift range ($z\sim 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\sim 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{\sim}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 ($f_{\mathrm{esc,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(\dfrac{1+z}{6}\right)^{0.5},1\right]$. This corresponds to quasars contributing 10 per cent of the total ionizing photons at $z\sim 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).

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 $M_{\mathrm{BH}}$, 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, $M_{1450}$ 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 ($\Delta{M}_{\mathrm{BH}}<M_{\mathrm{BH}}$, see Fig. 11), for the sake of calculation speed, we estimate the mean number of ionizing photons produced per black hole, $N_{\gamma,q}$ with the bolometric luminosity at the beginning of accretion. We calculate $N_{\gamma,q}$ as follows:

We calculate the UV magnitude, $M_{1450}$ using equations (22)-(27). 2. 2.

We calculate the UV flux with $M_{1450}$ in units of $\mathrm{erg\ s^{-1}Hz^{-1}}$ through

[TABLE] 3. 3.

We calculate the flux at Lyman limit following G15

[TABLE]

where $\alpha_{q,\mathrm{optical}}=0.44$ and $\alpha_{\mathrm{q}}=1.57$ correspond to a double power-law AGN SED. 4. 4.

We calculate the instantaneous emissivity by

[TABLE] 5. 5.

The duration of accreting mass, $\Delta M_{\mathrm{BH}}$ can be calculated through

[TABLE]

Therefore, 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

[TABLE]

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, $L_{\mathrm{bol}}$ increases exponentially (equation 22). Since the B-band bolometric correction, $k_{\mathrm{B}}$, 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 (40) 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 $k_{\mathrm{B}}$, $N_{\gamma,q}$ can be analytically calculated by integrating the AGN light curve. However, we note that since $k_{\mathrm{B}}\left(t\right)\leq k_{\mathrm{B}}|_{t=t_{\mathrm{acc}}}$, this approximation overestimates $N_{\gamma,q}$.

Fig. 12 presents the evolution of emissivity, neutral hydrogen fraction and optical depth for different QuasarReion models assuming constant $\dot{N}_{\mathrm{ion}}$ (QuasarReion) and $k_{\mathrm{B}}$ (QuasarReion_kB), respectively. Since the time interval between two snapshots is much smaller than the Eddington accretion time-scale ($t_{\mathrm{Edd}}\sim 450\mathrm{Myr}$), the black hole mass increment is still within the linear regime and therefore we see that the impact from the calculation of $N_{\gamma,q}$ is not significant.