# Constraining the Contribution of Galaxies and Active Galactic Nuclei to Cosmic Reionization

Shintaro Yoshiura, Kenji Hasegawa, Kiyotomo Ichiki, Hiroyuki Tashiro, Hayato Shimabukuro, Keitaro Takahashi

► To cite this version:

Shintaro Yoshiura, Kenji Hasegawa, Kiyotomo Ichiki, Hiroyuki Tashiro, Hayato Shimabukuro, et al.. Constraining the Contribution of Galaxies and Active Galactic Nuclei to Cosmic Reionization. Monthly Notices of the Royal Astronomical Society, 2017, 471 (3), pp.3713-3726. 10.1093/mn-ras/stx1754 . hal-01645173

HAL Id: hal-01645173

https://hal.science/hal-01645173v1

Submitted on 24 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Constraining the contribution of galaxies and active galactic nuclei to cosmic reionization

Shintaro Yoshiura,<fnref n="1" />* Kenji Hasegawa,<fnref n="2" /> Kiyotomo Ichiki,<fnref n="2" />,<fnref n="3" /> Hiroyuki Tashiro,<fnref n="2" /> Hayato Shimabukuro<fnref n="4" /> and Keitaro Takahashi<fnref n="1" />

Accepted 2017 July 11. Received 2017 July 11; in original form 2016 October 15

## ABSTRACT

Understanding the detailed process of cosmic reionization is one of the remaining problems in astrophysics and cosmology. Here we construct a model of cosmic reionization that includes contributions from high-z galaxies and active galactic nuclei (AGNs) and calculate reionization and thermal histories with the model. To keep the model general and realistic, we vary the escape fraction of ionizing photons, $f_{esc}$, and the faint-end slope of the AGN luminosity function at high redshifts, $\alpha_{hz}$, within constraints from the observed cosmic star formation history and observed bright-end UV luminosity functions at $z \le 6$. Additionally, we model the spectral energy distribution (SED) of AGNs, which depends on the Eddington ratio and the black hole mass. By comparing the computed reionization histories with the observed H I fractions and the optical depth for Thomson scattering from Planck, we find that $\alpha_{hz} > -1.5$ and $f_{esc} < 0.15$ are favoured when we employ the bright-end luminosity function obtained by Giallongo et al. Our result suggests that an AGN-dominated model with an abundance of faint AGNs as large as the estimate by Giallongo et al. is allowed only if the contribution from high-z galaxies is almost negligible, while a galaxy-dominant model is also allowed. We also find that the shape of the SED has a significant impact on the thermal history. Therefore it is expected that measurements of the thermal state of the intergalactic medium (IGM) will provide useful information on the properties of ionizing sources.

**Key words:** dark ages, reionization, first stars.

## 1 INTRODUCTION

It is generally believed that the Universe became neutral at $z \approx 1100$ (cosmic recombination) and was gradually ionized again by ionizing sources formed following cosmic recombination. The time evolution of this global transition, called cosmic reionization, has been widely explored. One observational approach to estimating neutral fractions at high redshifts is to measure absorption features imprinted on high-z quasar (QSO) spectra. The H I Lyman $\alpha$ (Ly$\alpha$) absorption lines imprinted on QSO spectra have indicated that hydrogen was almost reionized by $z \sim 6$ (Fan et al. 2006). Similar observations measuring He II absorption lines have also shown that helium reionization was mostly completed by $z \approx 2.7$ (Becker et al. 2011; Syphers & Shull 2014; Worseck et al. 2014a). Observing the redshift evolution of the Ly$\alpha$ luminosity function (Ly$\alpha$ LF) can also provide valuable information on neutral hydrogen fractions

during the epoch of reionization (EoR), because the detectability of Ly$\alpha$ radiation is very sensitive to the amount of neutral hydrogen. Recent observations for $z \ge 5.7$ Ly$\alpha$ emitters (LAEs) have shown that the number density of LAEs decreases with redshift (Hu et al. 2010; Ota et al. 2010; Ouchi et al. 2010; Kashikawa et al. 2011; Konno et al. 2014). This fact implies that the neutral hydrogen fraction increases from $z \sim 6$ to $z \sim 7$ (Choudhury et al. 2015). Furthermore, measurements of the cosmic microwave background (CMB) optical depth to Thomson scattering have imposed an additional constraint on the reionization history. The most recent observation by Planck has shown that the optical depth $\tau_e$ is $0.058 \pm 0.012$, which corresponds to an instantaneous reionization redshift of $z_t = 8.5 \pm 0.9$ (Planck Collaboration XLVII 2016).

Despite advancing understanding of the history of cosmic reionization, there is no crucial observation of the type of ionizing source responsible for cosmic reionization. Revealing the type of ionizing source is also important for understanding the process of reionization, since the mean free path of ionizing photons is essential for the reionization evolution.

<footnote n="1">¹Department of Physics, Kumamoto University, Kumamoto, Japan</footnote>

<footnote n="2">²Department of Physics, Nagoya University, Aichi, Japan</footnote>

<footnote n="3">³Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya 464-8602, Japan</footnote>

<footnote n="4">⁴Sorbonne Universites, UPMC, LERMA, Observatoire de Paris, PSL research university, CNRS, F-75014 Paris, France</footnote>

<footnote>* E-mail: 161d9002@st.kumamoto-u.ac.jp</footnote>
High-z galaxies are thought to be plausible candidates for ionizing sources, since a number of galaxies have already been discovered at $z > 6$ (Ouchi et al. 2009; Oesch et al. 2014; Bouwens et al. 2015). In addition, some recent observations have indicated that faint-end slopes of ultraviolet luminosity functions (UVLFs) during the EoR are steep and extended down to an absolute UV magnitude of $M_{UV} \sim -16$ (Atek et al. 2015; Ishigaki et al. 2015). These facts strengthen the possibility that high-z galaxies were the main ionizing sources of reionization. Robertson et al. (2015) have indeed shown that the Thomson scattering optical depth and H I fractions inferred from QSO spectra and LAE observations can simultaneously be achieved solely by high-z galaxies if the escape fraction of ionizing photons from the galaxies ($f_{esc}$) is as high as $\sim 0.2$. However, the typical escape fraction during the EoR is highly uncertain.

At intermediate redshift $z \sim 3$, Lyman continuum (LyC) photons from galaxies have been detected directly and the escape fractions have been evaluated as $f_{esc} \sim 0.01-0.1$ (Steidel, Pettini & Adelberger 2001; Iwata et al. 2009; Mostardi et al. 2015; Matthee et al. 2017; Naidu et al. 2016), although some of the LyC emitters may be low-z contaminants (Siana et al. 2015). In contrast to galaxies at intermediate redshifts, no ionizing photons from galaxies at the EoR have been detected, due to strong attenuation by the IGM. However, Dunlop et al. (2013) have claimed that spectral energy distributions (SEDs) of high-z galaxies imply an escape fraction $f_{esc} \approx 0.1-0.2$, with some assumptions. The escape fraction has also been assessed by numerical simulations (e.g. Gnedin, Kravtsov & Chen 2008; Razoumov & Sommer-Larsen 2010; Yajima, Choi & Nagamine 2011; Kimm & Cen 2014; Paadekooper et al. 2014; Wise et al. 2014). Although most of the simulations show a similar trend, i.e. the escape fraction decreases with the mass of the galaxy, there is a variance of $\sim 1-2$ orders of magnitude among the simulation results. Thus, as yet there is no consensus on the value of the escape fraction and this fact leads to significant uncertainty regarding the contribution of high-z galaxies to cosmic reionization.

Active galactic nuclei (AGNs) are also expected to have been another type of ionizing source. Since photons with energies above 54.4 eV are required for doubly ionizing helium, AGNs are the most promising ionizing sources at the epoch when helium in the IGM is fully ionized, i.e. $z \sim 2.7$ (Becker et al. 2011; Worseck et al. 2014a). However, it has traditionally been believed that the contribution of AGNs to reionization is less important, because of a rapid decrease in their abundance at $z > 3$ (Masters et al. 2012; Ueda et al. 2014). Very recently, Giallongo et al. (2015) have found faint AGNs at $z \sim 4-6$. This suggests that the contribution of AGNs to reionization could increase by orders of magnitude. Although it is still under debate, the contribution of AGNs might overcome that of galaxies if this claim were valid beyond $z \approx 6$. Unfortunately, while there is currently no direct observation for faint AGNs at $z > 6$, the contribution of AGNs could be imprinted on the H I 21-cm line signal (Baek et al. 2010). Hence the Square Kilometre Array (SKA) and its pathfinders will yield information on the abundance of AGNs at $z > 6$.

Madau & Haardt (2015, hereafter MH15) calculated the reionization history, extrapolating the evolution of AGN abundance reported by Giallongo et al. (2015) up to $z = 12$ and neglecting the contribution of galaxies. As a result, they have found that the ‘AGN-dominant scenario’ can satisfy the observed H I/He II fractions and Thomson scattering optical depth simultaneously. However, as described above, if $f_{esc}$ is high then the additional contribution of galaxies would result in the overproduction of ionized hydrogen and electrons at $z > 6$. Hao, Yuan & Wang (2015) also performed calculations including both galaxies and (faint) AGNs and found

that their model is consistent with the EoR measurements, Thomson scattering optical depth and helium reionization. However, they neglected redshifted high-energy photons, which ionize the IGM effectively at lower redshifts, in particular helium reionization.

Mitra, Choudhury & Ferrara (2016) investigated the contribution of AGNs to reionization history by using their semi-analytic reionization model. They have pointed out that their high AGN emissivity model can reionize H I. Their model also leads to more rapid evolution of Lyman-limit systems than in observations, while reionization evolution through the combination of lower AGN emissivity and galaxies with an escape fraction $\sim 10$ per cent shows good agreement with the observation results. Khaire et al. (2016) have shown numerically that abundant faint AGNs or $f_{esc}$ rising with redshift are required to explain the constant photoionization rate at $3 < z < 5$ reported by observations.

