Constraining the astrophysics and cosmology from 21cm tomography using deep learning with the SKA

TL;DR

This paper demonstrates that deep learning, specifically CNNs, can accurately extract key astrophysical and cosmological parameters from simulated 21cm maps, even with realistic SKA noise, aiding future cosmological studies.

Contribution

The authors develop CNN architectures capable of simultaneously estimating multiple astrophysical and cosmological parameters from 21cm maps, including realistic SKA-like noise conditions.

Findings

CNNs recover parameters with over 92% accuracy under noise

Accuracy improves to over 99% at low redshift and neutral fraction

Future 21cm observations can tightly constrain parameters with few frequency channels

Abstract

Future Square Kilometre Array (SKA) surveys are expected to generate huge datasets of 21cm maps on cosmological scales from the Epoch of Reionization (EoR). We assess the viability of exploiting machine learning techniques, namely, convolutional neural networks (CNN), to simultaneously estimate the astrophysical and cosmological parameters from 21cm maps from semi-numerical simulations. We further convert the simulated 21cm maps into SKA-like mock maps using the detailed SKA antennae distribution, thermal noise and a recipe for foreground cleaning. We successfully design two CNN architectures (VGGNet-like and ResNet-like) that are both efficiently able to extract simultaneously three astrophysical parameters, namely the photon escape fraction (f$_{\rm esc}$), the ionizing emissivity power dependence on halo mass ($C_{\rm ion}$) and the ionizing emissivity redshift evolution index ($D_{\rm ion}$), and three cosmological parameters, namely the matter density parameter ($\Omega_{m}$), the dimensionless Hubble constant ($h$), and the matter fluctuation amplitude ($\sigma_{8}$), from 21cm maps at several redshifts. With the presence of noise from SKA, our designed CNNs are still able to recover these astrophysical and cosmological parameters with great accuracy ($R^{2} > 92\%$), improving to $R^{2} > 99\%$ towards low redshift and low neutral fraction values. Our results show that future 21cm observations can play a key role to break degeneracy between models and tightly constrain the astrophysical and cosmological parameters, using only few frequency channels.

Tables

Table 1. Table 1: Summary of our assumed SKA array design.

Array design	866 compact core
Station diameter, $D$ [m]	35
Station area, $A\left[m^2\right]$	962.11 ${\left(\frac{110}{\nu \left[\mathrm{MHz}\right]}\right)}^2$
System temperature [K] ($T_{\mathrm{sys}}=T_{\mathrm{sky}}+T_{\mathrm{rcvr}}$)	1.1 $T_{\mathrm{sky}}$ + 40
Total observation time tint [h]	1000
Frequency resolution $\Delta \nu$ [kHz]	48
Redshift	10, 9, 8 ,7
Frequency [MHz]	129 , 142, 158, 178
FWHM [arcmin]	1.37, 1.24, 1.12, 0.99
Beam angle $\theta$ [rad]	0.066, 0.06, 0.054, 0.048
Default wedge slope $m$, Equation. (4)	0.27, 0.23, 0.19, 0.15

Table 2. Table 2: The architecture of Network I for this study.

	Layer	Output shape
1	Input	(1, 200, 200)
2	3$\times$3 Convolutional Layer	(32, 200, 200)
3	3$\times$3 Convolutional Layer	(32, 200, 200)
4	Batch Normalization	–
5	ReLU Activation	–
6	2$\times$2 Max Pooling	(32, 100, 100)
7	3$\times$3 Convolutional Layer	(64, 100, 100)
8	3$\times$3 Convolutional Layer	(64, 100, 100)
9	Batch Normalization	–
10	ReLU Activation	–
11	2$\times$2 Max Pooling	(64, 50, 50)
12	3$\times$3 Convolutional Layer	(128, 50, 50)
13	3$\times$3 Convolutional Layer	(128, 50, 50)
14	Batch Normalization	–
15	ReLU Activation	–
16	2$\times$2 Max Pooling	(128, 25, 25)
17	3$\times$3 Convolutional Layer	(256, 25, 25)
18	3$\times$3 Convolutional Layer	(256, 25, 25)
19	Batch Normalization	–
20	ReLU Activation	–
21	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(1024)
22	Batch Normalization	–
23	ReLU Activation	–
24	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(1024)
25	Batch Normalization	–
26	ReLU Activation	–
27	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(1024)
28	Batch Normalization	–
29	ReLU Activation	–
30	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(6)

Table 3. Table 3: The architecture of Network II for this study.

	Layer	Output shape
1	Input	(1, 200, 200)
2	Convolutional Layer	(16, 200, 200)
3	Batch Normalization	–
4	ReLU Activation	–
5	Residual Layer (3 Residual Blocks)	(16, 100, 100)
6	Residual Layer (6 Residual Blocks)	(32, 50, 50)
7	Residual Layer (6 Residual Blocks)	(64, 25, 25)
8	Residual Layer (3 Residual Blocks)	(128, 13, 13)
9	Inception Module	(240, 13, 13)
10	Max Pooling	(240, 7, 7)
11	Inception Module	(240, 7, 7)
12	Inception Module	(256, 7, 7)
13	Inception Module	(288, 7, 7)
14	Max Pooling	(288, 4, 4)
15	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(1024)
16	Batch Normalization	–
17	ReLU Activation	–
18	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(1024)
19	Batch Normalization	–
20	ReLU Activation	–
21	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(1024)
22	Batch Normalization	–
23	ReLU Activation	–
24	$\mathrm{Fully}\mathrm{Connected}\mathrm{Layer}$	(6)

Table 4. Table 4: The hyperparameters and optimisers used to train the algorithms.

	optimiser	learning rate	batch	cost function
I	Nesterov	0.005	128	${\ell }_1$ norm
II	Adam	0.01	128	rmse

Equations

$$\text{centering}{f}_{\mathrm{esc}}{R}_{\mathrm{ion}}\ge R_{\mathrm{rec}},\text{@add@centering}$$

$$\frac{R_{\mathrm{rec}}}{V}=9.85\times 10^{-24}(1+z{)}^{5.1}{\left[\frac{{\left(\Delta /1.76\right)}^{0.82}}{1+{\left(\Delta /1.76\right)}^{0.82}}\right]}^4,$$

