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

HAL
open science

# Analysing the 21 cm signal from the epoch of reionization with artificial neural networks

Hayato Shimabukuro, Benoit Semelin

► To cite this version:

Hayato Shimabukuro, Benoit Semelin. Analysing the 21 cm signal from the epoch of reionization with artificial neural networks. Monthly Notices of the Royal Astronomical Society, 2017, 468 (4), pp.3869-3877. 10.1093/mnras/stx734. hal-01554641

HAL Id: hal-01554641

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

Submitted on 27 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_2.png" /></figure>

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

# Analysing the 21 cm signal from the epoch of reionization with artificial neural networks

Hayato Shimabukuro* and Benoit Semelin*

Sorbonne Universités, UPMC, LERMA, Observatoire de Paris, PSL research university, CNRS, F-75014 Paris, France

Accepted 2017 March 22. Received 2017 March 22; in original form 2017 January 24

## ABSTRACT

The 21 cm signal from the epoch of reionization should be observed within the next decade. While a simple statistical detection is expected with Square Kilometre Array (SKA) pathfinders, the SKA will hopefully produce a full 3D mapping of the signal. To extract from the observed data constraints on the parameters describing the underlying astrophysical processes, inversion methods must be developed. For example, the Markov Chain Monte Carlo method has been successfully applied. Here, we test another possible inversion method: artificial neural networks (ANNs). We produce a training set that consists of 70 individual samples. Each sample is made of the 21 cm power spectrum at different redshifts produced with the 21cmFast code plus the value of three parameters used in the seminumerical simulations that describe astrophysical processes. Using this set, we train the network to minimize the error between the parameter values it produces as an output and the true values. We explore the impact of the architecture of the network on the quality of the training. Then we test the trained network on the new set of 54 test samples with different values of the parameters. We find that the quality of the parameter reconstruction depends on the sensitivity of the power spectrum to the different parameters at a given redshift, that including thermal noise and sample variance decreases the quality of the reconstruction and that using the power spectrum at several redshifts as an input to the ANN improves the quality of the reconstruction. We conclude that ANNs are a viable inversion method whose main strength is that they require a sparse exploration of the parameter space and thus should be usable with full numerical simulations.

**Key words:** intergalactic medium – dark ages, reionization, first stars – cosmology: theory.

## 1 INTRODUCTION

According to the standard hierarchical structure formation theory based on the ΛCDM cosmology model, density perturbations grow into small haloes where the first stars and galaxies are formed (e.g. Abel, Bryan & Norman 2002; Bromm, Coppi & Larson 2002; Yoshida et al. 2006; Bromm 2013). These astrophysical sources produce ultraviolet ionizing photons that escape into the intergalactic medium (IGM), which is eventually reionized. The epoch when this process occurs is called the ‘epoch of reionization’ (EoR). The detailed astrophysics of the EoR is currently poorly understood because no observation has yet probed the early and middle stages of the EoR although current observations have provided fruitful information on late stage of the EoR. For example, the absorption spectra of high-z quasars indicate that reionization was complete by $z \sim 6$ (e.g. Fan, Carilli & Keating 2006) and the number density of Ly α emitter galaxies at $z = 6–7$ implies that the neutral hydrogen

fraction increases at $z > 6$ (e.g. Konno et al. 2014). On cosmological scales, reionization induces Thomson scattering of CMB photons off free electrons. The optical depth of Thomson scattering measured by Planck is $\tau_e = 0.058 \pm 0.012$, which corresponds, using an instantaneous reionization toy model, to a reionization at $z = 8.8 \pm 0.9$ (Planck Collaboration XLVII 2016b).

To further improve our understanding of the EoR, observations targeted on the cosmological redshifted 21 cm signal from the EoR are ongoing. The 21 cm signal emission is due to the hyperfine structure of neutral hydrogen atoms and is expected to be a powerful tool to probe the neutral IGM, yielding both astrophysical and cosmological information such as matter density fluctuations and the ionization state and thermal history of the IGM at high redshift (Furlanetto et al. 2006; Pritchard & Loeb 2012).

Recently, some first-generation radio interferometers have been attempting to detect statistically the 21 cm signal from the EoR, such as the Murchison Wide field Array (Tingay et al. 2013), the LOw Frequency ARray (LOFAR) (Rottgering 2003) and the Precision Array for Probing the Epoch of Reionization (PAPER) (Jacobs et al. 2015). These observational efforts have resulted in upper limits of the 21 cm power spectrum (PS) and on the ionized state

<footnote>* E-mail: hayato.shimabukuro@obspm.fr (HS); benoit.semelin@obspm.fr (BS)</footnote>
of the IGM at $z \sim 8$ (e.g. Parsons et al. 2014; Ali et al. 2015; Pober et al. 2015; Beardsley et al. 2016; Greig, Mesinger & Pober 2016). Furthermore, future instruments such as the Square Kilometre Array (SKA) (Mellema et al. 2013) and Hydrogen Epoch of Reionization Array (HERA) (DeBoer et al. 2017) are designed to detect the 21 cm signal PS with higher signal-to-noise ratio and at higher redshift, during the cosmic dawn. The SKA should also be able to image the signal in 3D, which requires high sensitivity on small enough scale (that is sufficient collecting area in a large enough core).

So we expect a wealth of 21 cm data in the near future. Then we face the fundamental question: what we can learn from the data? From theoretical and numerical works, we have already some insights on the process of reionization. For instance, considering galaxies in relatively massive host haloes results in larger and more uniform ionized bubbles and imprints a peak at larger scales in the 21 cm PS (e.g. McQuinn et al. 2007; Iliev et al. 2012). On the other hand, abundant and small minihaloes, which serve as absorption systems by self-shielding from ionizing photons, result in small and disjointed ionized regions (e.g. Ciardi et al. 2006). To extract information on the EoR from the observed 21 cm signal, we need to be able to compute the 21 cm signal from the basic physics of reionization (e.g. Semelin & Iliev 2015). The most self-consistent method is to run numerical simulations. There are currently two approaches. The first is to implement the full radiation hydrodynamics (RHD). This type of simulations is relevant if small scales, where feedback effects such as photoheating play an important role, are resolved (e.g. Ciardi et al. 2000; Ricotti, Gnedin & Shull 2002; Hasegawa & Semelin 2013). The second approach is to run large-scale simulations where radiative transfer (RT) is computed in post-processing (e.g. Iliev et al. 2006; Trac & Cen 2007; Iliev et al. 2014). This type of simulations are able to account for large-scale fluctuations in the signal but they require sub-grid modelling to treat processes on galaxy scales (feedback, ionizing photons escape fraction etc.).