In those numerical studies, SEDs of AGNs have been assumed to be simple power-law spectra; however, AGNs at $z \sim 1$ show SEDs with double peaks at energy ranges in the UV ($\sim 10$ eV) and X-ray ($\sim 10$ keV) (e.g. Zheng et al. 1997; Laor et al. 1997), which can be explained theoretically as the radiation from a system composed of the accretion disc and coronal gas (e.g. Kawaguchi, Shimura & Mineshige 2001). It is well known that the shape of the SED controls the efficiency of the ionization and heating processes of the IGM. Therefore, the SED shape is an important factor in calculation of the EoR. In spite of its importance, the SED shape of AGNs in the reionization era ($z > 6$) is highly uncertain, especially in the X-ray regime, where the heating process effectively contributes to the thermal history. Therefore it is worth elucidating how much the SED shape impacts both the reionization and thermal histories.

In this theoretical work, we aim to constrain the contribution of galaxies and AGNs to cosmic reionization, as well as investigating the effect of the SED shape of AGNs on reionization and thermal histories. For this purpose, we first construct a model of the redshift evolution of AGN emissivity, in which we consider the redshift evolution of the AGN LF and SED depending on the central black hole mass $M_{BH}$ and Eddington ratio $\lambda_{Edd} \equiv L/L_{Edd}$, where $L_{Edd}$ is the Eddington luminosity. We then compute the H/He reionization and thermal histories with different parameter sets of $f_{esc}$ and AGN abundance. To clarify the impact of the SED shape on the reionization and thermal histories, we employ not only the newly developed SED model but also a simple power-law SED model. Finally, we search for the possible contribution of galaxies and AGNs to reionization through a comparison between the computed reionization histories and observations.

The structure of this article is as follows. In Section 2, we describe the method for computing the hydrogen and helium reionization histories. In Section 3, we show our numerical results and compare them with observed H I/He II fractions. We discuss the implication of our results and comment on the model dependence in Section 4. Finally, a summary is given in Section 5. Throughout this article, we assume $\Lambda$CDM cosmology with ($\Omega_m, \Omega_\Lambda, \Omega_b, H_0) = (0.308, 0.692, 0.02237, 67.80 \text{ km s}^{-1} \text{ Mpc}^{-1}$) (Planck Collaboration XXIII 2015).

## 2 METHODOLOGY

The evolution of the ionization fraction $f_i$ (for species $i = \text{H I, He II, He III}$) can be written as (e.g. Madau, Haardt & Rees 1999)

$$ \frac{d f_{\text{H II}}}{dt} = \frac{1}{n_{\text{H}}(1+z)^3} \frac{d N_{\gamma, \text{H I}}}{dt} - \alpha_{\text{B, H II}} C n_{\text{e}} (1+z)^3 f_{\text{H II}}, \quad (1) $$
$$ \frac{\mathrm{d}f_{\mathrm{He II}}}{\mathrm{d}t} = \frac{1}{n_{\mathrm{He}}(1+z)^3} \left( \frac{\mathrm{d}N_{\gamma, \mathrm{He I}}}{\mathrm{d}t} - \frac{\mathrm{d}N_{\gamma, \mathrm{He II}}}{\mathrm{d}t} \right) - C n_e (1+z)^3 (\alpha_{\mathrm{B, He II}} f_{\mathrm{He II}} - \alpha_{\mathrm{B, He III}} f_{\mathrm{He III}}), \quad (2) $$

$$ \frac{\mathrm{d}f_{\mathrm{He III}}}{\mathrm{d}t} = \frac{1}{n_{\mathrm{He}}(1+z)^3} \frac{\mathrm{d}N_{\gamma, \mathrm{He III}}}{\mathrm{d}t} - \alpha_{\mathrm{B, He III}} C n_e (1+z)^3 f_{\mathrm{He III}}, \quad (3) $$

where $\alpha_{\mathrm{B}, i}$ is the case B recombination coefficient for each species and $n_{\mathrm{H}}$, $n_{\mathrm{He}}$ and $n_e$ are the comoving number densities for hydrogen, helium and electrons, respectively. Here C is the clumping factor of the IGM, for which we use $C=3$, as shown by hydrodynamic simulations taking into account photo-heating of the IGM (Pawlik, Schaye & van Scherpenzeel 2009; Raicević & Theuns 2011). In equations (1)–(3), $dN_{\gamma, i}/dt$ is the photoionization rate per unit volume for the $i$th species. In this article, we consider two types of ionizing source: using the photoionization rate per unit volume of star-forming galaxies, $dN_{\mathrm{AGN}, i}/dt$, and that of AGNs, $dN_{\mathrm{AGN}, i}/dt$, the total photoionization rate per unit volume is expressed as $dN_{\gamma, i}/dt = dN_{\mathrm{AGN}, i}/dt + dN_{\mathrm{AGN}, i}/dt$. We will describe how we consider each contribution in the following subsections.

Ionizing photons emitted by stars and AGNs also heat up the IGM. We should note that previous studies (e.g. MH15) omitted thermal evolution and assumed a fixed temperature. However, it is expected that the thermal history is sensitive to the shape of the SED and AGN abundance. Therefore we calculate the thermal and ionization evolution of the IGM consistently. The evolution of the IGM temperature $T_e$ follows

$$ \frac{\mathrm{d}T_e}{\mathrm{d}t} = \frac{2}{3} \frac{\mu}{k_B n_B} \left( \sum_i \epsilon_{\mathrm{heat}, i} - \epsilon_{\mathrm{cool}} \right) - 2H T_e - \mu T_e \frac{d\mu^{-1}}{dt}, \quad (4) $$

where H is the Hubble parameter, $k_B$ is the Boltzmann constant, $\mu$ is the mean molecular weight, $n_B$ is the proper baryon number density and $\epsilon_{\mathrm{cool}}$ is the IGM cooling rate. We take into account the IGM cooling processes summarized in the appendix of Fukugita & Kawasaki (1994). In equation (4), $\epsilon_{\mathrm{heat}, i}$ describes the heating rate for each species $i$. We also decompose the heating rate into two components: from star-forming galaxies, $\epsilon_{\mathrm{*,i}}$, and from AGNs, $\epsilon_{\mathrm{AGN}, i}$, as $N_{\gamma}$ in the case of calculations of the ionization fraction.

When we solve equations (1)–(4), we employ some artificial treatments to bypass some numerical inconveniences. As described later, we assume that all ionizing photons from star-forming galaxies are consumed to ionize the local IGM. Due to this treatment, $f_{\mathrm{H I}}$ can exceed 1.0 and $f_{\mathrm{H II}}$ essentially becomes zero. Unfortunately, this unrealistic behaviour results in numerical oscillations of $f_{\mathrm{H I}}$ and $f_{\mathrm{H II}}$. In order to avoid such unrealistic behaviour, we set the lower limit of the neutral hydrogen fraction as $f_{\mathrm{H I, I}} = 10^{-4}$. However, this artificial limit leads to the continual heating of the IGM. Thus we also set the upper limit of the temperature as $T_{e, u} = 20\,000$ K.

## 2.1 Ionizing photons from stars

In most previous studies on reionization, star-forming galaxies have been considered as the main sources of ionizing photons. In a star-forming galaxy, a huge amount of photons is produced from stars and some fraction of these photons can escape from the galaxy and ionize the surrounding IGM. Thus, to evaluate the photoionization rate per unit volume by star-forming galaxies ($dN_{\mathrm{AGN}, i}/dt$), we need

to know the value of the escape fraction $f_{\mathrm{esc}}$. However, $f_{\mathrm{esc}}$ is highly uncertain, as described in Section 1. Therefore we take $f_{\mathrm{esc}}$ as a free parameter in our model. Since $f_{\mathrm{esc}}$ controls the contribution of high-z galaxies to the cosmic reionization process, we obtain the constraint on $f_{\mathrm{esc}}$ from a comparison with observational data of the EoR. For simplicity, we assume that $f_{\mathrm{esc}}$ does not depend on the halo mass and redshift. Using $f_{\mathrm{esc}}$, we assume that the proper number density of photons emitted per unit time from star-forming galaxies at redshift $z$ is proportional to the star formation rate density (SFRD),

$$ \dot{n}_{*\nu} = (1+z)^3 f_{\mathrm{esc}} \dot{\rho}_{*} \gamma_{\nu}, \quad (5) $$

where $\dot{\rho}_{*}$ is the comoving SFRD, for which we adopt the fitting formula by Madau & Dickinson (2014):

$$ \dot{\rho}_{*}(z) = 0.015 \frac{(1+z)^{2.7}}{1 + [(1+z)/2.9]^{3.6}} [\mathrm{M}_{\odot} \mathrm{yr}^{-1} \mathrm{Mpc}^{-3}]. \quad (6) $$

In equation (5), $\gamma_{\nu}$ is the intrinsic number of photons produced by stars per solar mass at a frequency $\nu$. For $\gamma_{\nu}$, we use the corresponding formula in Choudhury & Ferrara (2005) derived by a population synthesis calculation,

$$ \int_{v_1}^{v_2} \mathrm{d}\nu \, \gamma_{\nu} = \begin{cases} 5.43 \times 10^{60} / \mathrm{M}_{\odot}, & (v_1, v_2) = (v_{\mathrm{H I}}, v_{\mathrm{He I}}), \\ 2.61 \times 10^{60} / \mathrm{M}_{\odot}, & (v_1, v_2) = (v_{\mathrm{He I}}, v_{\mathrm{He II}}), \\ 0.01 \times 10^{60} / \mathrm{M}_{\odot}, & (v_1, v_2) = (v_{\mathrm{He II}}, v_{\mathrm{max,*}}), \end{cases} \quad (7) $$