$$\frac{R_{\mathrm{ion}}}{M_h}=1.1\times 10^{40}(1+z{)}^{D_{\mathrm{ion}}}{\left(\frac{M_h}{9.51\times 10^7}\right)}^{C_{\mathrm{ion}}}\exp {\left(\frac{-9.51\times 10^7}{M_h}\right)}^{3.0},$$

$$m=\frac{D{H}_{0}E(z)\sin \theta }{c(1+z)},$$

$$\sqrt{⟨|N{|}^2⟩}[\mathrm{Jy}]=\frac{2k_{\mathrm{B}}{T}_{\mathrm{sys}}}{A\sqrt{\Delta \nu t_{\mathrm{int}}}},$$

$$R^2=1-\frac{{\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2}{{\sum }_{i=1}^n(y_i-\bar{y}{)}^2},$$

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}\sum (y_{\mathrm{predicted}}-y_{\mathrm{true}}{)}^2},$$

Full text

††pubyear: 2015††pagerange: Constraining the astrophysics and cosmology from 21cm tomography using deep learning with the SKA–LABEL:lastpage

1 Introduction

The last global phase transition in the Universe, known as the Epoch of Reionization (EoR), marks the time at which the first stars gradually reionized the Inter-Galactic Medium (IGM) and the Universe transitioned from highly neutral-opaque to a highly ionized-transparent state (for a review, see e.g. Loeb & Barkana, 2001). This epoch represents a crucial period in the Universe’s history, particularly with regard to the formation and evolution of early galaxies.

Constraining the astrophysical and cosmological parameters has been the focus for most observational and theoretical studies. Several techniques have been developed to constrain the cosmological parameters (e.g matter density parameter $\Omega_{m}$ and Hubble constant $H_{0}$) such as using the Cosmic Microwave Background (CMB) anisotropies measurements (e.g. Hinshaw et al., 2013; Planck Collaboration et al., 2016), Sunyaev-Zel’dovich cluster surveys (e.g. Battye & Weller, 2003), galaxy clusters in optical and X-ray bands (e.g. Moscardini et al., 2001), gamma ray burst X-ray afterglow light curves (e.g. Cardone et al., 2010), lensed GW+EM signals (e.g. Li et al., 2019), Ly-$\alpha$ forest power spectrum and COBE-DMR (e.g. Phillips et al., 2001), large-scale clustering of SDSS luminous red galaxies (e.g. Padmanabhan et al., 2007), and a joint CMB and weak lensing analysis (e.g. Contaldi et al., 2003). On the other hand, several works have attempted to constrain the astrophysical parameters (e.g. the photon escape fraction, $f_{\rm esc}$, and ionizing emissivity evolution, $\dot{N}_{\rm ion}$), using Ly-$\alpha$ forest measurements (e.g. Becker & Bolton, 2013), Lyman continuum (LyC) radiation from local galaxies (e.g. Leitet et al., 2013), and inferred constraints by tuning different theoretical models to other measurements (e.g. Mitra et al., 2015; Finlator et al., 2015a).

While all these methods show different levels of success to place constraints on various parameters, tighter constraints are expected to come from the EoR through measurements of the 21cm fluctuations on cosmological scales. With its strong dependence on the ionization and density fields, the 21cm signal carries a wealth of information that is important in order to understand early stages of galaxy formation and evolution. In this light, many radio interferometer experiments, such as the Low Frequency Array (LOFAR; van Haarlem et al., 2013), the Precision Array for Probing the Epoch of Reionization (PAPER; Parsons et al., 2010), the Murchison Wide field Array (MWA; Bowman et al., 2013), the Giant Metrewave Radio Telescope (GMRT; Paciga et al., 2011), the Hydrogen Epoch of Reionization Array (HERA; DeBoer et al., 2017) and Square Kilometer Array (SKA; Mellema et al., 2013) are devoted to detecting reionization in the near future. These growing observational efforts require equivalent efforts in both the theoretical and statistical sides, in order to prepare for extracting all possible information and constrain the cosmological and astrophysical parameters from future 21cm surveys.

Several studies have already shown that combining the 21cm power spectrum with Markov Chain Monte Carlo (MCMC) analysis is a powerful technique to obtain tighter constraints and break degeneracy between models (e.g. Greig & Mesinger, 2015; Liu et al., 2016; Pober et al., 2016; Hassan et al., 2017; Park et al., 2019). Besides the power spectrum, future 21cm surveys, the SKA in particular, are also expected to generate huge imaging datasets for the 21cm fluctuations on large scales that will contain more information than the power spectrum. Going beyond the power spectrum has been the target of many studies  (e.g. Bharadwaj & Pandey, 2005; Barkana & Loeb, 2008; Watkinson & Pritchard, 2015; Majumdar et al., 2018), in which more information can be obtained through investigating the non-gaussian nature of the 21cm signal using higher-order statistics such as the bispectrum. To efficiently use the 21cm information stored in the 2-dimensional 21cm maps, Convolutional Neural Networks (CNNs) have been a very successful deep learning tool to recover the astrophysical parameters during reionization (Gillet et al., 2019), to learn the reionization history (La Plante & Ntampaka, 2019; Mangena et al., 2019), to emulate reionization simulations (Chardin et al., 2019), and to identify reionization sources from different models (Hassan et al., 2019). However, the astrophysical parameter recovery by Gillet et al. (2019) ignores the instrumental effects as an initial proof-of-concept study. Accounting for these effects such as the angular resolution, foreground cleaning, and thermal noise, are all crucial in order to add realism to the simulated 21cm images as we prepare for the 21cm era.