Although numerical study based on simulations is the most consistent method, it implies a large computational cost. This disadvantage is somewhat alleviated by using a seminumerical approach instead of RHD or RT simulations. In most seminumerical simulations, the production of ionizing photons is calculated based on the excursion set formalism, using analytically derived halo mass functions, and density fluctuations are obtained with linear perturbation theory (Mesinger & Furnaletto 2007; Santos et al. 2010; Mesinger et al. 2011). However, the growth of H II regions is simply evaluated by considering the balance between the production of ionizing photons and the number of neutral hydrogen atoms (Furlanetto, Zaldarriaga & Hernquist 2004). Recently, the effect of recombination is also taken into account (e.g. Sobacchi & Mesinger. 2014). The results produced by the seminumerical approach show good agreement with those by RT simulations on a large scale (≥ Mpc) (Zahn et al. 2007, 2011; Mesinger et al. 2011).

To maximize the scientific return of the upcoming observations, it is important to establish systematic procedures to derive constraints on the EoR modelling parameters from the observed 21-cm data. Recently, several works have studied how such constraints can be obtained by exploring the EoR parameter space with seminumerical simulations. For example, Fisher analysis (Pober et al. 2014; Kubota et al. 2016; Liu & Parsons 2016; Shimabukuro et al. 2017) and Bayesian parameter inference such as the Markov Chain Monte Carlo (MCMC) approach (Greig & Mesinger 2015; Bernardi et al. 2016) have been applied.

In this work, we suggest a new approach for parameter reconstruction based on a machine learning method. Machine learning is one of the hot topics in data science as a method to deal with big

data. It is currently applied to many fields such as pattern recognition or search engine. The main purpose of machine learning is to find approximate functions that, given the input, produce the desired outputs. This is achieved by 'learning' from training data sets with known inputs and outputs (e.g. Ball & Brunner 2010; Bloom & Richards 2012).

Recently, machine learning methods have been applied in the field of astronomy. For example, learning from huge galaxy image catalogues helps with morphological classification of galaxies (Folkes, Lahav & Maddox 1996; Lahav et al. 1996; Banerji et al. 2010). Machine learning can also help with the selection and classification of transients (Morii et al. 2016; Mukund et al. 2016). Applying machine learning to a large sample of spectroscopic and photometric galaxies data can improve the accuracy of estimates of photometric redshifts (Collister & Lahav 2004; Vanzella et al. 2004; Samui & Samui Pal 2017). Using simulated gravitational wave data with noise as a learning sample, machine learning can help us search for gravitational wave signals from noisy real data (Kim et al. 2015). In the context of cosmology, machine learning is used to model galaxy formation (Kamdar, Turk & Brunner 2016a,b) or to make templates of non-linear matter PS (Agarwal et al. 2012, 2014). Closer to our field of interest, there is a study applying a machine learning method to estimate the escape fraction of ionizing photons during the EoR (Jensen et al. 2016). In this study, they show how machine learning can estimate the Lyman continuum escape fraction by using mock spectroscopic simulation data. In our work, using a simple astrophysical parametrization of the EoR, we apply artificial neural network (ANN), which is one of the machine learning methods, to reconstruct the parameter values from the 21 cm PS data.

This paper is organized as follows. In Section 2, we introduce the cosmological 21 cm signal and the EoR parametrization we focus on. In Section 3, we describe our ANN and test the impact of its chosen architecture. In Section 4, we show our main results, and we give a summary and discussion in Section 5. Throughout this paper, we employ the best-fitting cosmological parameters obtained by Planck Collaboration XIII (2016a).

# 2 COSMOLOGICAL 21 CM SIGNAL

## 2.1 Introduction of 21 cm PS

The brightness temperature for the 21 cm signal is given by (e.g. Mellema et al. 2013)

$$ \delta T_b(v) = \frac{T_S - T_\gamma}{1+z} (1 - e^{-\tau_{v0}}) \\ \sim 27x_H(1 + \delta_m) \left( \frac{H}{dv_r/dr + H} \right) \left( 1 - \frac{T_\gamma}{T_S} \right) \\ \times \left( \frac{1+z \ \frac{0.15}{10} \ \Omega_b h^2}{\Omega_m h^2} \right)^{1/2} \left( \frac{\Omega_b h^2}{0.023} \right) \left( \frac{\Omega_b h}{0.031} \right) [\text{mK}]. \quad (1) $$

Here, $T_S$ and $T_\gamma$ represent the local spin temperature of the IGM and the CMB temperature, respectively. $\tau_{v0}$ is the local optical depth in the 21 cm rest-frame frequency $v_0 = 1420.4$ MHz, $x_H$ is the local neutral fraction of the hydrogen gas, $\delta_m(x, z) \equiv \rho/\bar{\rho} - 1$ is the evolved matter overdensity, $dv_r/dr$ is the local gradient of the gas velocity along the line of sight and $H(z)$ is the Hubble parameter. All quantities are evaluated at redshift $z = v_0/v - 1$. As we can see, the 21 cm signal includes both astrophysical and cosmological information. Thus, we can hope to use the 21 cm signal to disentangle and quantify them (Pritchard et al. 2015).
Let us now introduce the PS of the 21 cm fluctuations. We define the 21 cm PS as

$$ \langle \delta T_b(\mathbf{k}) \delta T_b(\mathbf{k}') \rangle = (2\pi)^3 \delta(\mathbf{k} + \mathbf{k}') P_{21}(\mathbf{k}). \quad (2) $$

In our context, we use the *dimensional* 21 cm power spectrum, $k^3 P(k)/2\pi^2$. The 21 cm PS would describe the statistical properties of the 21 cm fluctuations perfectly if they were a Gaussian random field. However, the 21 cm fluctuations are expected to deviate from a Gaussian behaviour because of astrophysical effects such as ionization and X-ray heating. Thus, it is useful to compute higher order statistics and one-point statistics such as the bispectrum, the variance and the skewness to estimate non-Gaussian features in the 21 cm fluctuations (e.g. Barkana & Loeb 2008; Harker et al. 2009; Baek et al. 2010; Patil et al. 2014; Shimabukuro et al. 2015, 2016).