where $v_{\mathrm{max,*}}$ is the upper limiting frequency of photons from galaxies in our calculation and is set to $h_p v_{\mathrm{max,*}} = 100$ eV. Here, $h_p$ is the Planck constant.

As mentioned above, the ionizing photons from star-forming galaxies are in the range 13.6 eV–100 eV. Since these photons have short mean free paths, they are consumed to ionize a local region of the IGM. This means that all photons from star-forming galaxies can be assumed to contribute to ionizing the local IGM. We note that this assumption is not reasonable for H I ionizing photons at lower redshift $z < 5$, because the mean free path is of the order of 100 proper Mpc at $z=3$ (shown in fig. 10 of Worseck et al. 2014b). Following previous studies (e.g. Hao et al. 2015), we work with this assumption but also employ an additional modification that ensures the conservation of the number of ionizing photons. With this assumption, a photon is inevitably absorbed by either an H I atom, He I atom or He II atom. Therefore, to calculate photoionization rates, it is necessary to evaluate how great an extent of photons from star-forming galaxies is consumed for the ionization of the $i$th species, otherwise double-counting of ionizing photons occurs. The relative number of photons consumed for the ionization of the $i$th species should be proportional to the relative opacity of the corresponding species (e.g. Susa 2006; Umemura et al. 2012). Taking into account the frequency dependence of the ionization cross-sections for hydrogen and helium, we can write the photoionization rate per unit volume for the $i$th species as

$$ \frac{\mathrm{d}N_{\mathrm{AGN}, i}}{\mathrm{d}t} = \int_{v_i}^{\infty} \mathrm{d}\nu \frac{n_i \sigma_i \dot{n}_{*\nu}}{n_{\mathrm{H I}} \sigma_{\mathrm{H I}} + n_{\mathrm{He I}} \sigma_{\mathrm{He I}} + n_{\mathrm{He II}} \sigma_{\mathrm{He II}}}, \quad (8) $$

where $n_i$ and $\sigma_i$ are the number density and cross-section for the $i$th species, respectively. The Lyman limit frequency of the $i$th species is expressed as $v_i$, i.e. $h_p v_{\mathrm{H I}} = 13.6$ eV, $h_p v_{\mathrm{He I}} = 24.5$ eV and $h_p v_{\mathrm{He II}} = 54.4$ eV.

Next, we consider the heating rate due to radiation from star-forming galaxies $\epsilon_{\mathrm{*,i}}$. When a photon with frequency $\nu$ ionizes a particle of the $i$th species, the energy $h_p(\nu - v_i)$ contributes to heat the IGM gas. Thus, the photo-heating rate per unit volume $\epsilon_{\mathrm{*,i}}$ for
the ith species is given by

$$ \epsilon_{*,i} = \int_{v_i}^{\infty} dv \frac{n_i \sigma_i h_p (v - v_i) \dot{n}_{*v}}{n_{H_i} \sigma_{H_i} + n_{He_i} \sigma_{He_i} + n_{He_{ii}} \sigma_{He_{ii}}} \quad (9) $$

As shown by equation (7), we only know the integrated $\gamma_v$. Therefore, we perform the integration in equations (8) and (9) by using the cross-sections averaged over each frequency range.

## 2.2 Ionizing photons from AGNs

Photons from AGNs are also important in the reionization process of the IGM in the early Universe. X-ray photons have much larger mean free paths than UV photons and can ionize and heat up the IGM at cosmological distances and hence the assumption of ‘local ionization’ adopted for star-forming galaxies is invalid. Therefore, to evaluate the ionization and heating rates at a redshift $z$, we need to take into account X-ray photons coming from AGNs at redshifts higher than $z$. Then, the photoionization and photo-heating rates by AGNs, respectively, are given by (e.g. Madau et al. 1999)

$$ \frac{dN_{\text{AGN},i}}{dt}(z) = \int_{v_i}^{v_{\text{max,AGN}}} \frac{dv}{h_p v} \sigma_i n_i (1+z)^3 F(z, v), \quad (10) $$

and

$$ \epsilon_{\text{AGN},i}(z) = \int_{v_i}^{v_{\text{max,AGN}}} \frac{dv}{h_p v} \sigma_i(v) n_i h_p (v - v_i) (1+z)^3 F(z, v), \quad (11) $$

where $F(z, v)$ represents the specific energy flux from AGNs at redshifts higher than $z$ at a frequency $v$. The upper bound of the integration, $v_{\text{max,AGN}}$, is set to be $2.4 \times 10^{20}$ Hz, which corresponds to an energy of $10^3$ keV. The specific energy flux at redshift $z$ is calculated as

$$ F(z, v) = \int_z d z' \epsilon_c(z', v') \frac{c(1+z)^2}{H} \exp(-\tau_v(z, z', v')), \quad (12) $$

where $c$ is the speed of light and $v' = (v(1+z'))/(1+z)$. Here, $\epsilon_c(z, v)$ is the comoving emissivity, i.e. the specific radiation energy density per unit time coming from AGNs at redshift $z$ and frequency $v$. As mentioned in Section 1, previous studies focusing on the contribution of AGNs to reionization (e.g. Madau & Haardt 2015; Hao et al. 2015; Khaire et al. 2016; Mitra et al. 2016) have assumed power-law spectra up to the hard X-ray regime ($\gtrsim 10$ keV) and the evolution of $\epsilon_c$ itself for simplicity. However, the shape of the SED and the comoving emissivity are strongly expected to depend on the evolution of AGNs. Therefore, we take into account the evolution of the AGN LF and SED depending on $M_{BH}$ and $\lambda_{Edd}$ to elucidate the AGN contribution at the EoR in more detail.

Using the SED $L_v$ and UVLF $\Phi(z, M_{UV})$ of AGNs, where $M_{UV}$ is the absolute magnitude at a wavelength of 1450 Å, we can write $\epsilon_c(z, v)$ as

$$ \epsilon_c(z, v) = \int_{M_{\text{UV}}^{\text{min}}}^{M_{\text{UV}}^{\text{max}}} L_v(M_{\text{UV}}) \Phi(z, M_{\text{UV}}) \, dM_{\text{UV}}. \quad (13) $$

The upper bound of the integration is fixed to be $M_{\text{UV}}^{\text{max}} = -34$, while the lower bound is set to $M_{\text{UV}}^{\text{min}} = -18$ as a fiducial value. The SED $L_v$ is a function of $M_{UV}$. We will describe how $M_{BH}$ and $\lambda_{Edd}$ are related to $M_{UV}$ for performing integration with respect to $M_{UV}$ (Section 2.2.2) and how the evolution of the UVLF $\Phi$ [Mpc$^{-3}$ mag$^{-1}$] is modelled (Section 2.2.1). Hereafter, ‘UV’ denotes a wavelength of 1450 Å, unless otherwise noted.

<table><caption>**Table 1.** Parameters of our AGN UVLF model at $z = 2.0, 4.25, 4.75, 5.75, 8.00$ and $10.0$. The faint-end slope $\alpha_{hz}$ at $z > 4.25$ is a parameter. The log $\Phi_*$ values are for $z_{AGN} = 6.0$.</caption><thead><tr><th>z</th><th>α</th><th>β</th><th>M<sub>*</sub></th><th>log Φ<sub>*</sub></th></tr></thead><tbody><tr><td>2.00</td><td>-1.5</td><td>-3.76</td><td>-25.4</td><td>-5.86</td></tr><tr><td>4.25</td><td>α<sub>hz</sub></td><td>-3.13</td><td>-23.2</td><td>-4.89</td></tr><tr><td>4.75</td><td>α<sub>hz</sub></td><td>-3.14</td><td>-23.6</td><td>-5.36</td></tr><tr><td>5.75</td><td>α<sub>hz</sub></td><td>-3.35</td><td>-23.4</td><td>-5.38</td></tr><tr><td>8.00</td><td>α<sub>hz</sub></td><td>-3.73</td><td>-23.0</td><td>-5.43</td></tr><tr><td>10.0</td><td>α<sub>hz</sub></td><td>-3.99</td><td>-22.8</td><td>-5.47</td></tr></tbody></table>

The optical depth for a photon with an observed frequency $v$ travelling from $z'$ to $z$, $\tau_v(z, z', v)$, is given by

$$ \tau_v(z, z', v) = \int_z^{z'} dz \frac{c(1+z)^2}{H} \left( f_{H_i} n_{H_i} \sigma_{H_i}(\hat{v}) + f_{He_i} n_{He_i} \sigma_{He_i}(\hat{v}) + f_{He_{ii}} n_{He_{ii}} \sigma_{He_{ii}}(\hat{v}) \right), \quad (14) $$

where $\hat{v} = v(1 + \hat{z})/(1 + z)$. We assume that the $H_i$ distribution in the Universe is uniform to obtain equation (14), since cosmological simulations usually show that the IGM tends to be uniform with increasing redshift (e.g. Iliev et al. 2006). However, we should note that the inhomogeneity of the IGM grows during the EoR and might affect the process of the EoR (Khaire et al. 2016). Unfortunately, the distribution of Lyman limit systems at the EoR, which contributes to the inhomogeneity of the IGM, is still too unknown to consider appropriately in this study (but see also Inoue et al. 2014). We will discuss the impact of IGM inhomogeneity on our results in Section 4.

### 2.2.1 AGN luminosity function

In this article, we construct an AGN UVLF model consistent with recent observations. According to Croom et al. (2009), we express the UVLF in the shape of a double power law as