In this work, we take a step further to design two different CNNs to simultaneously estimate several parameters from 21cm maps at several redshifts and different stages through reionization. We here that assume all observations at different redshifts are performed independently. We simply take maps from different redshifts and apply the instrumental noise directly on each map assuming a single frequency channel of a size $\sim$ 50 KHz (i.e. simulation resolution). We finally combine the maps from different redshifts to create our training datasets. We note that learning from light-cones is beyond the scope of the current work. Our aim is to provide a network that is able to predict parameters without requiring the redshift nor neutral fraction as inputs, which is a more flexible design. Three astrophysical parameters are evaluated: the photon escape fraction (fesc), the ionizing emissivity power dependence on halo mass ($C_{\rm ion}$) and the redshift evolution index ($D_{\rm ion}$). Additionally, we estimate three cosmological parameters: the matter density parameter ($\Omega_{m}$), the dimensionless Hubble constant ($h$), and the matter fluctuation amplitude ($\sigma_{8}$). To assess the ability of future 21cm tomography to constrain these parameters, we follow the recipe presented in Hassan et al. (2019) to add a physically motivated and realistic 21cm noise to large scale 21cm maps that are produced using our semi-numerical model, SimFast21 (Santos et al., 2008, 2010). This paper is organized as follows: we first describe our suite of simulations of the 21cm signal and noise in §2. We then present the two network designs in §3 and the training dataset in §4. We present the main results in §5, and draw our concluding remarks in §6.

2 simulations

2.1 Semi-Numerical Model, SimFast21

We use the Instantaneous model of our semi-numerical simulations SimFast21, that has been developed in Hassan et al. (2016), to improve over previous implementations of the ionizing source and sink populations in these semi-numerical simulations. In addition, it has been recently shown that this model is in a relatively good agreement with predictions from our radiative transfer simulation (ARTIST; Molaro et al., 2019), particularly in terms of the morphology and power spectrum of the ionization and 21cm fields. However, the reionization history can be quite different for the same photon escape fraction value. This arises from violation of photon conservation which is an intrinsic problem in the use of excursion set-formalism in semi-numerical simulations (Zahn et al., 2007; Paranjape et al., 2016; Hassan et al., 2017). As indicated by ARTIST, as a temporary solution all our photon escape fraction predictions can be adjusted by a factor of 20% to account for the photon conservation problem. We here briefly describe the simulation ingredients, and defer to Santos et al. (2010) for the full details of the simulation algorithm, and to Hassan et al. (2016) for the Instantaneous model development.

The dark matter density is generated in the linear regime from a Gaussian distribution using a Monte-Carlo approach. Evolving the density field to non-linear regime is performed through the Zel’dovich (1970) approximation. Halos are then generated using the excursion-set formalism (ESF). Ionized regions are identified using a similar form of the ESF that is based on a direct comparison between the instantaneous rates of ionization $R_{\rm ion}$ and recombination $R_{\rm rec}$ in spherical regions of decreasing sizes as specified by the ESF. Regions are flagged as ionized if:

[TABLE]

where $f_{\rm esc}$ is the escape fraction. The $R_{\rm rec}$ is obtained from a radiative transfer simulation (Finlator et al., 2015b), in order to account for the clumping effects below our cell size. The $R_{\rm rec}$ is parameterized as a function of overdensity $\Delta$ and redshfit $z$ as follows:

[TABLE]

where V refers to the cell volume. The $R_{\rm ion}$ parameterization is derived from a combination of the radiative transfer simulation (Finlator et al., 2015b), and a larger hydrodynamic simulation (Davé et al., 2013) that both have been shown to reproduce wide range of observations, including low-$z$ observations. The $R_{\rm ion}$ is parameterized as a function of halo mass $M_{\rm h}$ and redshift $z$ as follows:

[TABLE]

where the best fit values of the ionizing emissivity dependence on halo mass $C_{\rm ion}$ and redshift $D_{\rm ion}$ were found to be $C_{\rm ion}=0.41$ and $D_{\rm ion}=2.28$, respectively. Later, we will change these parameters to generate the training and testing datasets. Note that Equation (3) shows that $R_{\rm ion}$ scales as M${}_{\rm h}^{1.41}$, which is consistent with the SFR$-$Mh relation previously found by Finlator et al. (2011). We defer to Hassan et al. (2016) for the full details on the derivation of the $R_{\rm ion}$ and $R_{\rm rec}$ fitting functions and their effects on several reionization key observables.

2.2 21cm Instrument Simulation

We here describe the method used to account for various instrumental effects following the recipe developed in Hassan et al. (2019). We briefly review this method below and refer the interested readers to Hassan et al. (2019) for detailed information and complete steps of how we convert an 21cm simulated map into a mock map according to the assumed array design. In this work, we restrict our analysis to the SKA proposed design and leave a more detailed comparison between different arrays, such as HERA and LOFAR, for future works. The instrumental noise is applied separately on each redshift assuming a single frequency channel corresponding to the map size ($\sim$ 50 kHz). We leave to future works learning from the light-cones by considering many frequency channels over a large bandwidth in the analysis.

The 21cm Instrument simulation pipeline consists of three parts:

• •

Foreground cleaning: Foreground contaminated modes of the signal lie inside the foreground wedge in the $k_{\perp}-k_{\parallel}$ plane. The foreground wedge slope (m) is given by:

[TABLE]

where $H_{0}$ is the Hubble parameter, $c$ is the speed of light, $E(z)\equiv\sqrt{\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}}$, and $\theta$ is the beam angle. To clean foregrounds, we simply zero out all modes within the wedge, satisfying k${}_{\parallel}<$ m k*⟂*. For the same experiment, the slope increases with redshift, which means more modes are removed at higher redshifts. We quote exact wedge slope values for the SKA at our redshifts of interest in Table 1. This the first step of the noise pipeline to clean foregrounds from the 3-dimensional co-eval cubes.

• •

Angular resolution: we account for the angular resolution of a given array by exploiting its detailed baseline distribution, via the $uv$-coverage, which is a measure of the baseline intensity observing the signal modes in directions perpendicular to the sightline. The $uv$-coverage is computed using the 21cmSense package111https://github.com/jpober/21cmSense from our assumed SKA antennae distribution. We then Fourier transform the simulated 21cm map and set the signal to zero at k*⟂* modes whose $uv$-coverage is zero222Modes with zero $uv$-coverage lie outside the angular resolution of the experiment.. We additionally smooth down the simulated maps using a Gaussian filter whose full width half maximum (FWHM) is given by: FWHM = $\lambda_{21cm}$ (1+z)/B, whereas the maximum baseline length B=5,834 m for our assumed SKA design, and $\lambda_{21cm}$ is the rest frame wavelength of the 21cm signal. This sets the minimum angular resolution for our assumed SKA design. For instance, our simulated maps initially have an angular resolution of $\sim$ 0.3′ at z=7, that are smoothed to have a lower angular resolution of $\sim$ 1′ according to the FWHM at this redshift. Exact angular resolution values as a function of redshift are quoted in Table 1. The angular resolution recipe is applied on maps extracted from the 3-dimensional foreground filtered boxes from the previous step.