In order to generate the 21 cm PS for a given set of astrophysical parameters, we use the publicly available code **21cmFAST** (Mesinger et al. 2011). This code is based on a seminumerical model of cosmic reionization and thermal history of the IGM. It quickly generates maps and PS of the brightness temperature, matter density, velocity, spin temperature and ionization fraction at designated redshifts.

We performed simulations in a 200 Mpc³ comoving box with 300³ grid cells for a wide range of EoR parameters described in the next section. In our calculation, we use the 21 cm PS in the range $0.06\text{Mpc}^{-1} \le k \le 1.4\text{Mpc}^{-1}$ divided into 14 bins.

## 2.2 EoR model parameters

It is common to characterize EoR models with parameters and then examine the effect of changing the parameters on the 21 cm signal. We employ three key parameters that are often used. Let us briefly define these three parameters.

1. $\zeta$, the ionizing efficiency. $\zeta$ is the combination of several parameters related to ionizing photons escaping from high-redshift galaxies and is defined as $\zeta = f_{\text{esc}} f_{*} N_{\gamma} / (1 + \bar{n}_{\text{rec}})$ (Furlanetto et al. 2004, 2006). Here, $f_{\text{esc}}$ is the fraction of ionizing photons escaping from galaxies into the IGM and $f_{*}$ is the fraction of baryons locked into stars. These parameters are extremely uncertain at high redshift (Gnedin et al. 2008; Wise & Cen 2009). $N_{\gamma}$ is the number of ionizing photons produced per baryon in stars and $\bar{n}_{\text{rec}}$ is the mean recombination rate per baryon. In our calculation, we explore a range of $10 \le \zeta \le 60$.

2. $T_{\text{vir}}$, the minimum virial temperature of haloes producing ionizing photons. $T_{\text{vir}}$ parametrizes the minimum mass of haloes producing ionizing photons during the EoR. Typically, $T_{\text{vir}}$ is chosen to be $10^4$ K, corresponding to the temperature above which atomic cooling becomes effective. $T_{\text{vir}}$ parametrizes the physics of star formation in high-redshift galaxies. In haloes with virial temperature > $10^4$ K, atomic cooling is sufficient to trigger gravothermal instability and thus star formation. However, star formation is quenched if AGN or supernovae feedback is effective and the IGM is heated up. This leads to effective minimum virial temperature larger than $10^4$ [K]. In haloes with virial temperature < $10^4$ K, hydrogen molecule cooling is necessary. However, if stars begin to form in a halo, radiative feedback such as the photodissociation of H₂ by Lyman-Werner photons may become effective and prevent the gas from cooling (Fialkov et al. 2013; Sobacchi & Mesinger. 2013). Conversely, positive feedback, such as the enhancement of H₂ molecules formation due to an increase in the free electrons density, tends to push the minimum virial temperature to lower value because cooling becomes more effective. Thus $T_{\text{vir}}$ parametrizes the uncertainty in the efficiency of

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Typical architecture of an ANN. The architecture of the ANN consists of an input layer, a hidden and an output layer of neurons. Each neuron connects the neurons in the next layer.</figcaption></figure>

radiative feedback. In our work, we explore $T_{\text{vir}}$ ranging from $10^3$ K to $5 \times 10^5$ K.

3. $R_{\text{mfp}}$, the mean free path of ionizing photons. The propagation of ionizing photons through the ionized IGM strongly depends on the presence of absorption systems and the sizes of ionized regions are determined by the balance between sinks and sources of ionizing photons (e.g. McQuinn, Oh & Faucher-Giguère 2011). This process is modelled by the maximum mean free path of ionizing photons, $R_{\text{mfp}}$ (Sobacchi & Mesinger. 2014). Physically, the mean free path of ionizing photons corresponds to the typical distance travelled by photons within ionized regions before they are absorbed and is determined by the number density and the optical depth of Lyman-limit systems. In our calculation, we explore $R_{\text{mfp}}$ from 10 to 60 Mpc.

# 3 ARTIFICIAL NEURAL NETWORKS

## 3.1 Architecture of ANN

In this section, we introduce ANN. ANNs are one of the machine learning methods and are a mathematical model inspired by the natural neuron network in our brain. The main purpose of ANNs is to construct approximate functions that associate input data with output data. In order to construct such a function, the ANN has to learn from 'training data'. The architecture of a simple class of ANN consists of three layers: the input layer, the hidden layer and the output layer. Each of them has a number of neurons as shown in Fig. 1. In a more general case, we could choose the number of hidden layers and the number of neurons at each layer arbitrarily. In our study, we use one hidden layer. Note that it is mathematically proven that neural networks with only one hidden layer can approximate any function with any accuracy if we use a large enough number of neurons (Cybenko. 1989; Hornik et al. 1989).

Let us briefly describe the architecture of our ANN. The input data $x_j$ is fed to the jth neurons in the input layer. Each neuron in the input layer is connected to the ith neuron in the hidden layer and a weight $w_{ij}^{(1)}$ is associated with the connection. The input to the ith neuron in the hidden layer $s_i$ is a linear combination of all the input neurons with weight $w_{ij}^{(1)}$:

$$ s_i = \sum_{j=1}^{n} w_{ij}^{(1)} x_j. \quad (3) $$
Here, n is the number of input data. In the hidden layer, the ith neuron is activated by an activation function φ such as its output is $t_i = \phi(s_i)$. Generally, the activation function is a non-linear function. We use the sigmoid function $\phi(s) = 1/(1 + e^{-s})$. The properties of the sigmoid function are such that it saturates and returns a constant output when the absolute value of the input is large and that it is a smooth and differentiable function. Thanks to the non-linear activation function, a trained ANN can express any function.

In the output layer, we compute linear combinations of the activated outputs of the neurons in the hidden layer with weights $w_{ij}^{(2)}$ and obtain the output vector:

$$y_i = \sum_j^k w_{ij}^{(2)} t_j \quad (4)$$