$$ \Phi(z, M_{\text{UV}}) = \frac{\Phi_* \exp(-z/z_{\text{AGN}})}{10^{0.4(\alpha+1)(M_{\text{UV}}-M_*)} + 10^{0.4(\beta+1)(M_{\text{UV}}-M_*)}} \quad (15) $$

Here, the parameters $\Phi_*$, $M_*$, $\alpha$ and $\beta$, respectively, control the amplitude, characteristic $M_{UV}$ where the slope switches, faint-end slope and bright-end slope of the UVLF. These parameters are set to reproduce observed UVLFs at $z < 4$ (Croom et al. 2009; Giallongo et al. 2015). We use $\alpha = -1.5$ at $z = 2$ (Croom et al. 2009). In equation (15), the decrease in AGN number density at high $z$ is controlled by $z_{AGN}$, for which we assume $z_{AGN} = 6$.

The details of UVLFs at higher redshifts are still under debate, especially in their faint-end parts. Our aim in this article is to evaluate these contributions to the EoR process. Therefore, we introduce a parameter $\alpha_{hz}$ as the faint-end slope parameter at $z > 4.25$. We then interpolate $\alpha$ between $z = 2$ and $z = 4.25$. For other UVLF parameters, i.e. $\Phi_*, M_*$, and $\beta$ at high redshifts $z > 4.25$, we take $\log(1+z)$ interpolation using the values at redshifts $z = 4.25, 4.75$ and $5.75$ represented in Giallongo et al. (2015). We extrapolate these parameters further for $z > 5.75$. The UVLF parameters used in this article are listed in Table 1. We show our AGN UVLF models for different redshifts in Fig. 1. Our UVLF models are consistent with observations at each redshift. As shown by the right panel of Fig. 1, the contribution of faint AGNs to reionization increases with decreasing $\alpha_{hz}$ and increasing $M_{UV}^{\text{min}}$.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** AGN UVLFs at $z < 2$ (left), 4.25 (centre) and 5.75 (right). Observational data are taken from Jiang et al. (2008), Croom et al. (2009), Glikman et al. (2011) and Giallongo et al. (2015). In the left panel, the red solid line and the blue dashed line, respectively, represent models of LFs at $z = 0.5$ and $2.0$ with $\alpha_{hz} = -1.5$. In the central and right panels, modelled LFs with $\alpha_{hz} = 0.0, -1.5$ and $-2.2$, respectively, are shown by the red thick solid, blue thick dashed and green thin dashed lines.</figcaption></figure>

### 2.2.2 AGN spectral energy distribution

In this section, we derive an SED model depending on $M_{BH}$ and $\lambda_{Edd}$ based on results by Kawaguchi et al. (2001), who have computed the SED of a system composed of a standard accretion disc and coronal gas numerically. Hereafter we call this model the AC model (accretion disc + corona). We first introduce a shape function $\phi(E)$:

$$vL_{\nu} = L_{UV}\phi(E), \qquad (16)$$

$$\phi(E) = \left\{ \left( \frac{E}{E_1} \right)^{\gamma_a} + \left( \frac{E}{E_1} \right)^{\gamma_b} \right\}^{-1} + A_c \left\{ \left( \frac{E}{E_2} \right)^{\gamma_c} + \left( \frac{E}{E_2} \right)^{\gamma_d} \right\}^{-1}, \qquad (17)$$

where $L_{UV}, A_c, \gamma_a, \gamma_b, \gamma_c, \gamma_d, E_1$ and $E_2$ are the model parameters selected to reproduce the SED with $M_{BH} = 3 \times 10^9 M_{\odot}$ and $\lambda_{Edd} = 0.05$ obtained in Kawaguchi et al. (2001).<fnref n="1" /> Carrying out $\chi^2$ fitting, the parameters are found to be $A_c = 0.397, \gamma_a = -0.0191, \gamma_b = 0.685, \gamma_c = -0.255, \gamma_d = 0.211, E_1 = 9.44 \times 10^3$ keV, $E_2 = 1.22 \times 10^{-3}$ keV and $L_{UV} = 3.4 \times 10^{45}$ erg s$^{-1}$.

We employ the theoretically motivated scaling relations shown in Kawaguchi et al. (2001) to take the dependence of the SED on $M_{BH}$ and $\lambda_{Edd}$ into account. According to the standard accretion disc model, the peak frequency of the big blue bump (the bump at $E \sim 10$ eV in Fig. 2) follows the scaling relation

$$v_p \propto M_{BH}^{-1/4} \lambda_{Edd}^{1/4}. \qquad (18)$$

In addition, the scaling of the total luminosity is simply given by

$$L = \lambda_{Edd} L_{Edd} \propto \lambda_{Edd} M_{BH}. \qquad (19)$$

Kawaguchi et al. (2001) have indeed shown that these simple relations explain approximately the behaviour of their computed SEDs.

Owing to the scaling relations, we can uniquely convert $L_{\nu}(M_{BH}, \lambda_{Edd})$ into $L_{\nu}(M_{UV})$ appearing in equation (13). For the redshift evolution of $\lambda_{Edd}$, we adopt $\lambda_{Edd}(z) = \min[0.0091 \times (1+z)^{1.9}, 1.0]$, based on observations by Shen & Kelly (2012).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: weighted average SED for AGNs at $z = 1$. The red thick solid and black thin curves represent the AC and PL models, respectively. The SED with $M_{BH} = 3 \times 10^9 M_{\odot}$ and $\lambda_{Edd} = 0.05$ computed by Kawaguchi et al. (2001, green dotted curve) and the observed composite SED (Laor et al. 1997; Zheng et al. 1997, blue dashed curve) are also plotted for reference. Bottom: weighted average SED for AGNs at $z = 6$ with $\alpha_{hz} = -1.5$. The AC and PL models are shown by red thick and black thin solid curves, respectively. The red solid line is obtained by performing the integrations in equation (21) over $-32 < M_{UV} < -15$. The blue dashed, green dot-dashed and orange dotted curves represent the contributions from AGNs with $M_{UV} < -23, -23 \le M_{UV} < -18$ and $-18 \le M_{UV} < -15$, respectively. The sum of these three curves corresponds to the red solid curve. We note that the contribution from faint AGNs is controlled by $\alpha_{hz}$.</figcaption></figure>

<footnote n="1">¹ To be exact, there is a degeneracy between $\lambda_{Edd}$ and the radiation efficiency $\eta$. We implicitly assume $\eta = 0.1$ for the fiducial model.</footnote>
In order to discuss the impact of the shape of the SED on the reionization and thermal histories, we also perform calculation with a simple power-law SED (hereafter the PL model) proposed by Lusso et al. (2015):

$$L_{\nu} \propto \begin{cases} \nu^{-0.61}, & h_{\text{p}}\nu < 13.6 \text{ eV}, \\ \nu^{-1.70}, & h_{\text{p}}\nu \ge 13.6 \text{ eV}, \end{cases} \quad (20)$$

which is employed in Madau & Haardt (2015). Also, Khaire et al. (2016) assumed a power-law SED that has different power. The amplitude of the PL model is determined to match its luminosity at a wavelength of 1450 Å to that of the AC model. We use equation (20) up to the X-ray regime (>0.5 keV) for the PL model, as in MH15, but we should mention that broken power-law spectra with power-law index switching to ≈−1 at the X-ray regime are also proposed (e.g. Hasinger et al. 1993). With this type of broken power-law spectra, a larger amount of hard X-ray photons can be produced, similarly to the AC model, and thus a broken power-law spectrum should be regarded as an alternative to the AC model (see e.g. fig. 1 in Laor et al. 1997). At present, it is still unclear which of these spectra, i.e. the AC model or the broken spectra, is suitable for reproducing the SEDs of high-z AGNs. We adopt the AC model for comparison with the simple power-law model used in previous studies, nevertheless, because the AC model and the broken power-law model are expected to yield similar thermal histories if the hard X-ray luminosities are almost identical.

To compare our SED models with observed SEDs, we calculate the weighted average SED as

$$\overline{L}_{\nu}(\nu, z) = \frac{\int \Phi(M_{UV}, z) L_{\nu}(M_{UV}) dM_{UV}}{\int \Phi(M_{UV}, z) dM_{UV}}. \quad (21)$$

In the top panel of Fig. 2, we show the SEDs for the AC/PL models at z = 1, the SED computed by Kawaguchi et al. (2001) and observed composite SED (Laor et al. 1997; Zheng et al. 1997). The figure shows that all SEDs are satisfactorily consistent with each other at E ~ 10⁻² keV, though the composite SED proposed by Zheng et al. (1997) is slightly lower than the others, due to less IGM correction (Lusso et al. 2015). On the other hand, in the X-ray regime (E > 0.5 keV) the AC model shows excess coming from coronal gas compared with the PL model. Additionally, in the bottom panel of Fig. 2, we find that the AC model has significantly higher ionizing photon emissivity at z = 6, due to the dependence of the peak frequency on M<sub>BH</sub> and λ<sub>Edd</sub> shown in equation (18).

### 2.2.3 AGNs emissivity at 13.6 eV