• •

Thermal noise: the thermal noise is uncorrelated between measurements, and can be drawn from a Gaussian distribution of unit mean and standard deviation (Zaldarriaga et al., 2004) given by:

[TABLE]

where tint here is the integration time to observe a single visibility at a frequency resolution $\Delta\nu$, and $k_{\rm B}$ is the Boltzmann constant. The total system temperature Tsys and other parameters are summarized in Table 1. Having generated the thermal noise in 2D grid using the above equation in Fourier space, we further suppress the noise by the amount of the $uv$-coverage $N_{uv}$ by a factor of $\sim 1/\sqrt{N_{uv}}$. We finally inverse Fourier transform the noise map and add it to the angular resolution - foreground filtered signal map to form our mock 21cm map.

Using this pipeline with parameters listed in Table 1, the rms brightness temperature (noise level) is about $\sim$ 3 mK at z=8, consistent with previous estimates (e.g. see Furlanetto et al., 2006; Kakiichi et al., 2017; Giri et al., 2018). This pipeline is used to add realism to our simulated training and testing dataset, in order to assess the ability of future SKA 21cm surveys to constrain the astrophysical and cosmological constraints. In Figure 1, we show an example of four randomly selected 21cm maps (top) with their mock versions (bottom) from our training datatset for different set of astrophysical and cosmological parameters as quoted in the subtitles. These maps are generated from different simulations realizations of a box size of L=150 Mpc, number of cells N=200, resulting in a resolution of 0.75 Mpc. We find that most of the large and small scale ionized bubbles are still present after adding the instrumental effects. This is due to the high angular resolution of our assumed SKA design as well as the high $uv$-coverage that extends down to a very small scales ($\sim$ 3.5 h/Mpc) during these epochs.

However, fully ionized ($x_{\rm HI}<0.01$) and fully neutral maps ($x_{\rm HI}>0.99$), as described later in §4, are already excluded from the training sample, since they are identical for different set of parameters. Distinguishing identical maps is challenging for neural networks, where more information, such as redshift evolution, is required to assist parameter recovery at these extreme limits. When the Universe is highly ionized (e.g. $x_{\rm HI}\sim 0.1-0.2$, see column 4 in Figure 1), the noise dominates but nevertheless the residual neutral patches can still be seen and recognized. These residual patches are usually different for different set of parameters, which might help the network to distinguish between maps and parameters. On the other hand, in the beginning of reionization, the ionized regions are very small due to the small number of sources and ionizing photons. The noise then contaminates and fills these small ionized regions (e.g. $x_{\rm HI}\sim 0.8-0.9$, see column 1 in Figure 1), and hence maps might look similar to those from a fully neutral Universe. This makes recognizing the prominent signal features more challenging, and many of the reionization realizations for a highly neutral Universe become approximately indistinguishable. This might impact the parameter recovery from a highly neutral IGM, which basically exists at high redshifts where the noise is stronger.

3 Network architecture

We consider two types of network in this work. It is worth reiterating that our main objective is to able to infer both astrophysical $\{f_{\rm esc},C_{\rm ion},D_{\rm ion}\}$ and cosmological $\{\Omega_{\rm m},\textit{h},\sigma_{8}\}$ parameters simultaneously from their corresponding 21cm maps. To this end, our main focus is to simply explore different network designs with different layout in width and depth as an attempt to achieve our goal.

The first architecture (network I) considered for our investigation is given in Table. 2. It is slightly similar to VGG-Net (Simonyan & Zisserman, 2014) in terms of chaining convolutional layers before downsampling, however the key difference here is that each stage333which we refer to as mapping the input x without reducing the dimensions ($\rm{weight}\times\rm{height}$)., we have two convolutional layers with same number of feature maps in a row followed by batch normalisation and ReLU activation (Conv+Conv+BN+ReLU) as shown in Figure 2. We note that the representation N x N + M(S) denotes the kernel size (N x N) and the stride (M) of a convolutional layer. In total, we have four stages, each with a Conv+Conv+BN+ReLU layer followed by a Max Pooling with stride = 2 to reduce the dimensions of the inputs444In other words, the ouputs from the previous stage., and four dense layers each with 1024 units with the exception of the output layer, which has only 6 units corresponding to the number of inferred parameters. This network design is also similar to the previous design used in our reionization models classifier (Hassan et al., 2019), except that the convolutional layers used here are wider and no dropout seems to be needed. This is consistent with the disharmony observed between batch normalization and dropout (Li et al., 2018). Similar to our previous works in the classifier, we initialize the network weights using a generalized form of Xavier initializer (Glorot & Bengio, 2010) that is also called the Variance Scaling initializer, in which the random numbers are drawn from a zero mean Gaussian distribution whose variance is equal to the inverse of the average of the number of input and output neurons. This initializer ensures that the variance of the input data is preserved as it propagates through the network layers.

Our second architecture, which we simply name network II, is based on a combination of residual layers (He et al., 2016) and inception modules (Szegedy et al., 2015) as shown in Table. 3. The inputs, as described in §2, are first fed into a convolutional layer followed by a batch normalization (Ioffe & Szegedy, 2015) before a ReLU activation (Conv+BN+ReLU). This is then followed by four residual layers, each composed of three, six, six and three residual blocks respectively. It was shown in He et al. (2016) that the resulting error (both training and testing) of deeper architecture tends to be larger than that of shallower architecture. Therefore they proposed a residual layer which allowed them to increase the depth of the model in order to gain better performance. In contrast with network I, instead of using simple convolutional layers we stack residual blocks, which are achieved with the schematic shown in Figure. 3 (right panel) where the residual learning is constructed using a Conv+BN+ReLU+Conv+BN+ReLU+Conv+BN layer. Depending on whether there is downsampling (Fig. 3 left panel) through the chain of convolutional layers in a residual block, the input needs to be downsampled using a Conv+BN layer to match the dimension of the output of the chain of convolutional layers.

In each residual layer, the downsampling occurs at the first residual block. There are variants of deep residual networks, but in essence what we consider here is such that the network performance is optimized for our specific task.