Here, k is the number of neurons in the hidden layer. Note that we do not activate the output value. The aim of training the ANN is to find a set of weights that ensures that the output vectors produced by the ANN for a set of input vectors are sufficiently close to the desired output vectors. Once we adjust the weights to reach this goal on a training sample, we can make predictions for output vectors for arbitrary input vectors outside of the training sample (e.g. new observational data).

A popular algorithm to compute the trained weights is the 'Back propagation algorithm' (Rumelhart et al. 1986). We will describe this algorithm briefly in the following section.

## 3.2 Back propagation algorithm

In this section, we present the back propagation algorithm for the one hidden layer case. We also show the back propagation algorithm for multiple hidden layers case in the Appendix.

In order to quantify how well the output obtained by the ANN approximates the desired output for the training data set, we define the (total) cost function as

$$E = \sum_{n=1}^{N_{\text{train}}} E_n = \sum_{n=1}^{N_{\text{train}}} \left[ \frac{1}{2} \sum_{i=1}^{m} (y_{i,n} - d_{i,n})^2 \right], \quad (5)$$

where $N_{\text{train}}$ is the number of training input vectors and m is the number of neurons at the output layer. y and d are outputs of the ANN and the (desired) training output data, respectively. Our purpose is to find the weight set that minimises the cost function. In order to find this weights set, we need to compute the partial derivative of E with respect to the individual weights $w_{ij}^{(l)}$ ($l = 1, 2$) and find the local minimum of E using gradient descent.

The weights are updated by gradient descent following the formula:

$$\Delta w_{ij}^{(l)} = -\eta \frac{\partial E}{\partial w_{ij}^{(l)}} = -\eta \sum_{n=1}^{N_{\text{train}}} \frac{\partial E_n}{\partial w_{ij}^{(l)}} \quad (6)$$

Here, η is a learning coefficient that controls how fast the weights are updated. We used η = 0.01. We only need to calculate the derivative of the cost function for each training input vector and then sum over all input vectors as shown in equation (6). First, let us consider the derivative with respect to the weights between output layer and hidden layer. In this case ($l = 2$), we can simply calculate the derivative of E as

$$\begin{aligned} \frac{\partial E_n}{\partial w_{ij}^{(2)}} &= \frac{\partial E_n}{\partial y_{i,n}} \frac{\partial y_{i,n}}{\partial w_{ij}^{(2)}} \\ &= (y_{i,n} - d_{i,n})t_j \\ &= (y_{i,n} - d_{i,n})\phi(s_j). \end{aligned} \quad (7)$$

Next, we calculate the derivative of E with respect to the weights between the hidden layer and the input layer. In this case ($l=1$), the derivative of E is