The AGN hydrogen ionizing emissivity ε<sub>912</sub>, namely the emissivity of AGNs at E = 13.6 eV, is a good indicator to show how much the photons emitted from AGNs contribute to cosmic reionization. Thus, we first compare our model with other models in previous studies by using ε<sub>912</sub>. We plot the redshift evolution of ε<sub>912</sub> for different α<sub>hz</sub> in Fig. 3. The corresponding values in MH15 and Khaire & Srianand (2015) are shown by the thick solid and thin solid lines, respectively. The luminosity of AGNs over a wide redshift range is obtained through many independent observations. Using equations (13) and (20), we also calculate ε<sub>912</sub> for different observations. Palanque-Delabrouille et al. (2013) provided g-band luminosity functions in the redshift range 0.68 < z < 4 from SDSS-III and Multiple Mirror Telescope data. We convert the g-band LFs into i-band LFs by using equation (2) in Croom et al. (2009) and the i-band LFs into UVLFs by using equation (B1) in

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Redshift evolution of AGN emissivity at 912 Å. The plotted values are converted from the emissivity at 1450 Å by using equation (20). The dot-dashed, dotted and dashed lines represent the emissivity for our UVLF model with α<sub>hz</sub> = −0.5, −1.5 and −2.2, respectively. For comparison, the model in MH15 is shown by the thick solid line and the emissivity obtained by optical measurements in Khaire & Srianand (2015) by the thin solid line. We also show AGN ionizing emissivity inferred from other observations: Masters et al. (2012) in red small squares, Giallongo et al. (2015) in green large squares and Palanque-Delabrouille et al. (2013) in blue circles.</figcaption></figure>

Lusso et al. (2015):

$$ \begin{align} M_g(z=2) &= M_i - 1.25 \log \left( \frac{4670 \text{ Å}}{7471 \text{ Å}} \right), \nonumber \\ M_i(z=2) &= M_{UV} + 1.28, \tag{22} \end{align} $$

where M<sub>g</sub> is the g-band magnitude and M<sub>i</sub> is the i-band magnitude. It should be mentioned that different power-law indexes are assumed in Croom et al. (2009) and Lusso et al. (2015). However, the difference in M<sub>UV</sub> caused by employing the different power-law indexes is almost negligible (ΔM<sub>UV</sub> = 0.056). We show ε<sub>912</sub> inferred from Palanque-Delabrouille et al. (2013) by blue circles in Fig. 3. The red small squares in Fig. 3 show ε<sub>912</sub> at z > 3 inferred from quasar luminosity functions obtained from the Cosmic Evolution Survey (COSMOS: Masters et al. 2012), by which a characteristic decrease of the quasar number density from z = 3 to 4 is suggested. The green large squares in Fig. 3 are obtained from the UVLFs of faint AGNs at 4 < z < 6.5 found by Giallongo et al. (2015).

As shown in Fig. 3, the AGN emissivity at z ≈ 4 has a large uncertainty; although Giallongo et al. (2015) used only 22 AGN candidates in the CANDLES GOODS–South field as samples, their result is obviously different from that by Masters et al. (2012). This inconsistency is probably caused by the different sample sets for each study. Furthermore, the uncertainty in the AGN emissivity increases at z > 6, since the number of observed high-z AGNs decreases with increasing redshift. Due to the large uncertainty, measurement of the AGN ionizing emissivity cannot put a strong constraint on the evolution factor α<sub>hz</sub>.

Before the comparison with EoR observation data in the next section, we show the AGN emissivity at z < 15 for different α<sub>hz</sub> with the AC and PL models. Fig. 4 shows that the thick lines corresponding to emissivity for the AC model are larger than the thin lines corresponding to emissivity for the PL model. For reference, we also show the AGN-dominant model in MH15 by the dot-dashed line. Fig. 4 also indicates that ε<sub>912</sub> in the extremely
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Redshift evolution of AGN emissivity at 912 Å for AC model (thick lines) and PL model (thin lines) with $\alpha_{hz} = -0.5$ (dot-dashed), $-1.5$ (dotted) and $-2.2$ (dashed). The solid line shows the model in MH15. This figure indicates that $\epsilon_{912}$ in the extremely AGN-abundant model ($\alpha_{hz} = -2.2$) is more than an order of magnitude higher than that in MH15.</figcaption></figure>

AGN-abundant model ($\alpha_{hz} = -2.2$) is more than an order of magnitude higher than that in MH15. This suggests that photoionization rates in the extreme models exceed observationally measured values at $2 < z < 6$ (e.g. Calverley et al. 2011; Becker & Bolton 2013), which are reproduced consistently by MH15. We discuss this point later in Section 4.3.

## 3 RESULTS

Since the AGN abundance has a large uncertainty at $z > 4$, hereafter we treat $\alpha_{hz}$ as a free parameter. The escape fraction $f_{esc}$ is also highly uncertain, because direct measurement of $f_{esc}$ is considerably difficult. Thus, we compare our results with observational data of H I/He II ionizing histories and Thomson optical depth of the CMB in order to obtain constraints on these parameters.

### 3.1 Dependence of the reionization history on $f_{esc}$ and $\alpha_{hz}$

Before comparison with EoR observations, we first show the dependence of the reionization evolution on the parameters $f_{esc}$ and $\alpha_{hz}$. Figs 5 and 6, respectively, show the H I and He II reionization histories for different sets of $f_{esc}$ and $\alpha_{hz}$. In these figures, the reionization histories calculated with the AC and PL models are represented by thick and thin lines, respectively. We will describe the impact of the SED shape on the reionization history in Section 3.3. As we naturally expect, the H I reionization tends to be completed earlier as $f_{esc}$ ($\alpha_{hz}$) increases (decreases), while the He II reionization history has a stronger dependence on $\alpha_{hz}$ than on $f_{esc}$, because of the existence of X-ray photons that can ionize He II more effectively.

The CMB optical depth to Thomson scattering is a good probe for the evolution of the EoR. The optical depth is obtained as (e.g. Planck Collaboration XLVII 2016)

$$ \tau_e(z) = \int_0^z c \sigma_T n_e (1+z)^3 \frac{dt}{dz} dz, \qquad (23) $$

where $\sigma_T$ is the Thomson cross-section. Fig. 7 shows the dependence of $\tau_e$ on the parameter set of $f_{esc}$ and $\alpha_{hz}$. The optical depth $\tau_e$ does not change after H I reionization is completed, because almost all electrons are produced by hydrogen ionization. Therefore, $\tau_e$ is

<figure><img src="image_7.png" /><figcaption>**Figure 5.** H I reionization histories computed with various sets of $f_{esc}$ and $\alpha_{hz}$. The thick lines represent H I reionization histories computed with the AC model and the thin lines represent those with the PL model. Top panel: $\alpha_{hz}$ is fixed to be $-1.0$. The results with $f_{esc} = 0.05, 0.1$ and $0.15$ are shown by solid, dashed and dot-dashed lines, respectively. Bottom panel: $f_{esc}$ is fixed to be 0.1. The results with $\alpha_{hz} = -1.5, -1.0$ and $0.5$ are shown by solid, dashed and dot-dashed lines, respectively. In both panels, the grey vertical line corresponds to $z = 5.7$, which is the completion redshift of the H I reionization inferred from observed QSO spectra (Fan et al. 2006). Two data points with error bars represent the constraints imposed by LAE observations (Ouchi et al. 2010; Konno et al. 2014). The H I reionization history tends to be completed earlier with increasing $f_{esc}$ and decreasing $\alpha_{hz}$. Additionally, for a fixed parameter set of $f_{esc}$ and $\alpha_{hz}$, H I reionization with the AC model is always completed earlier than with the PL model.</figcaption></figure>

sensitive to the hydrogen reionization history and thus to both $f_{esc}$ and $\alpha_{hz}$.

### 3.2 Constraints on the model parameters

There are several observations that provide constraints on the evolution of the H I and He II fractions. From comparison with our numerical results, we can obtain constraints on the parameter set of $f_{esc}$ and $\alpha_{hz}$. We summarize the constrained parameter set for the AC model in Fig. 8. Observations of the spectra of high-$z$ quasars suggest that the cosmic reionization of H I is completed by $z = 5.7$ (Fan et al. 2006). Therefore, we set the criterion for the completion of H I reionization as $f_{H1} < 10^{-4}$ at $z_H = 5.7$. Ouchi et al. (2010) provided a constraint on the neutral fraction of $f_{H1} < 0.2 \pm 0.2$ at $z = 6.6$ from the observation of LAEs. Further, Konno et al. (2014) indicated $0.3 < f_{H1} < 0.8$ at $z = 7.3$ from a drastic decrease in the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** He II reionization histories computed with various sets of $f_{esc}$ and $\alpha_{hz}$. The parameter sets in both panels are the same as those in Fig. 5; the thick lines represent He II reionization histories computed with the AC model and the thin lines represent those with the PL model. Top panel: $\alpha_{hz}$ is fixed to be -1.0. The results with $f_{esc} = 0.05$, 0.1 and 0.15 are shown by solid, dashed and dot-dashed lines, respectively. Bottom panel: $f_{esc}$ is fixed to be 0.1. The results with $\alpha_{hz} = -1.5$, -1.0 and 0.5 are shown by solid, dashed and dot-dashed lines, respectively. Similarly to the H I reionization history, the He II reionization history tends to be completed earlier with increasing $f_{esc}$ and decreasing $\alpha_{hz}$, although the dependence on $f_{esc}$ is significantly weak. This result indicates straightforwardly that AGNs contribute dominantly to the He II reionization. In addition, it turns out that the SED shape of AGNs makes a strong impact on the He II reionization history.</figcaption></figure>

number density of LAEs. Thus, observations of LAEs and QSOs give an upper limit on the H I fraction and set minimum amounts of ionizing sources. As we can see in top panels of Fig. 8, they provide the lower (upper) bound on $f_{esc}$ ($\alpha_{hz}$). On the other hand, the lower bound on the H I fraction at $z = 7.3$ rules out models with large values of $f_{esc}$ and lower values of $\alpha_{hz}$, as shown in the middle left panel of Fig. 8.