As proposed by Szegedy et al. (2015), in order to improve the recognition of more complex features at the higher levels of the network, we make use of four inception modules after the residual layers. The prescription suggested in Szegedy et al. (2015) is to deal with the computational complexity related to the depth of the network, i.e. increasing the size of the network while maintaining the computational cost. The inception module used in this network design is shown in Figure. 4. The idea behind convolving the inputs with a $1\times 1$ filter before the convolutional layers with $3\times 3$ and $5\times 5$ kernels is to reduce the number of feature maps from the inputs as computations are more expensive with larger kernel size. The features at different scales – captured by different kernel sizes $1\times 1$, $3\times 3$ and $5\times 5$ – can be learned simultaneously (Szegedy et al., 2015). It is worth noting that we opt for He initialization (He et al., 2016) for all layers in network II.

For training, as usual for any machine learning tasks, one needs to fine-tune the hyperparameters in order for the algorithm to generalize well and hence achieve the best possible performance, where the distance between the ground truth and network predictions is minimum. To that end, as shown in Table 5 in Appendix B, we use two completely different approaches in terms of optimisation for the two architectures. Although the two networks produce similar results, as will be presented in §5, network I converges faster. This can be explained by the capacity of network I with its number of trainable parameters of about 167 millions which translates to $\sim 2.05\times 10^{9}$ floating point operation per second (flops) at inference time, whereas network II has $\sim 10$ millions of trainable parameters corresponding to $\sim 1.39\times 10^{9}$ flops at inference time.

For reproducibility, we have used TensorFlow package (Abadi et al., 2015) to develop network I which has been trained for 50,000 training steps (about 100 training steps per epoch) which spend $\sim$ 15 hours on a single GPU. Each training step with a batch size of 128 images takes about $\sim$ 1 second. The network converges from the first few epochs but reaches minimum (RMSE$\sim$ 0.001) at epoch = 40 (see Figure 9 in Appendix A). This indicates that same results can be obtained in about 6 hours with a single GPU. For network II, we have used PyTorch (Paszke et al., 2019) resorting to three GPUs to speed up the convergence. Each epoch, in which each GPU processes in parallel a batch of 128 images at a time, takes about 2 minutes which is translated to 40 hours for 1200 epochs.

4 Training dataset

We generate the training dataset from a large simulation box of a size L=150 Mpc with N=2003 cells. We run 1,000 different reionization simulations realizations with 1,000 different random seeds for the initial density field fluctuations. The prior range assumed to our parameters of interest is as follows:

• •

Cosmology:

• –

Matter density parameter: $0.2\leq\Omega_{m}\leq 0.4$.

• –

Hubble constant: $0.6\leq{\it h}\leq 0.8$.

• –

Matter fluctuation amplitude: $0.7\leq\sigma_{8}\leq 0.9$.

• •

Astrophysics:

• –

Photon escape fraction: $0.01\leq f_{\rm esc}\leq 1$.

• –

Rion-Mh power dependence: $0\leq C_{\rm ion}\leq 1$.

• –

Rion redshift evolution index: $0\leq D_{\rm ion}\leq 2$.

The ranges considered for the astrophysical parameters are motivated from our previous MCMC estimates to reproduce various reionization observables (Hassan et al., 2017), and those of the cosmology are inspired by the recent parameters estimates from the Planck Collaboration 2018 (Aghanim et al., 2018). From these priors, we select 1,000 values for each parameter using Latin Hypercube Sampling (LHS) in order to efficiently explore our 6-dimensional parameter space and ensuring that the simulation doesn’t run twice using the same set of parameters. From each simulation run, we store the 21cm brightness temperature at several redshifts $z$ = 10, 9, 8, 7 in order to have a sufficiently large number of maps to ensure training. Balancing the training data set is important to ensure equal learning at all redshifts and all neutral fraction values. This can be achieved by flattening the distribution according to the neutral fraction at each redshift. Flattening the distributions has previously been used in learning cluster masses (Ntampaka et al., 2015) to reduce the bias towards low mass clusters. To flatten our distribution, we take the following steps: First, we ignore highly ionized ($x_{\rm HI}<0.01$) and highly neutral ($x_{\rm HI}>0.99$) 21cm boxes. Second, at each redshift, we bin the boxes according to their neutral fraction. Since the neutral fraction changes strongly with different parameters at different redshifts, the number of boxes in each neutral fraction bin is also different. One has to choose a fixed number to select boxes from bins to flatten the distribution. We here choose 20 boxes. These 20 boxes are randomly selected from each bin. If the neutral fraction bin has less than 20 boxes, then we consider all boxes in this specific bin. If all neutral fraction bins have 20 boxes at the 4 different redshifts, then the total number of all boxes is 800. However, few bins have less than 20 boxes which reduces the total number of boxes to 763. Each selected box has 200 different maps along each of $x,y,z$ directions. This means each redshift has $200\times 3\times 1000=600000$ possible different 21cm maps. However, close maps in the same box would contain similar features, and from our own experience (e.g. Hassan et al., 2019), we have found that $\sim$ 2 Mpc separation between maps is sufficient to obtain distinct maps. To be more conservative, we consider $\sim$ 4 Mpc separation between maps to only select 40 slices along two directions (e.g. x, y) for training, and take 10 slices on the third direction (e.g. z) for testing and validation each. This results in 763 $\times$ ( 40$\times$2 + 10 + 10) = 76,300 total number of images, in which 80% is used for training, 10% for validation and the remaining 10% for testing. In Figure 5, we show the histogram of the training data set as a function of neutral fraction at each redshift. The distribution is approximately flat by construction and includes all possible neutral fraction values at each redshift. We here ignore all current reionization constraints to allow specific neutral fraction values at each redshift such as allowing only high neutral fractions at high redshifts and vice versa. This is an important initial test when constraining parameters, which is to verify the method viability to recover these parameters without imposing any priors and constraints. It is worthwhile to mention that the time-dependence between maps is included through the following:

• •

Each set of the six parameters corresponds to four boxes at redshifts z=10,9,8,7. This shows that the network sees the same six parameters for four different maps from four different redshifts.

• •

The redshift information is encoded in each box through the density field contribution on small scales. This shows that the network sees four different levels of density field contribution in the neutral regions in all maps.