$$\begin{aligned} \frac{\partial E_n}{\partial w_{ij}^{(1)}} &= \frac{\partial E_n}{\partial s_i} \frac{\partial s_i}{\partial w_{ij}^{(1)}} \\ &= \left( \sum_{p=1}^{m} \frac{\partial E_n}{\partial y_p} \frac{\partial y_p}{\partial s_i} \right) x_j \\ &= \left( \sum_{p=1}^{m} (y_{p,n} - d_{p,n}) w_{pi}^{(2)} \phi'(s_i) \right) x_j \end{aligned} \quad (8)$$

In the second line, we use the chain rule for derivative because E depends on the activated neuron $\phi(s)$ in the hidden layer only through the output neuron y. Here, $\phi'$ denotes the derivative of the activation function with respect to s. Using equations (6)–(8), we can iterate on the gradient descent until the outputs obtained by the ANN converge to the desired outputs (minimum of the cost function).

The back-propagation algorithm can be summarized as follows.

(i) Starting with random weights, compute the output of the ANN using equations (3) and (4) for all input vectors in the training set (forward propagation).

(ii) Compute the cost function.

(iii) Compute the derivative of the cost function with respect to the weights between output layer and hidden layer with equation (7) and then the derivative with respect to the weights between input layer and hidden layer with equation (8) (back-propagation).

(iv) Update the weights with equation (6).

(v) Go back to (i) and iterate until the cost function converges to a minimum.

## 3.3 Training data sets and architecture of the ANN in our work

In our case, we prepared 70 training data sets. Each set consists of the 21 cm PS P(k) obtained with 21cmFast as input data and the corresponding EoR parameters used in the simulation, $\theta(= R_{\text{mfp}}, T_{\text{vir}}, \zeta)$ as output data. The architecture of the ANN is the following: (i) 14 neurons in the input layer, (ii) 14 neurons in the hidden layer and (iii) three neurons in the output layer. As we mentioned in Section 2.1, we use 14 bins for the 21 cm PS, and we use three EoR parameters. This is why the numbers of neurons in the input layer and in the output layers are 14 and 3, respectively. Note that the number of redshifts used to train will always match the number of redshifts being fit.

To analyse in details how the ANN performs we look at partial cost functions: we use the normalized root-mean-square error (RMSE) defined by

$$\text{RMSE} = \sqrt{\frac{1}{N_{\text{train}}} \sum_{i=1}^{N_{\text{train}}} X^2} \quad (9)$$

Here, $X = (\theta_{\text{ANN}} - \theta_{\text{data}})/\theta_{\text{data}}$ (with $\theta$ one of $R_{\text{mfp}}, T_{\text{vir}}$ or, $\zeta$). $\theta_{\text{ANN}}$ and $\theta_{\text{data}}$ are the EoR parameters evaluated by the ANN and from the training data, respectively.

First, we study how the RMSE depends on the number of iteration for each EoR parameter. The result is plotted in Fig. 2. In this figure, we fix the number of neurons in the hidden layer at 14 and perform calculation with $10^6$ iterations. After initial fluctuations, the RMSE decreases and converges for >100 000 iterations. Next, we show
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Training error (RMSE, see the main text for definition) for $R_{mfp}$ (red), $\zeta$ (green) and $T_{vir}$ (blue) as functions of the number of iterations for the learning process, for a network with 14 neurons.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Training error (RMSE, see the main text for definition) for $R_{mfp}$ (solid line), $\zeta$ (dashed line) and $T_{vir}$ (dot-dashed line) as functions of the number of neurons for $10^5$ iterations in the learning process.</figcaption></figure>

how the RMSE depends on the number of neurons in the hidden layer after $10^5$ iterations. As we can see in Fig. 3, the RMSE is only weakly dependent on the number of neurons in the explored range. In the following section, we perform $10^5$ iterations for the back-propagation algorithm with 14 neurons in the hidden layer unless specifically mentioned.

## 4 RESULT

We apply the trained ANN to 54 test data sets. As for the training data sets, these test data sets consist of the 21 cm PS and the corresponding EoR parameters. We use the 21 cm PS data as input for the ANN and obtain estimated values for the EoR parameters. Then, we compare the estimated values with the true values that were used to compute the PS in the 21cmFast simulation, and evaluate how the ANN performs.

First, we apply the ANN to the 21 cm PS at a single redshift without including the effect thermal noise or sample variance. In Figs 4 and 5, we plot the EoR parameters estimated with the ANN from the 21 cm PS data against the true values. In Fig. 4, we show the result at $z = 12$. As we can see, $\zeta$ and $T_{vir}$ estimated with the ANN are in relatively good agreement with the true values. On the other hand, for $R_{mfp}$, the correct values are not recovered. This is because the 21 cm PS is not sensitive at $z = 12$ to $R_{mfp}$. On the other

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The EoR model parameter values computed by the ANN from the PS against the values used in the simulation at $z=12$. Note that the result for the Virial temperature is plotted in log scale. It is the same for the following figures.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The EoR model parameter values computed by the ANN from the PS against the values used in the simulation at $z=9$.</figcaption></figure>

hand, the 21 cm PS becomes sensitive to $R_{mfp}$ at lower redshift when reionization advances much more (since $R_{mfp}$ physically expresses the maximum bubble size, the effect of $R_{mfp}$ on the 21 cm PS is remarkable mostly at lower redshift). Thus the cost function is insensitive to changes in $R_{mfp}$ values; the learning process of the ANN is incomplete and systematic deviations remain. Note that the reason why the estimated value of $R_{mfp}$ seems to be constant at ~30 is that $R_{mfp}=30$ occurred more often than other values in the training data sets.

In Fig. 5, we show the same plots obtained from the PS at $z = 9$. The agreement between recovered and expected values extends over a larger range than at $z=12$, in particular for $R_{mfp}$. Indeed, the 21 cm PS at $z = 9$ is sensitive to $R_{mfp}$ and thus the ANN learning process works better.

In the previous two cases, we considered the 21 cm PS without any source of noise. We will now consider both the contribution of thermal noise and sample variance. We model the thermal noise as a Gaussian random field characterized by its PS. In an annulus in Fourier space with $N_c$ cells, the thermal noise PS is (e.g. Morales 2005; McQuinn et al. 2006; Yoshiura et al. 2015):

$$ \delta P_{\Delta T}^{21}(k, \theta) \approx \sqrt{\frac{1}{N_c} \frac{A_e x^2 y}{\lambda^2 B^2}} C^N(k, \theta), \quad (10) $$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The EoR model parameter values computed by the ANN from the PS against the values used in the simulation at z=12. In this case, we include both thermal noise and sample variance.</figcaption></figure>

where $\theta$ is the angle between **k** and line of sight, $A_e$ is the effective area of antenna, $\lambda$ is the observed wavelength and $x, y$ are factors converting uv space units into comoving wavenumber units and are determined by cosmology. $C^N$ is expressed by

$$C^N(k, \theta) = \left( \frac{\lambda^2 B T_{\text{sys}}}{A_e} \right)^2 \frac{1}{B t_k}, \qquad (11)$$

where $T_{\text{sys}}$ is the total system temperature, $B$ is the bandwidth and $t_k$ is the effective observing time. The spherically averaged noise can be obtained by summing over $\theta$ in a shell with radius $k$ shown in equation (12).

$$\delta P_{\Delta T}^{21}(k) = \left\{ \sum_{\theta} \left[ \frac{1}{\delta P_{\Delta T}^{21}(k, \theta)} \right]^2 \right\}^{-1/2} . \qquad (12)$$

We recommend the readers to read McQuinn et al. (2006). In our case, we assume SKA specification Dewdney. (2013). We estimate the sample variance from 10 simulations using different realizations of the initial conditions.

In the general case, we would include the noise at the level of the visibilities and compute the resulting noisy PS. Assuming that the noise is a random Gaussian field, we simply generate a noisy PS by using the following formula:

$$P_N^{ij}(k) = P^i(k) + N^j(k) \qquad (13)$$

where $P_N^{ij}(k)$ is the noisy PS, $P^i(k)$ is the PS produced by the simulation for a given set of parameters generically labelled by the index $j$ and $N^i(k)$ is a random draw from a Gaussian probability distribution with variance equal to the thermal noise PS or the sample variance. A different and independent draw is required for each value of $k$. For the learning sample we include only the sample variance. For each parameter set (70 possible values for $i$), we generate 50 realizations of the noise (labelled by $j = 1...50$). This means that we use 3500 different PS for the learning sample. This is necessary as we are not aware of a standard technique for directly including an uncertainty in the inputs of an ANN. For the test data sets, we add both thermal noise and sample variance to the 21 cm PS, using the same procedure.

In Fig. 6, we show the EoR parameters found by the ANN as functions of the values used in the simulations, using the 21 cm signal PS at $z=9$, including both thermal noise and sample variance. The difference between Figs 5 and 6 is not obvious at a glance.

<table><caption>**Table 1.** The mean chi-square values for the EoR model parameter values computed by the ANN from the PS as functions of the values used in the simulations. We show mean chi-square values for four types of data. The first column of this table is the chi-square values when we do not include any noise, and in second column we take thermal noise and sample variance into account. We use the 21 cm PS at $z=9$ for both cases. The third column is when we take redshift evolution ($z=9, 10, 11$) into account and the fourth column is when we take both thermal noise and sample variance into account and we reduce the number of training data samples.</caption><thead><tr><th></th><th>χ<sup>2</sup><sub>wo/noise</sub></th><th>χ<sup>2</sup><sub>w/noise</sub></th><th>χ<sup>2</sup><sub>w/noise, zevolution</sub></th><th>χ<sup>2</sup><sub>w/noise, reduced</sub></th></tr></thead><tbody><tr><td>R<sub>mfp</sub></td><td>0.228</td><td>0.258</td><td>0.172</td><td>0.262</td></tr><tr><td>ζ</td><td>0.271</td><td>0.288</td><td>0.168</td><td>0.290</td></tr><tr><td>log (T<sub>vir</sub>)</td><td>0.027</td><td>0.038</td><td>0.019</td><td>0.029</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparing the EoR model parameter values computed by the ANN against the values used in the simulation using the PS at $z=9, 10$ and $11$ as input data and including both thermal noise and sample variance.</figcaption></figure>

order to quantify the difference, we compute the mean chi-square value, $\chi^2$:

$$\chi^2 = \sqrt{\frac{1}{N_{\text{test}}} \sum_{i=1}^{N_{\text{test}}} \frac{(\theta_{\text{ANN}} - \theta_{\text{true}})^2}{\theta_{\text{true}}^2}} \qquad (14)$$

$\theta_{\text{ANN}}$ is the $R_{\text{mfp}}$, $T_{\text{vir}}$, $\zeta$ reconstructed by the ANN and $\theta_{\text{true}}$ is the value of the corresponding parameter used in the simulation. In Table 1, we compare the $\chi^2$ values for each of the parameters, with and without noise. As we can see, the $\chi^2$ values for the 21 cm PS with noises are worse than those without noises. Noise alters the efficiency of the learning process for the ANN.

Until now we have used as input data for the ANN the 21 cm PS at a single redshift. We now expand our input data to take the redshift evolution of the 21 cm PS into account. We use the 21 cm PS, which includes thermal noise and sample variance, at three redshifts $z=9, 10$ and $11$ as input data for the training and test data sets. Thus each set as $14 \times 3$ input values and the number of neurons in the input layer is also $14 \times 3=42$. We also increase the number of neurons from 14 to 42 in the hidden layer. Note that the number of training data sets is unchanged. In Fig. 7, we show parameter values computed by the ANN as functions of the values used in the simulations when the ANN is using the 21 cm PS at three different redshifts. Compared with the single-redshift case, it seems that the scatter is smaller. We quantify this in Table 1, and check that the scatter is indeed smaller for all three parameters. The reason why using the 21 cm PS at multiple redshifts improves the accuracy of
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The EoR model parameter values computed by the ANN for the PS against the values used in the simulation at $z=9$, using only 20 training data sets during the learning process.</figcaption></figure>

the ANN results is simply that we increase the information the ANN has to work with. At each redshift, the dependence and sensitivity of the 21 cm PS on the parameters are different. Therefore, the trained network can learn much more if we consider the redshift evolution of the 21 cm PS. We also notice that the chi-square values obtained by using multiple redshifts and including noise are smaller than those for a single redshift without any noise. We next check the influence of the size of the training sample. If we reduce the number of training data sets, we expect that the learning efficiency of the ANN will decrease. To check this, we train the ANN with a subsample of the training data sets ($N_{train}=20$). We show the results in Fig. 8. As expected, we can see that the scatter of reconstructed parameter values is larger when using fewer training data sets. This is quantified in Table 1. In particular, the scatter of $R_{mfp}$ increases more than for the other parameters. This is connected to the fact already established that the PS is not very sensitive to $R_{mfp}$. Thus, the size of the learning sample is crucial to make accurate EoR parameter predictions from observed data.

## 5 SUMMARY AND DISCUSSION

In this paper, we applied ANN to estimate EoR parameters from the 21 cm PS. We used 21cmFAST to produce 21 cm PS for different values of the following parameters: $R_{mfp}$, $\zeta$ and $T_{vir}$. We ran 70 simulations, which provided us with 70 training data sets made of the 21 cm PS and the corresponding EoR parameter values. With these data sets, we trained the ANN, which consists of an input layer, one hidden layer and an output layer. We achieved the training by applying the back-propagation algorithm.

In order to estimate the quality of the learning process of the ANN, we first show that the RMSE for the training data sets depends on the number of iterations of the back-propagation algorithm and on the number of neurons in the hidden layer. We found that the back-propagation algorithm converges as the number of iteration increases and weakly depends on the number of neurons in the hidden layer.

We next applied the trained network to test data sets made of the 21 cm PS at a single redshift including or not thermal noise and sample variance. As expected, we found that the chi-square values for the ANN-recovered values of the parameters as functions of the true (simulation) values are larger when including thermal noise and sample variance. We also found that the $\zeta$ and $T_{vir}$ values

<table><caption>**Table 2.** We compare the mean chi-square value obtained by the ANN with the $1\sigma$ fractional error obtained by 21CMMC for HERA and SKA observations at $z=9$. For comparison, we adopt the noise used in 21CMMC for HERA.</caption><thead><tr><th></th><th>χ<sup>2</sup><sub>SKA</sub></th><th>χ<sup>2</sup><sub>HERA</sub></th><th>1σ<sub>SKA</sub></th><th>1σ<sub>HERA</sub></th></tr></thead><tbody><tr><td>R<sub>mfp</sub></td><td>0.258</td><td>0.278</td><td>0.178</td><td>0.184</td></tr><tr><td>ζ</td><td>0.288</td><td>0.354</td><td>0.167</td><td>0.220</td></tr><tr><td>log (T<sub>vir</sub>)</td><td>0.038</td><td>0.040</td><td>0.024</td><td>0.033</td></tr></tbody></table>

obtained by the ANN were a better match to the true values than for $R_{mfp}$ due to the different sensitivity of the 21 cm PS to the different parameters.

Then we increased the information in the inputs given to the ANN by including the PS at three different redshifts. We found that the accuracy of parameters computed by the ANN was improved compared with those obtained based on the PS at a single redshift. This is simply that the neural network takes advantage in its learning process of the fact that the input holds more information and that the 21 cm PS is sensitive to different parameters at different redshifts.

One of the remarkable points of our study is that we can recover the EoR parameters from the PS using the ANN with good accuracy although we use only 70 training data sets. This number is quite small compared with that required for MCMC (e.g. $10^5$ realizations of the PS are used in 21CMMC; Greig & Mesinger 2015). In Table 1, we show the mean chi-square value for each parameter calculated by the ANN. To enable a quantitative comparison with 21CMMC, we perform a parameter estimation with the ANN for HERA and SKA observations and compare to the $1\sigma$ uncertainties obtained with 21CMMC (for HERA, we actually use the noise data provided in 21CMMC package). We show the result in Table 2. Using a single redshift, the accuracy of our reconstruction is of the same order as that of 21CMMC, although somewhat worse. Remember however that it does improve when using several redshifts, as does the accuracy of 21CMMC. However, note that the MCMC method finds the local shape of the likelihood space, while the ANN method only finds the average error over a swath of parameter space. So, while the ANN method requires a far sparser exploration of the parameter space, 21CMMC does provide more detailed information. In addition, we need to mention the effect of the foregrounds. The foregrounds wedge may hinder the use of part of the PS. Thus, taking this into account would likely affect the constraints on the $R_{mfp}$ because it depends heavily on the available range of k-space. We defer an ANN analysis incorporating the effect of foregrounds to a future work. Both the training process for the ANN and the parameter estimation from the observed data using the trained ANN have a very small computing cost compared to that required for producing training data sets.

In this work, we used the 21 cm PS generated by 21cmFAST. However, it is worth applying the ANN to other EoR modelling methods including full simulations in order to check its ability to reconstruct the astrophysical parameters. In our model, we focused on the EoR and used the 21 cm PS. Thus, a possible development is to expand the model beyond the EoR and include the cosmic dawn. During the cosmic dawn, X-ray heating of the IGM plays an important role. Thus, it is important to incorporate X-ray heating models and study their impact on the 21 cm (Baek et al. 2010; Mesinger et al. 2013; Ewall-Wice et al. 2016; Fialkov et al. 2017). Recently, Cohen et al. (2016) explored the behaviour of the 21 cm global signal in a large parameter space from $z=40$ to $z=6$. We believe it is worth applying the ANN method using the 21 cm global
signal as an input and try to reconstruct parameter characteristic of the cosmic dawn and EoR both.

Another possible application of ANNs is to use the 21 cm tomography as an input. Then the learning sample is made of light cone data from 21 cm simulations instead of the 21 cm PS. Although the 21 cm PS is a powerful tool, the 21 cm signal is not Gaussian, and thus the PS does not contain all information. The 21 cm light cone data on the other hand holds all the information that will ever be accessible through observations. It also contains irrelevant information that the ANN will have to filter out: the random phases of the initial density field.

If we consider a large survey field of view and depth, that is a large box size for the learning sample simulation, we will need to treat a large amount of data with the ANN. Machine learning methods have been developed to treat such big data. One of the methods that have received attention is deep learning (e.g. Schmidhuber 2015; LeCun, Bengio & Hinton 2015). In our case, we used only one hidden layer in the ANN. A deep learning neural network consists of one input layer, multiple hidden layers and one output layer. Deep learning has been applied successfully to image recognition. Therefore, one of the our future works will be to apply deep learning to 21 cm light cones from full simulations data.

## APPENDIX

Here, we show the back propagation algorithm in the case of multiple hidden layers. By analogy with equation (8), we can express the derivative of $E_n$ with respect to the weight in the lth layer, $w_{ij}^{(l)}$ ($l=1, 2, ...L-1$) as

$$ \frac{\partial E_n}{\partial w_{ij}^{(l)}} = \frac{\partial E_n}{\partial s_i^{(l)}} \frac{\partial s_i^{(l)}}{\partial w_{ij}^{(l)}}, \quad (15) $$

where $s_i^{(l)} = \sum_j w_{ij}^{(l)} t_j^{(l-1)}$. Since the changes in $s_i^{(l)}$ are transmitted to $E_n$ through each $s_k^{(l+1)}$ in the (l+1)th layer, the derivative of $E_n$ with respect to $s_i^{(l)}$ can be expressed as

$$ \frac{\partial E_n}{\partial s_i^{(l)}} = \sum_k \frac{\partial E_n}{\partial s_k^{(l+1)}} \frac{\partial s_k^{(l+1)}}{\partial s_i^{(l)}}, \quad (16) $$

Remember that $s_k^{(l+1)}$ can be expressed with the activation function $\phi$ as $s_k^{(l+1)} = \sum_j w_{kj}^{(l+1)} t_j^{(l)} = \sum_k w_{kj}^{(l+1)} \phi(s_j^{(l)})$,

$$ \frac{\partial s_k^{(l+1)}}{\partial s_i^{(l)}} = w_{ki}^{(l+1)} \phi'(s_i^{(l)}). \quad (17) $$

Here, we define $\delta_i^{(l)} \equiv \partial E_n / \partial s_i^{(l)}$, and then we can re-write equation (16) as

$$ \delta_i^{(l)} = \sum_k \delta_k^{(l+1)} w_{ki}^{(l+1)} \phi'(s_i^{(l)}). \quad (18) $$

Combining equation (18) with $\partial s_i^{(l)} / \partial w_{ij}^{(l)} = t_j^{(l-1)}$, equation (15) can be re-written simply as

$$ \frac{\partial E_n}{\partial w_{ij}^{(l)}} = \delta_i^{(l)} t_j^{(l-1)}. \quad (19) $$

This form tells us that we can easily obtain the derivative of the cost function with respect to $w_{ij}^{(l)}$, which connects the neuron j in the $(l-1)$th layer to the neuron i in the lth layer, as the product of $\delta_i^{(l)}$ and $t_j^{(l-1)}$. As shown in equation (18), we start to calculate $\delta_i^{(l)}$ from the output layer ($l=L$) to the input layer. This is why this algorithm is called 'back propagation'. If we use equation (5) as the cost function, we easily calculate $\delta_i^{(L)} = y_i - d_i$.

## ACKNOWLEDGEMENTS

This work is benefited from a grant from the French ANR funded project ORAGE (ANR-14-CE33-0016). We thank G.Mellema, A. Fialkov, S. Majumdar, S.Giri and K. Hasegawa for their useful comments and thank S. Yoshiura for providing the thermal noise data.

## REFERENCES

Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93

Agarwal S., Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2012, MNRAS, 424, 1409

Agarwal S., Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2014, MNRAS, 439, 2102

Ali Z. S. et al., 2015, ApJ, 809, 61

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

Ball N. M., Brunner R. J., 2010, Int. J. Mod. Phys. D, 19, 1049

Banerji M. et al., 2010, MNRAS, 406, 342

Barkana R., Loeb A., 2008, MNRAS, 384, 1069

Beardsley A. P. et al., 2016, ApJ, 833, 102

Bernardi G. et al., 2016, MNRAS, 461, 2847

Bloom J. S., Richards J. W., 2012, in Way M. J., Scargle J. D., Ali K. M., Srivastava A. N., eds, Advances in Machine Learning and Data Mining for Astronomy. CRC Press, Boca Raton, p. 89

Bromm V., 2013, Rep. Prog. Phys., 76, 112901

Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23

Ciardi B., Ferrara A., Governato F., Jenkins A., 2000, MNRAS, 314, 611

Ciardi B. et al., 2006, MNRAS, 366, 689

Cohen A., Fialkov A., Barkana R., Lotem M., 2016, MNRAS, preprint (arXiv:1609.02312)

Collister A. A., Lahav O., 2004, PASP, 116, 345

Cybenko G., 1989, Math. Control Signals Syst., 2, 303

DeBoer D. R. et al., 2017, PASP, 129, d5001

Dewdney P., Turner W., Millenaar R., McCool R., Lazio J., Cornwell T., 2013, SKA1 System Baseline Design, Document number SKA-TEL-SKO-DD-001 Revision 1. Available at: https://www.skatelescope.org/wp-content/uploads/2012/07/SKA-TEL-SKO-DD-001-1_BaselineDesign1.pdf

Ewall-Wice A. et al., 2016, MNRAS, 460, 4320

Fan H. X., Carilli L. C., Keating B. G., 2006, ARA&A, 44, 415

Fialkov A., Barkana R., Visbal E., Tseliakhovich D., Hirata C. M., 2013, MNRAS, 432, 2909

Fialkov A., Cohen A., Barkana R., Silk J., 2017, MNRAS, 464, 3498

Folkes S. R., Lahav O., Maddox S. J., 1996, MNRAS, 283, 651

Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1

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

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

Greig B., Mesinger A., 2015, MNRAS, 449, 4246

Greig B., Mesinger A., Pober J. C., 2016, MNRAS, 455, 4295

Harker G. J. A. et al., 2009, MNRAS, 393, 1449

Hasegawa K., Semelin B., 2013, MNRAS, 428, 154

Hornik K., Stinchcombe M., White H., 1989, Neural Netw., 2, 359

Iliev I. T. et al., 2006, MNRAS, 369, 1625

Iliev I. T. et al., 2012, MNRAS, 423, 2222

Iliev I. T. et al., 2014, MNRAS, 439, 725

Jacobs D. C. et al., 2015, ApJ, 801, 51

Jensen H. et al., 2016, ApJ, 827, 5

Kamdar H. M., Turk M. J., Brunner R. J., 2016a, MNRAS, 455, 642

Kamdar H. M., Turk M. J., Brunner R. J., 2016b, MNRAS, 457, 1162

Kim K. et al., 2015, Class. Quantum Gravity, 32, 245002

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

Kubota K., Yoshiura S., Shimabukuro H., Takahashi K., 2016, PASJ, 68, 61

Lahav O., Naim A., Sodré L., Jr, Storrie-Lombardi M. C., 1996, MNRAS, 283, 207

LeCun Y., Bengio Y., Hinton G., 2015, Nature, 521, 436
Liu A., Parsons A. R., 2016, MNRAS, 457, 1864

McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815

McQuinn M. et al., 2007, MNRAS, 377, 1043

McQuinn M., Oh S. P., Faucher-Giguère C.-A., 2011, ApJ, 743, 82

Mellema G. et al., 2013, Exp. Astron., 36, 235

Mesinger A., Furlanetto S., 2007, ApJ, 669, 663

Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955

Mesinger A., Ferrara A., Spiegel D. S., 2013, MNRAS, 431, 621

Morales M. F., 2005, ApJ, 619, 678

Morii M. et al., 2016, PASJ, 68, 104

Mukund N., Abraham S., Kandhasamy S., Mitra S., Sajeeth Philip N., 2016, preprint (arXiv:1609.07259)

Parsons A. R. et al., 2014, ApJ, 788, 106

Patil A. H. et al., 2014, MNRAS, 443, 1113

Planck Collaboration XIII, 2016a, A&A, 594, A13

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

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

Pober J. C. et al., 2015, ApJ, 809, 62

Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys., 75, 086901

Pritchard J. et al., 2015, Proc. Sci., Cosmology from EoR/Cosmic Dawn with the SKA. SISSA, Trieste, PoS#012

Ricotti M., Gnedin N. Y., Shull J. M., 2002, ApJ, 575, 49

Rottgering H., 2003, New Astron. Rev., 47, 405

Rumelhart D. E., Hinton G. E., Williams R. J., 1986, Nature, 323, 533

Samui S., Samui Pal S., 2017, New Astron., 51, 169

Santos M. G., Ferramacho L., Silva M. B., Amblard A., Cooray A., 2010, MNRAS, 406, 2421

Schmidhuber J., 2015, Neural Networks, 61, 85

Semelin B., Iliev I., 2015, Proc. Sci., The physics of Reionization: processes relevant for SKA observations. SISSA, Trieste, PoS#013

Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., Ichiki K., 2015, MNRAS, 451, 4986

Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., Ichiki K., 2016, MNRAS, 458, 3003

Shimabukuro H., Yoshiura S., Takahashi K., Yokoyama S., Ichiki K., 2017, MNRAS, 468, 1542

Sobacchi E., Mesinger A., 2013, MNRAS, 432, 3340

Sobacchi E., Mesinger A., 2014, MNRAS, 440, 1662

Tingay S. J. et al., 2013, PASA, 30, 7

Trac H., Cen R., 2007, ApJ, 671, 1

Vanzella E. et al., 2004, A&A, 423, 761

Wise J. H., Cen R., 2009, ApJ, 693, 984

Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 652, 6

Yoshiura S. et al., 2015, MNRAS, 451, 4785

Zahn O. et al., 2007, ApJ, 654, 12

Zahn O. et al., 2011, MNRAS, 414, 727

This paper has been typeset from a TeX/LATeX file prepared by the author.