Next we focus on He II reionization. The effective optical depth for He II Lyman $\alpha$ absorption suggests that almost all helium in the IGM is fully ionized at $z \sim 2.7$ (Kriss et al. 2001; Shull et al. 2004, 2010; Zheng et al. 2004; Fechner et al. 2006; Furlanetto & Dixon 2010). However, with our LF model based on the evolution of $\Phi_*$ derived by Giallongo et al. (2015), the AGN emissivity is sufficiently large to ionize He II by $z = 2.7$ even if $\alpha_{hz} = 0.0$. Therefore, observations of He II optical depth rule out high-emissivity AGN models (see also Mitra et al. 2016).

<figure><img src="image_9.png" /><figcaption>**Figure 7.** CMB Thomson optical depth with various sets of $f_{esc}$ and $\alpha_{hz}$. The lines and the parameter sets in both panels are the same as those in Fig. 5; the thick lines represent $\tau_e$ computed with the AC model and the thin lines $\tau_e$ computed with the PL model. Top panel: $\alpha_{hz}$ is fixed to be -1.0. The results with $f_{esc} = 0.05$, 0.1 and 0.15 are shown by solid, dashed and dot-dashed lines, respectively. Bottom panel: $f_{esc}$ is fixed to be 0.1. The results with $\alpha_{hz} = -1.5$, -1.0 and 0.5 are shown by solid, dashed and dot-dashed lines, respectively. In each panel, a constraint imposed by the CMB observation, i.e. $\tau_e = 0.058 \pm 0.012$, is shown as a grey region (Planck Collaboration XLVII 2016).</figcaption></figure>

The Planck observation (Planck Collaboration XLVII 2016) measured the CMB optical depth for Thomson scattering and provided $\tau_e = 0.058 \pm 0.012$. In Fig. 7, we plot the lower and upper bounds on $\tau_e$. Since the optical depth depends on $f_{esc}$ and $\alpha_{hz}$, we obtain the constraint on these parameters. We plot the upper and lower bounds on $f_{esc}$ and $\alpha_{hz}$ in the bottom left panel of Fig. 8.

Combining all of the above constraints, we obtain the bottom right panel of Fig. 8, which shows the parameter sets satisfying all constraints. Again, we emphasize that our model is based on the evolution of $\Phi_*$ derived from Giallongo et al. (2015). As we can see, the evolution factor of the faint-end slope of the AGN LF should be $\alpha_{hz} < -1.5$ and the escape fraction must satisfy $f_{esc} < 0.15$. The total constraint is contributed mostly from observations of the H I fraction at $z = 5.7$ and 7.3. Employing the AC model, this result implies that the high AGN abundance at high redshifts reported in Giallongo et al. (2015), i.e. $\alpha_{hz} \sim -1.5$ (see Fig. 3), is not allowed unless $f_{esc}$ is considerably small. The allowed values of $f_{esc}$ are relatively smaller than those predicted in previous works, because AGNs with the AC model have non-negligible contribution to ionization, as shown by Fig. 2.
<figure><img src="image_10.png" /><figcaption>**Figure 10.** Top: evolution of the H II fraction for our constrained models with some typical sets of parameters. The results obtained with ($f_{esc}$, $\alpha_{hz}$) = (0.05, -1.00) and (0.10, 0.00) for the AC model are represented by solid and dot-dashed lines, respectively, while the results obtained with ($f_{esc}$, $\alpha_{hz}$) = (0.03, -1.90) and (0.15, 0.00) for the PL model are shown by dashed and dotted lines, respectively. The result with ($f_{esc}$, $\alpha_{hz}$) = (0.10, 0.00) and $z_{AGN}$ = 10 are also shown by the thin dotted line. Bottom: thermal history for the same models. All of the reionization histories are almost the same, but the thermal histories are obviously different, depending on the abundance of AGNs at high redshifts and the SED type of AGNs. The thermal history is more sensitive to SED type.</figcaption></figure>

an indicator of the existence of AGNs at high redshifts and a probe of the SED shape. Here, note that the thermal history is much more sensitive to the SED shape than $\alpha_{hz}$. Increasing $z_{AGN}$ also provides more efficient heating at higher redshifts. The IGM temperature in the case $z_{AGN}$ = 10 is higher than that of the fiducial model with the same values of $f_{esc}$ and $\alpha_{hz}$, because the larger value of $z_{AGN}$ leads to high AGN emissivity at higher redshifts.

The thermal history is expected to be imprinted on 21-cm line signals from neutral hydrogen through the spin temperature. The 21-cm line signals are a major target of the SKA and its pathfinders, such as the Murchison Widefield Array (MWA) and Low-Frequency Array (LOFAR). In particular, statistical quantities such as the power spectrum, bispectrum, variance and skewness are expected to be able to probe fluctuations in 21-cm brightness temperature (Madau, Meiksin & Rees 1997; Tozzi et al. 2000; Furlanetto, Oh & Briggs 2006; Pritchard & Furlanetto 2007; Baek et al. 2010; Pober et al. 2014; Shimabukuro et al. 2015; Yoshiura et al. 2015).

# 4 DISCUSSION

## 4.1 Model uncertainty

In this subsection, we discuss the uncertainties in our model. First, let us consider the SED and LF of AGNs. We assumed that the SED depends on the Eddington ratio $\lambda_{Edd}$ and all AGNs have the same value for a fixed redshift. In spite of the fact that some previous studies reported dependences of $\lambda_{Edd}$ on AGN luminosity (e.g. Collin & Kawaguchi 2004; Nobuta et al. 2012), we ignored the dependence, since it has large variance. For example, when faint (bright) AGNs have relatively small (large) values of $\lambda_{Edd}$, their ionization efficiency would become smaller (larger). Hence the constraint on $\alpha_{hz}$ would become weaker if we were to take this effect into account. Additionally, we need to mention that there are other studies reporting less abundant AGNs at $z > 3$ than Giallongo et al. (2015) (e.g. Georgakakis et al. 2015; Weigel et al. 2015). If we employ the luminosity functions obtained in their study, the allowed parameter regions in Figs 8 and 9 shift systematically towards the upper right-hand side.

Next, concerning star-forming galaxies, we used a fixed SFRD model. However, the SFRD depends on the faint-end shape and dust correction of UV luminosity functions. Although a recent observation indicated an observational faint-end limit of UV magnitude $M_{UV} = -17$ and a faint-end slope $\alpha = -2$ at $z > 6$ (Ishigaki et al. 2015; Bouwens et al. 2015), UV luminosity functions during the EoR are still less than that of at lower redshift. For instance, Ishigaki et al. (2015) showed that the SFRD integrated down to $M_{UV} = -10$ is larger than that integrated down to $M_{UV} = -17$ by a factor of 3, assuming $\alpha = -2$. Our SFRD model (Madau & Dickinson 2014) was obtained by integration of the luminosity function down to $0.03L^*$, where $L^*$ is the characteristic luminosity of the Schechter function. Further, dust correction causes a variation by a factor of 2 to the SFRD (Hopkins & Beacom 2006). Additionally, a rapid decrease in the SFRD at $z > 8$ is indicated from a recent observation (Bouwens et al. 2015). Thus, the contribution of galaxies to the H I reionization at $z > 8$ may be overestimated in our calculations. For these reasons, the SFRD has an uncertainty of a factor of ~2–6. In our calculation, $f_{esc}$ should be considered to incorporate this uncertainty.

Finally, let us see the effect of inhomogeneity of the IGM. We treated the IGM as being uniform in our calculations, however the observed Universe is obviously inhomogeneous. The effective optical depth $\tau_{eff}$ through an inhomogeneous medium is estimated as an integration of the H I distribution function and the continuum optical depth with respect to H I column density (Haardt & Madau 2012). Inoue et al. (2014) reported a distribution function of intergalactic absorbers (IGA) $F_{IGA}(z) = \partial^2 n/\partial z \partial N_{H I}$ at $z < 6$ and gave an analytic model. Khaire et al. (2016) have shown that an H I distribution function derived from a cosmological hydrodynamic simulation corresponds well to the analytic model by Inoue et al. (2014) at $z = 5-6$ and used it for computing the reionization history. Despite the efforts of these studies, there is still large uncertainty in the H I distribution function during the EoR; we estimate it to be consistent with the reionization process that affects the H I distribution function. Since constructing the complete model is very tough work, we evaluate the impact of inhomogeneity with a simple model. We adopt the analytic model by Inoue et al. (2014) in the post-reionization epoch and extrapolate to higher redshifts as $F_{IGA}(z) = F_{IGA}(z = 5.7) (f_{H I}(1 + f_{H II}^{ion})/f_{H II}^{ion}(1 - f_{H II}^{ion})) ((1+z)/(1+5.7))^2$. Here, $f_{H II}^{ion} = 3 \times 10^{-3}$ corresponds to the mass-weighted ionized
<figure><img src="image_11.png" /><figcaption>Figure 11. H I ionization history for a homogeneous (solid line) and inhomogeneous (dashed line) IGM.</figcaption></figure>

fraction at $z = 5.7$ (Fan et al. 2006). For simplicity, we assume that the distribution of He I and He II follows the H I distribution. Therefore, we estimate the column density of He I and He II as $N_{\text{HeI}} = (f_{\text{HeI}}/(f_{\text{H}} + f_{\text{HII}}^{\text{ion}})) (n_{\text{He}}/n_{\text{H}}) N_{\text{H}}$ and $N_{\text{HeII}} = (f_{\text{HeII}}/(f_{\text{H}} + f_{\text{HII}}^{\text{ion}})) (n_{\text{He}}/n_{\text{H}}) N_{\text{H}}$. As a result, as we can see in Fig. 11, a model with $f_{\text{esc}} = 0.05$ and $\alpha_{\text{hz}} = -1.00$ completes H I reionization at $z_{\text{H}} = 5.58$ in a clumpy Universe, while completion is at $z_{\text{H}} = 6.02$ for a uniform IGM. Therefore, the appropriate value of $f_{\text{esc}}$ ($\alpha_{\text{hz}}$) may become larger (smaller) than our results if we incorporate the inhomogeneity of the IGM.