However, an explicit inclusion of this effect is through creating light-cones for each set of the six parameters to account for ionized bubbles growth along the sightline (i.e. the k*∥* modes), redshift-space distortion and angular scales. However, we here assume that all observations at different redshifts are performed independently and apply the noise using a single frequency channel (resolution) corresponding to the map size. The light-cone is more relevant when the full bandwidth is considered. This study sets the baseline for a more detailed analysis to compare the results from 2-dimensional maps (single frequency channel) versus 3-dimensional light-cones (full bandwidth). Indeed, it is expected that more information exists in the reionization window (e.g. Liu et al., 2014) which contains the bubble evolution along the frequency axis through constructing the light-cone. This we leave for future works as its beyond the scope of the current paper.

While our box size, 150 Mpc, might be small to capture the large scale fluctuations and cosmic variance (Iliev et al., 2014), we have previously found that our simulation produces a convergent 21cm power spectrum with respect to the volume (see Figure 8 in Hassan et al., 2016), such that the 150 Mpc volume produces similar power to that from 300 Mpc volume. This indicates the ionized bubbles distribution in 150 Mpc volumes is similar, on average, to those in large volumes. In addition, it has also been found that such a volume produces a convergent reionization history (Iliev et al., 2014). However, the resolution (number of cells) is more important for the neural network performance, since higher resolution maps contain more information and structures. Our maps are composed of 200x200 pixels which are able to resolve the small and large scale fluctuations reasonably well. We leave investigating the network performance in terms of box size and resolution for future works.

5 Results

To assess how well the algorithms perform in terms of predicting the parameters from learning the input features, we use the coefficient of determination, also known as $R^{2}$ score, which is given by

[TABLE]

where $\hat{y}_{i}$, $y_{i}$ and $\bar{y}$ are the predicted value, the actual value and the average of all the actual values in the test sample respectively. The numerator of the second term in Equation 6 – residual sum of squares – quantifies the variation of the predicted values $\hat{y}_{i}$ around the actual values $y_{i}$, and the denominator accounts for the variation of actual values $y_{i}$ around their mean $\bar{y}$. This metric quantifies the strength of the correlation between the inferred and true values of the parameters, in other words unity $R^{2}$ indicates that the network predictions are identical to the ground truth. The $R^{2}$ also quantifies the fraction by which the predicted variance is less than the true variance. For each architecture, we carry out two types of training depending on the input features that the regressors are meant to learn from in order to infer our astrophysical and cosmological parameters (${f_{\rm esc}},{C_{\rm ion}},{D_{\rm ion}},{\Omega_{\rm m}},h,\sigma_{8}$):

• •

feature extraction from a simulated 2D 21cm map (clean/noiseless map hereafter), this involves training and testing using clean maps

• •

feature extraction from a simulated 2D 21cm map which was convolved with a simulated SKA like noise (noisy/mock map hereafter, see §2), this consists of training and testing the networks using noisy maps.

It is worthwhile to mention that we train our networks to predict standardized parameters, meaning that we first subtract the mean and divide by the standard deviation for each array of the parameters. After training, we scale back the predictions to the prior range. Standardizing parameters is important, particularly when the parameters range is different, to prevent the highest parameter range from dominating the loss function. This step is commonly used in multi-parameters regression deep learning tasks.

5.1 Learning from clean maps

The top two rows in Figure 6 show the test results when training the networks with the clean maps. The red and blue areas encompass the 15.9% and 84.1% percentiles (i.e. $\sim$ 1-$\sigma$ level) of the true values given the predictions at each bin for network I and network II respectively. Overall, the constraints on each parameter are very tight. The high value of the $R^{2}$ score ( $\geq 99\%$ for both network I and network II) corresponding to each parameter fitting denotes very strong correlation between the true and inferred parameter, suggesting that the algorithms are able to learn the salient features from the data. On comparing the performance of the two architectures, their $R^{2}$ score for each fitting suggests that they are in a fairly good agreement, and hence perform equally well.

5.2 Learning from noisy maps

The test results after training the algorithms on the noisy maps are presented in the bottom two rows in Figure 6. The constraints on all parameters are slightly weaker as compared to those obtained from training the fitters using the clean maps. The overall decrease in performance denoted by the lower values of $R^{2}$ score corroborates that finding. This can be accounted for by the fact that the relevant features are in this case convolved with noise, therefore extracting them is a bit more challenging.

Despite being convolved with noise, which essentially causes the quality of their features to degrade, all parameters are successfully recovered with an accuracy of $R^{2}\geq 92\%$ for network I and $R^{2}\geq 95\%$ for network II, which is remarkably promising. In contrast to training the algorithms with the clean maps, it can be noticed that, overall, network II outperforms network I, as demonstrated by the $R^{2}$ scores of the former, which are a bit higher than those of the latter on all parameters.

5.3 Dependence on redshift and neutral fraction

In real observations, both foreground contamination and the thermal noise become stronger with increasing redshift. One would then expect some form of dependence of the constraints on redshift. To investigate that possibility, we bin the maps according to their redshift in the test sample and do the predictions by considering each bin separately using the regressors trained with the noisy maps. The results presented in the left panel of Figure 7 suggest the parameter recovery improves with decreasing redshift. While network II tends to have a slightly higher accuracy for each parameter as a function of redshift than network I, the dependence on redshift is fairly mild. This weak dependence is due to the fact that there are all possible neutral fraction values at each redshift, without imposing any prior knowledge to the training dataset by allowing certain neutral fraction values for each redshift, following the current reionization history constraints.

As mentioned and seen earlier, the observed features in a 21cm map are more dependent on $x_{\rm HI}$. To address this effect on the performance of the algorithms, we now bin the slices according their value of $x_{\rm HI}$. It is noticeable in Figure 7 right panel that the performance of each fitter on all parameters declines with increasing value of the neutral fraction. This is expected, as previously seen in Figure 1, the noise always dominates the ionized regions. When the Universe is highly ionized, the prominent features, which are probably the recombining clumps of the remaining dense gas, can still be seen in the presence of noise. This is in contrast to the case where the Universe is highly neutral, and the bubbles are small. Here, the ionized bubbles extend to much smaller scales where the noise dominates, and hence recognition of the bubbles becomes challenging. At this limit, different realizations (with different sets of parameters) of a highly neutral Universe would look similar. This also explains the rapid increase of the accuracy of parameter recovery towards low neutral fraction values. Similar trends have been recently found with using deep learning to constrain the reionization history (e.g. Mangena et al., 2019). It is worthwhile mentioning that this interesting dependence on the neutral fraction cannot be used in future observations since the exact neutral fraction is not known prior observations, although some constraints can be obtained independently from Ly$\alpha$ forest observations (e.g. Fan et al., 2006). However, it is beyond the scope of current networks to use this dependence to constrain parameters. It is rather an interesting theoretical finding and consistent with redshift dependence trends since the Universe is highly ionized at low redshifts.