## 4.2 Remaining parameters

In this work, we have mainly focused on $f_{\text{esc}}$ and $\alpha_{\text{hz}}$. However, our model has two other parameters, $z_{\text{AGN}}$ and the faint end of UV luminosity $M_{\text{UV}}^{\text{min}}$.

First, we discuss $z_{\text{AGN}}$, which determines the critical redshift of AGN abundance decrease at higher redshifts. This parameter relates to many interesting questions on the origin of the first AGNs.

The value of $z_{\text{AGN}}$ would be larger than 6 where AGNs have been observed, but it is highly uncertain. When we set $z_{\text{AGN}} = 10$ rather than our fiducial value $z_{\text{AGN}} = 6$, AGNs ionize H I at higher redshifts. In the left panel of Fig. 12, we show the allowed parameter region for calculations with $z_{\text{AGN}} = 10$. Consequently, we find that the result is not significantly different from the fiducial one.

Secondly, we discuss the faint end of the AGN UV luminosity function. Although a faint AGN with $M_{\text{UV}} = -19$ has been found at $z = 5.75$ in Giallongo et al. (2015), the faint end has not been well determined observationally even at lower redshifts. However, we did not treat $M_{\text{UV}}^{\text{min}}$ as a main parameter in this work, because faint AGNs ($-18 < M_{\text{UV}} < -15$) do not contribute significantly to the total AGN luminosity, in as far as we consider a faint-end slope of $\alpha_{\text{hz}} > -2.0$ (the bottom panel of Fig. 2). Therefore, the results are expected to be rather insensitive to the value of the faint end. This can be confirmed in the right panel of Fig. 12, where the faint end is set to be $M_{\text{UV}}^{\text{min}} = -15$.

## 4.3 H I photoionization rate

As we mentioned in Section 2.2.3, an extreme abundant AGN model provides a higher H I photoionization rate than observationally measured values. There are several measurements at $z \sim 2-6$, such as the Lyman $\alpha$ opacity of the IGM (Becker & Bolton 2013) and the proximity effect (Dall’Aglio, Wisotzki & Worseck 2008, 2009). Here we set a conservative upper limit, $3 \times 10^{-12}$, at $z \sim 2-4$, because these observed values have a large uncertainty and conflict with each other (see e.g. fig. 13 in Calverley et al. 2011). Then, we evaluate the photoionization rate of AGNs as

$$ \Gamma_{\text{H I}}^{\text{AGN}} = \int_{v_{\text{H I}}} dv \frac{F(z, v)}{h_{\text{p}} v} \sigma_{\text{H I}}(v) [s^{-1}]. \quad (24) $$

Thus, we find that this constraint provides a lower limit of $\alpha_{\text{hz}} > -1.6$ for the AC model and $\alpha_{\text{hz}} > -1.9$ for the PL model. This additional constraint hardly changes the lower limit of $\alpha_{\text{hz}}$ for the AC model. However, for the PL model, the constraint rules out AGN-driven reionization. Despite these additional constraints, the existence of abundant faint AGNs at high redshifts, as suggested in Giallongo et al. (2015), is allowed.

In Fig. 13, we show the redshift evolution of $\Gamma_{\text{H I}}^{\text{AGN}}$ obtained from AC models with three different parameter sets satisfying all observational constraints in Fig. 8. We also plot $\Gamma_{\text{H I}}^{\text{AGN}}$ obtained from observations for comparison (Dall’Aglio et al. 2009; Calverley et al. 2011; Wyithe & Bolton 2011; Becker & Bolton 2013). The

<figure><img src="image_12.png" /><figcaption>Figure 12. Constraints on the parameters $f_{\text{esc}}$ and $\alpha_{\text{hz}}$ for two kind of models with $z_{\text{AGN}} = 10$ (left) and $M_{\text{min}} = -15$ (right). The black regions are combined allowed regions.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 13.** H I photoionization rate of AGN for three kinds of parameter sets with the AC model. The results with ($\alpha_{hz}$, $f_{esc}$) = (-1.5, 0.0), (-1.0, 0.07) and (-0.05, 0.10) are shown by a red solid line, blue dashed line and green dot-dashed line. It should be noted that the H I photoionization rate of galaxies is not included in these lines. We plot H I photoionization rates inferred from the literature with $1\sigma$ error: Becker & Bolton (2013) in yellow open squares, Calverley et al. (2011) in cyan open circles, Dall'Aglio et al. (2009) in magenta filled circle and Wyithe & Bolton (2011) in green filled squares.</figcaption></figure>

three lines of our models are consistent with constraints within $3\sigma$ error.

We note that $f_{H1}$ is fixed artificially after reionization in our analysis. However, in practice, it evolves at $z < 6$ and the photoionization rate has a variation proportional to $f_{H1}$. In addition, we do not consider the contribution to the H I photoionization rate from galaxies. Thus, the H I photoionization rate is increased when including the galactic contribution. According to Khaire et al. (2016), the H I photoionization at $2 < z < 4$ requires $f_{esc} < 0.1$, which conflicts with our constraints, except for an abundant AGN model with $\alpha_{hz} < -1$. Thus, a galaxy-dominant reionization model needs a decrease in $f_{esc}$ at $z < 4$.

## 4.4 Cosmic X-ray background

As we showed, high-redshift AGNs can contribute significantly to the reionization of intergalactic hydrogen and helium. Hence there is a possibility that they contribute to the unresolved X-ray background (XRB) as well. In this subsection, we estimate the unresolved XRB coming from AGNs and discuss whether the unresolved XRB can place an additional constraint on the parameter of AGN abundance.

We simply assume the flux from AGNs at $z > 5$ to be the unresolved XRB at the preset epoch. It should be noted that the intensity includes bright AGNs and, in practice, there is a contribution from faint AGNs at $z < 5$. In Fig. 14, we show the dependence of the unresolved XRB at 2 keV on $\alpha_{hz}$ for both AC and PL models. The unresolved XRB at 2 keV is estimated as

$$J[2 \text{ keV}] = \frac{F(z = 0, \nu_{2\text{keV}})}{4\pi} \qquad (25)$$

By comparing the observationally proposed value $J[2 \text{ keV}] = 2.9^{+1.6}_{-1.3} \times 10^{-27} [\text{erg cm}^{-2} \text{ s}^{-1} \text{ Hz}^{-1} \text{ sr}^{-1}]$ (Moretti et al. 2012), we find that $\alpha_{hz} < -2.1$ is prohibited for the AC model and hence the unresolved XRB does not work as an additional constraint on $\alpha_{hz}$ (see Fig. 8).

<figure><img src="image_14.png" /><figcaption>**Figure 14.** Unresolved X-ray background (XRB) spectrum at 2 keV for the AC model (solid line) and PL model (dashed line) as a function of $\alpha_{hz}$. The unresolved XRB derived by Moretti et al. (2012) is indicated by the shaded region. For the AC model, the numerically computed unresolved XRB exceeds the observed value if faint AGNs are extremely abundant by as much as $\alpha_{hz} < -2$, although this constraint does not dominate the other constraints shown in Fig. 8.</figcaption></figure>

## 4.5 Radio-loud quasars

If AGNs were too abundant and bright in the radio frequency range at high redshift, their emission would serve as a background of the 21-cm line of neutral hydrogen. In this subsection we estimate, using our constrained AGN model, the number of photons that radio-loud QSOs would emit in the radio frequency range and compare it with that of the CMB. We base our estimation on the result of Elvis et al. (1994) and Prieto et al. (2010), in which the authors found the SED of radio-loud QSOs in the radio frequency range to be 2–3 orders of magnitude smaller than the SED at $E \approx 10$ eV. Remembering that the intensity ratio between 1.25 µm and 1 keV is approximately $10^{-0.4}$, it is expected that the intensity ratio between 21 cm and 1 keV is about $\sim 10^{-2}$. Setting the ratio as $10^{-\alpha}$ ($\alpha = 2$ is the fiducial value), we write the SED at the 21-cm line frequency as

$$\frac{F_{\nu_{21}}}{4\pi} = \frac{10^{\alpha}}{4\pi} \left( \frac{\nu_{1\text{keV}}}{\nu_{21\text{cm}}} \right) F_{1\text{keV}} \qquad (26)$$

The conversion formula from intensity to brightness temperature is given by

$$\left(\frac{T_b}{K}\right) = 3.255 \times 10^{-5} \left(\frac{\nu}{\text{GHz}}\right)^2 \left(\frac{F_\nu/4\pi}{\text{Jy sr}^{-1}}\right) \qquad (27)$$

We evaluate $F_{\nu_{1\text{keV}}}/4\pi = 1.79 \times 10^{-3} \text{ Jy sr}^{-1}$ at $z = 6$ for the AC model with $\alpha_{hz} = -1.5$ and $T_b < 10 \times 10^{-\alpha} [\text{K}]$ in terms of brightness temperature. If we assume that AGNs are all radio-loud QSOs ($\alpha = 2$), this estimate gives $T_b < 0.1 \text{ K}$. Therefore, the contribution of radio-loud AGNs to the background of the 21-cm line should be negligible compared with the CMB ($T_\nu = 2.7(1+z)$).

# 5 CONCLUSION