Having established that the constraints are tighter at lower redshift and lower neutral fraction i.e. ionised case, we now apply some conditions on the test sample as follows

• •

select examples with $x_{\rm HI}<0.5$,

• •

select examples with lower neutral fraction at lower redshift, $x_{\rm HI}<0.5$ and $z<9$.

We show in the top two rows of Figure 8 the resulting constraints on all parameters when considering maps with $x_{\rm HI}<0.5$. The results show how the constraints greatly improve by selecting examples with lower neutral fraction from the test sample. For this specific case, the coefficient of determination value $\geq 0.94$ for network I and $\geq 0.96$ for network II on all parameters indicates that the performance of the algorithms, despite considering mock maps for their training, is comparable to their performance when trained with noiseless maps.

Restricting ourselves to lower $z$, together with only selecting maps with low neutral fraction, in the test dataset further improves the predictions on all parameters as indicated by $R^{2}\geq 0.99$ for both network I and network II (see Figure 8 bottom two rows). Inferring parameters from noisy maps at higher redshift is more challenging, since the noise is stronger (see Figure 7 left panel). Therefore one would expect further improvement of the predictions by combining the two criteria $x_{\rm HI}<0.5$ and $z<0.9$. This is an exciting result for future 21cm surveys that tighter constraints can be obtained from low redshift observations ($z\sim$ 6, 7), where the Universe is highly Ionized. This finding is supported by the fact that the noise is higher at high redshifts, and further confirmed by our additional tests in Appendix B. Using this technique, the SKA will be able to place their first constraints on the astrophysical and cosmological parameters in the near future and from the first cycle of imaging.

5.4 Generalisation error

For the sake of completeness and in order to able to compare our results to other similar studies, we compute, for each parameter, the resulting Root Mean Square Error, rmse, as follows:

[TABLE]

where the summation runs through the whole test dataset. This metric, among others, tells us about the generalisation error inherent to our parameter estimation, in other words the level of accuracy555Not to be confused with accuracy, the metric used in classification tasks. the fitters can achieve on average when estimating the parameters from encoding the inputs. We show the rmse values obtained for each parameter when considering both noiseless and noisy maps in Figure 6 (two top rows and two bottom rows respectively). By comparing the rmse values resulting from training on clean maps and those obtained from training on mock maps, we find that in the idealised scenario the prediction is subject to smaller average error for each parameter in contrast to the realistic one. This finding is expected and consistent with our results based on the $R^{2}$ metric in that the inference is more challenging for each parameter when the data considered for training/testing are contaminated by noise. Although the results based on the two metrics are found to be consistent, it is tempting to expect that for any two different parameters, irrespective of the case (noisy or clean maps), if the $R^{2}$ score of one of them is higher than that of the other, it implies that its rmse must be lower. This trend is seen for all parameters as quoted in the legends.

Gupta et al. (2018), by training a convolutional neural network with $\sim 26$ millions parameters to predict cosmological parameters from simulated noiseless convergence maps, arrived at a generalisation error of $35\times 10^{-3}$ for $\Omega_{{\rm m}}$ and $40.3\times 10^{-3}$ for $\sigma_{8}$. Ribli et al. (2018) improved the constraints with a different neural network architecture of about $\sim 1.4$ millions parameters, by also using simulated lensing maps, achieving rmse = $5.5\times 10^{-3},13.5\times 10^{-3}$ for $\Omega_{{\rm m}}$ and $\sigma_{8}$ respectively. In terms of encoding features from a 2D map using convolutional neural network to infer cosmological parameters, our results – rmse = $5\times 10^{-3}$ and $3\times 10^{-3}$ on $\Omega_{{\rm m}}$ and $\sigma_{8}$ from network I and network II respectively – are comparable to those obtained in these previous works. More importantly, our results corresponding to the realistic case, with/without imposing constraints (see Figure 6 bottom row and Figure 8), where the input maps are noisy are very promising and exciting for future 21cm surveys.

6 Conclusions

We have demonstrated in this work the feasibility of simultaneously inferring both astrophysical and cosmological parameters (${f_{\rm esc},{C_{\rm ion}},{D_{\rm ion}},{\Omega_{\rm m}},h,\sigma_{8}}$) using 21cm maps from the EoR, considering future Hi surveys with the SKA. To this end, we have generated thousands of realizations each with a different set of parameters using simfast21, then compiled a dataset composed of 2D maps (see §4 for details). The approach is to train our two proposed algorithms – convolutional neural network-based – to extract the features from the maps in order to predict the underlying astrophysical and cosmological parameters. We have considered an optimistic case where we train the networks with noiseless simulated maps and a real world-mimicking scenario in which the networks are trained with simulated maps contaminated by simulated SKA-like noise. We have used $R^{2}$ – coefficient of determination – as a performance metric.

We summarize our findings as follows:

• •

The overall results for the idealised case, with $R^{2}$ $\geq 99\%$ for both network I and network II on all parameters, suggest that the algorithms considered in this work are capable of learning the salient features from the maps in order to constrain the corresponding parameters with a remarkably excellent accuracy.

• •

In a more realistic setup, where maps from observations are subject to noise contamination, the constraints on all parameters are slightly weaker, with an accuracy of $R^{2}\geq 92\%$ for network I and $\geq 95\%$ for network II. This is expected since disentangling the relevant information from noise is more challenging. It has been found that network II, leveraging the combination of residual layers at lower level and inception module at higher level of the architecture, outperforms network I despite the former’s lower capacity. This then points towards deploying similar architectures to network II in a real world scenario.

• •