In this work, we studied H I reionization with a simple model of star-forming galaxies and AGNs. We then employed a theoretically motivated SED model and a power-law SED model. We compared observational data on the H I fraction and CMB optical depth with our model and obtained constraints on the escape fraction of ionizing photons, $f_{esc}$, and the faint-end slope of the AGN luminosity
function at $z > 4.25$, $\alpha_{hz}$. As a result, we have found that models with $\alpha_{hz} > -1.5$ and $f_{esc} < 0.15$ are consistent with existing observations and constraints come mostly from the H1 fraction history. Our result suggests that an AGN-dominated model with AGN abundance as large as the estimation in Giallongo et al. (2015) is allowed, while a galaxy-dominated model is also allowed. It should be noted that we model the redshift evolution of the AGN luminosity function so as to let $M_*$ and $\Phi_*$ be consistent with the values in Giallongo et al. (2015). This keeps AGNs to ionize the Universe even with $\alpha_{hz} = 0$ and requires lower $f_{esc}$. We have also calculated the thermal evolution of the IGM at $z > 6$ and found that it is sensitive to the abundance of AGNs and the shape of the SED. Therefore, forthcoming 21-cm line observations will be useful to obtain more powerful constraints on properties, i.e. the abundance and SED, of ionizing sources. Recently, Kulkarni et al. (2017) have shown that the AGN-dominant model enhances the 21-cm power spectrum on large scales.

We have discussed the H1 photoionization rate of AGNs and the contribution from AGNs to the unresolved XRB and found that extremely abundant AGN models conflict with the observed H1 photoionization rate and unresolved XRB measurements. Furthermore, even if we assume that radio-loud AGNs exist at high redshifts, the radiation is too small to serve as background emission of the 21-cm line signal of neutral hydrogen.

## ACKNOWLEDGEMENTS

We are grateful to A. K. Inoue for his suggestions regarding model of AGNs, to T. Kawaguchi for kindly providing us with the theoretical and observational SED data and to H. Shirakata for providing simulation results. We also thank the members of SKA-JP EoR science working group for fruitful discussion and comments. This work is supported by Grant-in-Aid from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, Nos 26610048, 15H05896, 16H05999 (KT), 24340048 (KT, KI), 25-3015 (HS), 15K17646 (HT), 16J01585 (SY) and a grant from NAOJ. HT also acknowledges support by the MEXT Program for Leading Graduate Schools PhD professional, 'Gateway to Success in Frontier Asia'.

## REFERENCES

Atek H. et al., 2015, ApJ, 800, 81

Baek S., Semelin B., Di Matteo P., Revaz Y., Combes F., 2010, A&A, 523, A4

Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

Becker G. D., Bolton J. S., Haehnelt M. G., Sargent W. L. W., 2011, MNRAS, 410, 1096

Bouwens R. J. et al., 2015, ApJ, 803, 34

Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543

Choudhury T. R., Ferrara A., 2005, MNRAS, 361, 577

Choudhury T. R., Puchwein E., Haehnelt M. G., Bolton J. S., 2015, MNRAS, 452, 261

Collin S., Kawaguchi T., 2004, A&A, 426, 797

Croom S. M. et al., 2009, MNRAS, 399, 1755

Dall'Aglio A., Wisotzki L., Worseck G., 2008, A&A, 491, 465

Dall'Aglio A., Wisotzki L., Worseck G., 2009, preprint (arXiv:0906.1484)

Dunlop J. S. et al., 2013, MNRAS, 432, 3520

Elvis M. et al., 1994, ApJS, 95, 1

Fan X. et al., 2006, AJ, 132, 117

Fechner C. et al., 2006, A&A, 455, 91

Frontera F. et al., 2007, ApJ, 666, 86

Fukugita M., Kawasaki M., 1994, MNRAS, 269, 563

Furlanetto S. R., Dixon K. L., 2010, ApJ, 714, 355

Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Phys. Rep., 433, 181

Georgakakis A. et al., 2015, MNRAS, 453, 1946

Giallongo E. et al., 2015, A&A, 578, A83

Glikman E., Djorgovski S. G., Stern D., Dey A., Jannuzi B. T., Lee K.-S., 2011, ApJ, 728, L26

Gnedin N. Y., Kravtsov A. V., Chen H.-W., 2008, ApJ, 672, 765

Haardt F., Madau P., 2012, ApJ, 746, 125

Hao J.-M., Yuan Y.-F., Wang L., 2015, MNRAS, 451, 1875

Hasinger G., Burg R., Giacconi R., Hartner G., Schmidt M., Trumper J., Zamorani G., 1993, A&A, 275, 1

Hopkins A. M., Beacom J. F., 2006, ApJ, 651, 142

Hu E. M., Cowie L. L., Barger A. J., Capak P., Kakazu Y., Trouille L., 2010, ApJ, 725, 394

Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625

Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2015, ApJ, 799, 12

Iwata I. et al., 2009, ApJ, 692, 1287

Jiang L. et al., 2008, AJ, 135, 1057

Kashikawa N., 2011, ApJ, 734, 119

Kawaguchi T., Shimura T., Mineshige S., 2001, ApJ, 546, 966

Khaire V., Srianand R., 2015, MNRAS, 451, L30

Khaire V., Srianand R., Choudhury T. R., Gaikwad P., 2016, MNRAS, 457, 4051

Kimm T., Cen R., 2014, ApJ, 788, 121

Konno A. et al., 2014, ApJ, 797, 16

Kriss G. A. et al., 2001, Science, 293, 1112

Kulkarni G., Choudhury T. R., Puchwein E., Haehnelt M. G., 2017, MNRAS, 469, 4283

Laor A., Fiore F., Elvis M., Wilkes B. J., McDowell J. C., 1997, ApJ, 477, 93

Lusso E., Worseck G., Hennawi J. F., Prochaska J. X., Vignali C., Stern J., O'Meara J. M., 2015, MNRAS, 449, 4204

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Madau P., Haardt F., 2015, ApJ, 813, L8

Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429

Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648

Masters D. et al., 2012, ApJ, 755, 169

Matthee J., Sobral D., Best P., Khostovan A. A., Oteo I., Bouwens R., Röttgering H., 2017, MNRAS, 465, 3637

Mitra S., Choudhury T. R., Ferrara A., 2016, preprint (arXiv:1606.02719)

Moretti A., Vattakunnel S., Tozzi P., Salvaterra R., Severgnini P., Fugazza D., Haardt F., Gilli R., 2012, A&A, 548, A87

Mostardi R. E., Shapley A. E., Steidel C. C., Trainor R. F., Reddy N. A., Siana B., 2015, ApJ, 810, 107

Naidu R. P. et al., 2016, preprint (arXiv:1611.07038)

Nobuta K. et al., 2012, ApJ, 761, 143

Oesch P. A. et al., 2014, ApJ, 786, 108

Ota K. et al., 2010, ApJ, 722, 803

Ouchi M. et al., 2009, ApJ, 706, 1136

Ouchi M. et al., 2010, ApJ, 723, 869

Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2009, MNRAS, 429, L94

Palanque-Delabrouille N. et al., 2013, A&A, 551, A29

Pawlik A. H., Schaye J., van Scherpenzeel E., 2009, MNRAS, 394, 1812

Planck Collaboration XXIII, 2015, A&A, 594, A13

Planck Collaboration XLVII, 2016, A&A, 596, A108

Pober J. C. et al., 2014, ApJ, 782, 66

Price L. C., Trac H., Cen R., 2016, preprint (arXiv:1605.03970)

Prieto M. A., Reunanen J., Tristram K. R. W., Neumayer N., Fernandez-Ontiveros J. A., Orienti M., Meisenheimer K., 2010, MNRAS, 402, 724

Pritchard J. R., Furlanetto S. R., 2007, MNRAS, 376, 1680

Raicević M., Theuns T., 2011, MNRAS, 412, L16

Razoumov A. O., Sommer-Larsen J., 2010, ApJ, 710, 1239

Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19

Shen Y., Kelly B. C., 2012, ApJ, 746, 169
Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., Ichiki K., 2015, MNRAS, 451, 467

Shull J. M., Tumlinson J., Giroux M. L., Kriss G. A., Reimers D., 2004, ApJ, 600, 570

Shull J. M., France K., Danforth C. W., Smith B., Tumlinson J., 2010, ApJ, 722, 1312

Siana B. et al., 2015, ApJ, 804, 17

Steidel C. C., Pettini M., Adelberger K. L., 2001, ApJ, 546, 665

Susa H., 2006, PASJ, 58, 445

Syphers D., Shull J. M., 2014, ApJ, 784, 42

Tozzi P., Madau P., Meiksin A., Rees M. J., 2000, ApJ, 528, 597

Ueda Y., Akiyama M., Hasinger G., Miyaji T., Watson M. G., 2014, ApJ, 786, 104

Umemura M., Susa H., Hasegawa K., Suwa T., Semelin B., 2012, Prog. Theor. Exp. Phys., 1, 01A306

Weigel A. K., Schawinski K., Treister E., Urry C. M., Koss M., Trakhtenbrot B., 2015, MNRAS, 448, 3167

Wise J. H., Demchenko V. G., Halicek M. T., Norman M. L. Turk M. J., Abel T., Smith B. D., 2014, MNRAS, 442, 2560

Worseck G., Prochaska J. X., Hennawi J. F., McQuinn M., 2014a, preprint (arXiv:1405.7405)

Worseck G. et al., 2014b, MNRAS, 445, 1745

Wyithe J. S. B., Bolton J. S., 2011, MNRAS, 412, 1926

Yajima H., Choi J.-H., Nagamine K., 2011, MNRAS, 412, 411

Yoshiura S., Shimabukuro H., Takahashi K., Momose R., Nakanishi H., Imai H., 2015, MNRAS, 451, 266

Zheng W., Kriss G. A., Telfer R. C., Grimes J. P., Davidsen A. F., 1997, ApJ, 475, 469

Zheng W. et al., 2004, ApJ, 605, 631

This paper has been typeset from a *TEX/LATEX* file prepared by the author.