In the case of learning from the noisy maps, the predictions are dependent on both the underlying neutral fraction of the map and its distance from an observer. The performance of the methods improves with decreasing neutral fraction and, as foreground contamination is more important at higher redshift, the constraints are tighter at lower $z$. The results obtained from the test sample is $R^{2}\geq 94\%$ with network I and $R^{2}\geq 96\%$ with network II when only selecting maps with $x_{\rm HI}<0.5$. Recovery improves with imposing constraints on both neutral fraction and redshift ($x_{\rm HI}<0.5$ at lower $z$ $<9$), resulting in $R^{2}>99\%$ with both network I and network II. This indicates that even in the presence of noise in maps, our methods can still estimate the relevant parameters to an excellent level of precision, which is indeed quite promising.

• •

We have computed the prediction error on average – rmse – on each parameter in both optimistic and realistic cases. It has been found that the rmse is smaller in the former, in good agreement with the results when using the coefficient of determination as a performance metric. Compared to other previous works, our approach has also shown a great potential for inferring the underlying parameters of what is observed in future cosmological experiments, such as Hi intensity mapping.

We here considered a redshift range that is consistent with an early reionization scenario, which has been increasingly favoured by galaxy-dominated models of reionization, although more recent work by Kulkarni et al. (2019) shows that galaxies can produce low optical depth and a late reionization scenario. However, late reionization as usually favoured by AGN-dominated scenarios is currently disfavoured (e.g. see Qin et al., 2017; Hassan et al., 2018; Mitra et al., 2018; Parsa et al., 2018). Regardless of the redshift range, the main result, that the accuracy increases with decreasing redshift and neutral fraction, would qualitatively remain valid if lower redshifts are included in this study, such as z=$5,6$, since the instrumental effects are always higher at a higher redshift.

In our analysis, we have generated our training samples based on 1,000 different reionization simulations to constrain 6 parameters. For instance, Gupta et al. (2018), La Plante & Ntampaka (2019), and Gillet et al. (2019) have used 96, 1,000, 10,000 model evaluations to constrain 2, 1, 3 parameters. Schmit & Pritchard (2018) further have shown that 100 model evaluations is sufficient to constrain 3 parameters. In comparison to these works, the number of simulations used in this study is low, which limits the presented results. The prior range assumed in this study is also small (i.e. 0.2 - 2) which places additional limitation to our results. Higher accuracy than reported in this study is expected with larger training samples and more model evaluations, which we will explore in future works.

Our results are entirely limited to the set of assumptions and approximation implemented in our 21cm instrument simulation. A more refined and sophisticated recipe to account for all of the implemented instrumental effects, such as the angular resolution, foreground cleaning and thermal noise, might alter our concluding remarks. The approximation and assumptions implemented in the semi-numerical simulations, through the use of the excursion set formalism to identify the ionized regions, as well as the choice of our dynamic range and resolution, place additional limitations to the presented results. While limited to the SKA, our analysis can be easily extended to include instrumental effects from other 21cm surveys such as HERA and LOFAR, which we leave for future works to perform a detailed comparison between different array designs and different observing strategies. Inferring parameters from the 3D light cones might improve recovery in the presence of noise without the need to impose constraints on the neutral fraction or redshift. Our analysis also can be easily extended to include all of the astrophysical parameters from the source and sink models, and all cosmological parameters, which we leave for future works.

This study has not only highlighted the constraining power of our methods, probing deep into EoR in the near future with the arrival of more advanced Hi instruments like SKA, but also shown how future 21cm surveys and Hi intensity mapping can help break the degeneracy between models by combining them with other experiments, such as Planck, to better the constraints on cosmological parameters in an era of precision cosmology.

Acknowledgements

The authors acknowledge helpful discussions with Laura Boucheron, Kristian Finlator, Jon Holtzman, Tumelo Mangena, James McAteer, Mario Santos, and Rene Walterbos. We particularly thank the anonymous referee for their constructive comments which have improved the paper quality significantly. Simulations and analysis were performed at UWC’s Pumbaa, IDIA/Ilifu cloud computing facilities, and NMSU’s DISCOVERY supercomputers. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562, and computational resources (Bridges) provided through the allocation AST190003P. SA acknowledges financial support from the South African Radio Astronomy Observatory (SARAO) and is grateful to Sean February and Martin Slabber from the Science Data Processing (SDP) team at SARAO for their strong technical support with regards to the computing resources. CCD thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; Their participation in the program has benefited this work. CCD acknowledges funding from the New Mexico Space Grant Consortium Grant #NNX15AL51H.

Appendix A The loss function evolution during training

Figure 9 shows the loss evolution of network I (left) and II (right) for training (blue) and validation (red) samples as a function of training epoch for the case of training on the noisy data set. In both cases, the loss is decreasing as training progresses, which indicates a reduction in the error rate and predictions are approaching the target labels. The fluctuations in the training and validation curves are due to random selection of batches during training. Regardless of these fluctuations, the loss evolution for validation converges and stays constant on average, which indicates that the networks are not over fitting. It is worth noting that the sudden drop in training/validation error while training network II is owing to the fact that the learning rate is updated to 10% of its initial value in order to escape the plateau.

Appendix B Redshift evolution

Our main result which suggests that the accuracy increases with decreasing redshift has been derived from a model trained on mixed maps from all redshifts. To confirm whether accuracy increases towards low redshift, we here perform additional learnings by restricting the training sample to have maps only from the minimum (z=7) or maximum (z=10) redshifts considered in this study (referred to as only). We also compare with predictions at these redshifts from training with the whole dataset, including all other redshifts (referred to as whole), for the case of noisy maps as reported in Table 5. In all cases, we find that the accuracy at z=7 is always higher than that of at z=10. This shows, regardless of training with whole mixed maps or maps at a given redshift, the qualitative trend that the accuracy increases towards low redshifts is still seen as summarized in Table 5. In addition, the quantitative results are also similar with a minimal difference of about $\lesssim$ 2% of accuracy for some parameters as summarised in Table 5 for training with maps at individual redshifts (only) versus those derived using a trained model on mixed maps (whole). Such a minimal difference is expected due to the different number of samples used in the case of “only” versus “whole” tests. This shows that our networks are successful to recover the same qualitative and quantitative results without explicitly including the redshift information as an input to the network (e.g. fitting parameters to four maps from the four different redshifts, $z=10-7$).