HInet: Generating Neutral Hydrogen from Dark Matter with Neural Networks

One of the most popular theoretical techniques used to make mock baryonic simulations is called the halo occupation distribution (HOD). HOD was first used to probabilistically model the number of galaxies residing in a host halo (Seljak 2000; Peacock & Smith 2000; Scoccimarro et al. 2001; Berlind & Weinberg 2002). The HOD technique assumes that the properties of various baryonic structures inside a halo are governed solely by the halo mass, and it ignores all other halo properties. The HOD technique therefore assumes a simple parametric relation between the halo mass and the baryonic properties and uses hydrodynamical simulations (and observations, if available) to calibrate the parameters in this relation. Recent applications of the HOD technique to H i have been in initial field reconstruction and testing the UV background effect on H i maps (Modi et al. 2019a, 2019b) based on the HOD model of Villaescusa-Navarro et al. (2018, hereafter VN18).

However, HOD ignores all environmental effects on the H i abundance and clustering. It is important to note that numerical simulations have shown that cosmological properties like the clustering of halos and galaxies are affected by properties other than halo mass, such as halo environment, halo concentration, spin, and velocity anisotropy (Wechsler et al. 2006; Dalal et al. 2008; Paranjape et al. 2018; Hadzhiyska et al. 2020). This phenomenon is referred to as assembly bias or secondary bias (Sheth & Tormen 2004; Gao et al. 2005).

Other techniques can be used to make mock baryonic simulations, such as subhalo abundance matching (SHAM) and semianalytic models (SAMs). SHAM involves assigning the highest H i mass to the most massive halos and vice versa, but it relies on multiple assumptions, such as that more massive baryonic structures are hosted by the most massive halos and that a monotonic relation, which is free of scatter, exists between masses of baryonic structures and halo masses. SAM, on the other hand, uses a set of simplified equations to model the key baryonic processes in the hydrodynamic simulations (see, e.g., Benson 2012).

It is worth mentioning that, apart from using simulations, there are also proposed perturbative forward-model techniques to evolve the tracer fields directly from linear initial conditions by using biasing schemes (Schmittfull et al. 2019; Modi et al. 2020), or modeling tracer fields by combining biasing schemes and N-body simulations (Sinigaglia et al. 2020).

Convolutional neural networks (CNNs) have been recently applied to numerous areas of physical research (Carleo et al. 2019) and have a lot of potential applications to cosmology (Ravanbakhsh et al. 2017; Modi et al. 2018; Giusarma et al. 2019; He et al. 2019; Ntampaka et al. 2019; Zamudio-Fernandez et al. 2019; Zhang et al. 2019; Kodi Ramanah et al. 2020). An important property of CNNs that makes them useful for cosmological applications is that they are translationally equivariant. In this paper, we use a modified version of a neural network architecture called U-net to generate mock 3D H i fields from a given DM field. Although we focus on H i in the post-reionization regime in this paper, it is worth mentioning some of the recent work on using neural networks to study the H i field during the epoch of reionization (Shimabukuro & Semelin 2017; Chardin et al. 2019; Gillet et al. 2019; Hassan et al. 2020; Kwon et al. 2020; List & Lewis 2020; Mangena et al. 2020; Villanueva-Domingo & Villaescusa-Navarro 2021).

Let us now briefly highlight some of the advantages of neural networks over the above traditional approaches. The HOD formalism typically assumes that the spatial H i density ρ_HI(x) at a particular point inside a halo only depends on the mass of the halo and the distance to its center:

$$\begin{eqnarray}&&{\rho }_{\mathrm{HI}}({\bf{x}})=f({M}_{\mathrm{halo}},| {\bf{x}}-{{\bf{x}}}_{\mathrm{center}}| ).\end{eqnarray} \tag{ 1 }$$

However, VN18 showed that including only the halo mass is not enough for precisely modeling the clustering of the H i field. A more general model for the H i field should also include the information on the environment of the halo and can be roughly intuited as

$$\begin{eqnarray}&&{\rho }_{\mathrm{HI}}({\bf{x}})=g({\rho }_{{\text{}}m}({\bf{x}}),{\rho }_{{\text{}}m}({\bf{x}}^{\prime} )),\end{eqnarray} \tag{ 2 }$$

where ${\rho }_{m}({\bf{x}}^{\prime} )$ is the matter density at points in the neighborhood ¹⁷ of x. Neural networks are universal fitting functions (Hornik et al. 1989) and can be used to accurately approximate the function g in Equation (2); this is the goal of our paper.

The paper is organized as follows. In Section 2, we briefly describe the hydrodynamical simulations that we have used. We then present the specific architecture of our machine-learning model and the method used in Section 3. We discuss the parameters of our benchmark HOD model in Section 4. We present our results in Section 5. Finally, we discuss our results in Section 6 and conclude in Section 7.

The overdensity in the H i field, ${\delta }_{\mathrm{HI}\,}={\rho }_{\mathrm{HI}}/{\bar{\rho }}_{\mathrm{HI}}-1$, varies in the TNG100-1 simulation across ∼9 orders of magnitude. Because the resolution of the TNG100-1 simulation is much higher than the one expected from upcoming surveys, we smooth the H i data with a top-hat filter with a smoothing radius of 300 h⁻¹ kpc. This has a twofold advantage: First, the grid resolution for the H i field is lowered to 140 h⁻¹ kpc, which reduces the size of the data set. Second, the dynamical range over which the H i density field varies is reduced: δ_HI varies over three orders of magnitude, and this reduces the sparsity problem, which we later discuss in Section 3. Note that we did not change the resolution of the input DM field, in order to use as much information in the DM field as possible. Because the training of deep learning models is facilitated when the input data are in the ${ \mathcal O }(1)$ range, we further perform the scaling:

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{\delta }}_{\mathrm{HI}\,}({\bf{x}}) & \equiv & \displaystyle \frac{1}{2}{\left(1+{\delta }_{\mathrm{HI}\,}({\bf{x}})\right)}^{0.2}\\ {\tilde{\delta }}_{\mathrm{DM}}({\bf{x}}) & \equiv & \displaystyle \frac{1}{5}{\left(1+{\delta }_{\mathrm{DM}}({\bf{x}})\right)}^{0.1},\end{array}\end{eqnarray} \tag{ 3 }$$

where δ_HI is the smoothed H i field. The above rescaling gets both ${\tilde{\delta }}_{\mathrm{HI}\,}$ and ${\tilde{\delta }}_{\mathrm{DM}}$ to be nearly in the range [0, 3]. We used a power law instead of a logarithm because the power-law distribution has a flatter tail for high values of δ_HI . We discuss why having a flatter tail is important in Section 3.2.

The deep neural network architecture used in this paper is inspired by the Deep Density Displacement Model (D³M) of He et al. (2019). D³M is the generalization of the standard U-net, which was first proposed by Ronneberger et al. (2015) for use in medical applications. Convolution neural networks, like the U-net, naturally provide properties that are relevant for our problem, such as translational invariance. Variations of the D³M model have been used for large-scale structure applications like learning galaxy modeling and neutrino effects in cosmology (Giusarma et al. 2019; Yip et al. 2019; Zhang et al. 2019). The network architecture we use in this work is shown in Figure 1, and further details are presented in Appendix.

Zoom In Zoom Out Reset image size

Let us discuss an important challenge we face when working with HI data from simulations. The most relevant summary statistics for 21 cm IM, e.g., the HI power spectrum or the HI probability density function (pdf), are dominated by the voxels with the highest HI density. The reason for this is that 21 cm IM is sensitive to the mass-weighted HI, rather than the volume-weighted HI of, e.g., the Lyα forest. Unfortunately, those dense voxels are rare in the simulation. For instance, there are ∼10⁵ halos with M_halo ≥ 10¹⁰ h⁻¹ M_⊙ in our data set, which translates into a very small subset (∼1 in 10³) of voxels of our training sample having a nonnegligible HI density (see the distribution of the voxel H i masses in Figure 7). Because of such sparsity in our data distribution, our model could easily achieve a high accuracy by predicting the low-mass H i voxels, ignoring the high-mass H i voxels; this makes our model harder to train.

Zoom In Zoom Out Reset image size

Modeling the fields in Lagrangian space rather than Eulerian space can in principle reduce the sparsity problem. This is because in Lagrangian space we use the displacement field, which is distributed over a larger region of space, as compared to the density field, which is largely concentrated inside the halo boundaries (He et al. 2019). However, modeling the H i field is not possible in Lagrangian space because, unlike DM, the number of gas particles is not fixed in the simulation. ¹⁹

A similar sparsity challenge exists when predicting the galaxy positions from a 3D DM field. Some of the previous neural-network-based studies have tried to overcome this challenge by using a combination of two neural networks (two-phase model; Modi et al. 2018; Yip et al. 2019; Zhang et al. 2019): the first phase predicts the halo/galaxy position, and the second phase predicts the mass of the halo/ number of galaxies. However, H i, unlike galaxies, is scattered over a wide volume of the universe, not just at the centers of large DM halos. Ignoring the low-mass voxels would remove the Lyα forest from our data, which is not ideal, as the Lyα forest is a powerful cosmological probe by itself, and its contribution to the 21 cm signal at high redshift becomes more important (VN18).

We therefore implement a different kind of two-phase model: the first (second) network is geared toward predicting the low-mass (high-mass) H i voxels. We have used the same U-net architecture shown in Figure 1 for each of the two phases of our model. We provide further details of the two-phase model in Appendix.

As the memory of GPUs is limited, we split the 75 h⁻¹ Mpc TNG100 volume into smaller subboxes. We train the U-net to predict H i boxes of side length ∼1.17 h⁻¹ Mpc and containing 8³ voxels using input DM boxes of side length ∼2.34 h⁻¹ Mpc and 64³ voxels. As illustrated in Figure 1, we want to predict the H i in the subbox residing at the center of a larger DM input. The motivation of using the larger DM box is to account for environmental information for voxels near the boundaries of the HI volume.

We obtain ∼2.5 × 10⁵ nonoverlapping H i subboxes when splitting the TNG100-1 volume. We divide these subboxes into three chunks: ∼60% of the subboxes are used for training the network, 12.5% for validation, and 27.5% for testing the network. We have constructed the test set such that it comprises all the subboxes that correspond to a larger box of side length ∼48 h⁻¹ Mpc. The total number of trainable parameters in the U-net shown in Figure 1 is 2.1 × 10⁷, which, although seemingly gigantic, can be optimized using the technique of automatic differentiation (gradient descent). The gradients are calculated based on the following loss function:

$$\begin{eqnarray}&&{\mathbb{L}}=\sum _{i}^{\mathrm{Voxels}}{\left({\tilde{\delta }}_{i}^{\mathrm{pred}}-{\tilde{\delta }}_{i}^{\mathrm{Illustris}}\right)}^{2}\times \exp [{\left({\tilde{\delta }}_{i}^{\mathrm{Illustris}}-\beta \right)}^{\alpha }],\end{eqnarray} \tag{ 4 }$$

where $\tilde{\delta }$ is the scaled H i density field from Equation (3). Notice that our loss function is different from the traditional mean square error (MSE) loss; we use the additional hyperparameters α and β to add more weight to the high-mass voxels in order to alleviate the aforementioned sparsity problem. ²⁰ We find that α = 2 and β = 0.7 give the best results. We provide further details on the training of the network in Appendix. We have also made use of the techniques of data augmentation and imbalanced sampling to tackle the sparsity problem, and we provide further details in Appendix A.3.

The most widely used summary statistic in cosmology is the power spectrum, which is the Fourier transform of the two-point correlation function. In 21 cm IM, the quantity that is directly observed is the 21 cm power spectrum, which is related to the HI power spectrum via

$$\begin{eqnarray}&&{P}_{21\mathrm{cm}}(k)={\bar{T}}_{b}^{2}\,{P}_{\mathrm{HI}}(k),\end{eqnarray} \tag{ 7 }$$

where ${\bar{T}}_{b}$ is the mean brightness temperature of the 21 cm line at redshift z. Let us now compare the two terms on the right-hand side separately. The mean brightness temperature scales as T_b ∝ Ω_HI, where Ω_HI is the ratio between the density of H i at redshift z and the universe's critical density at z = 0. The network and the HOD predict values of 10⁴ × Ω_HI(z = 1) in our test set to be 5.77 and 6.43, respectively, while that value is 5.82 for the IllustrisTNG simulation.

The second term to compare in Equation (7) is P_HI(k), and we show the results from different methods in Figure 3. We find that the network is able to reproduce the HI power spectrum from the simulations up to a deviation of ≲5% all the way to nonlinear scales k ≲ 1 h Mpc⁻¹, whereas the HOD deviates by ≲20% at low k, although it becomes slightly more accurate at high k. The increase in the HOD accuracy at high k arises mainly as a result of the one-halo term given in Equation (6). It is important, however, to note that the one-halo term in Equation (6) only aids in the accuracy of the power spectrum (amplitude of fluctuations) at high k but does not accurately model the phases of the H i fluctuations, which are relevant for the cross-correlation coefficient r_HI (compared in the right panel of Figure 3). VN18 argued that one possible cause of the discrepancy in the HOD power spectrum at low k may be the fact that they do not explicitly fit the total Ω_HI in the simulation volume, which would change the H i bias. However, even if the bias is changed, their r_HI should remain unchanged.

We do not show error bars arising from sample variance in the testing volume because both the DM and the H i fields from IllustrisTNG are evolved from the same initial conditions. We chose the particular range of k in Figure 3 because of the two following constraints: the low-k limit is set by the largest mode in the test set (making an even larger test set is possible, but at the expense of reduction in the training data), and the high-k limit is conservatively set such that the scales affected by the smoothing of the H i maps (∼5 × 300 h⁻¹ kpc) are removed (see Figure 9 for P_HI plotted until k = 10 h Mpc⁻¹).

Zoom In Zoom Out Reset image size

It is worth noting that the accuracy of P_HI(k) is dependent on modeling of H i in high-mass halos (see, e.g., Figure 8). As mentioned in Section 3.2, the training of our network is challenging in the high-mass end, because of the sparsity associated with the abundance of voxels in that regime. Our results can be further improved if we train the U-net on a simulation with a larger volume, ²¹ which would have a larger number of high-mass halos, or using a set of zoom-in simulations focused on galaxy clusters (Thiele et al. 2020).

Zoom In Zoom Out Reset image size

Let us now estimate, very crudely, the accuracy of the power spectrum model required for future 21 cm surveys using some fisher forecasts in the literature. Obuljen et al. (2018) showed that the peak signal-to-noise ratio (S/N) could be ∼10 for upcoming surveys like HIRAX and CHIME (note that S/N will be much larger for future surveys like SKA), and therefore the accuracy of the method to produce 21 cm mocks should roughly be ≲10% for the power spectrum.

Large regions of future H i surveys will overlap with regions sampled by galaxy spectroscopic surveys like DESI or Euclid. One of the most important summary statistics in such overlapping regions is the H i−galaxy cross-power P_{HI − Galaxy}. Using the cross-correlations will boost the S/N and mitigate the effects of foregrounds (e.g., Villaescusa-Navarro et al. 2015b predict the S/N for the cross-correlation between future 21 cm surveys and Lyman break galaxies to be larger than 21 cm autocorrelation by a factor of ∼10 (∼4) at large (small) scales; see also Padmanabhan et al. 2020; Modi et al. 2021). Furthermore, unlike the auto-power spectrum of 21 cm, which is yet to be detected, there have already been multiple detections of the P_{HI − Galaxy} signal (Chang et al. 2010; Masui et al. 2013; Anderson et al. 2018).

In Figure 4, we compare P_{HI − Galaxy} from our approach. We have included all the galaxies in the TNG100-1 sample with the stellar mass M_* > 10¹⁰ M_⊙ (which roughly corresponds to a number density of n = 10⁻³ h³ Mpc³) for this measurement. Similar to the P_HI(k) in Figure 3, the P_{HI − Galaxy} exhibits a bias for the HOD at low k.

Zoom In Zoom Out Reset image size

Another important statistic along the same lines is the cross-correlation coefficient of H i with halos, which is more sensitive to the way H i mass is distributed across halos and to the one-halo term. We see in Figure 5 that the U-net outperforms the HOD up to k ≲ 1 h Mpc⁻¹ for all bins of halo masses that we have considered. It is worthwhile to note that if we extend Figures 4 and 5 for k > 4 h Mpc⁻¹, the green HOD curve diverges farther away from the orange IllustrisTNG curve. It is therefore likely a coincidence that the HOD outperforms the U-net for 1 h Mpc⁻¹ ≲ k ≲ 4 h Mpc⁻¹, as the green curve is merely crossing the orange curve to the other side. However, we do not show the scales k > 4 h Mpc⁻¹, as they are affected by the smoothing of H i maps and a high-resolution test has to be performed to be definitive.

Zoom In Zoom Out Reset image size

If all the fluctuations in the universe were perfectly Gaussian, the field could be perfectly characterized by its two-point correlation function or its power spectrum. However, even if the primordial fluctuations were Gaussian, late-time gravitational clustering causes significant leakage of Gaussian information in the nonlinear regime (Scoccimarro et al. 1999; Takada & Jain 2004; Villaescusa-Navarro et al. 2020; Wadekar & Scoccimarro 2020). To recover this information, the lowest-order statistic that one needs to compute in Fourier space is the bispectrum. Unlike the power spectrum, the bispectrum is sensitive to the shape of the structures generated by gravitational instability and has the promise to break important degeneracies in the bias and cosmological parameters (Scoccimarro 2000; Sefusatti et al. 2006; Chudaykin & Ivanov 2019; Yankelevich & Porciani 2019; Hahn et al. 2020; Hahn & Villaescusa-Navarro 2021; Kamalinejad & Slepian 2020). The post-reionization 21 cm signal is expected to have significant information in the nonlinear regime, and there have been recent attempts at theoretical modeling of the H i bispectrum (Sarkar et al. 2019). One of the toughest parts in a bispectrum analysis is calculation of the error due to cosmic variance, and a fast technique to generate mock H i fields is therefore essential.

In Figure 6 we show the HI bispectrum of the IllustrisTNG simulation, together with the predictions of the HOD and neural network. We show two particular cases (which are representative of all possible triangle configurations): the first case is for equilateral triangles with different side lengths, and the second case is for triangles with a varying angle between two of its sides whose lengths are kept fixed. We checked that the results are similar for other triangle configurations. The residuals of the bispectrum of equilateral triangles from the network compared to the target are ≲20% for all scales we compared (0.2 h Mpc⁻¹ < k < 4 h Mpc⁻¹), and the residuals for HOD, on the other hand, are ≲45%.

Zoom In Zoom Out Reset image size

The upcoming H i surveys will observe systems ranging from low column densities (Lyα forest) to very high column densities (damped Lyman absorbers); the pdf of H i is a statistic that is sensitive to its distribution over the wide density range. The H i pdf can also be used for probing non-Gaussian information and for constraining luminosity functions (Breysse et al. 2017; Leicht et al. 2019). We show in Figure 7 the comparison of the U-net prediction over four orders of magnitude of H i voxel masses. The comparison is difficult in the high-mass end owing to sample variance: that regime is dominated by very massive halos, which are rare. Note that there is a slight discrepancy in the HOD approach for low H i mass voxels, as it does not take into account filaments and also misses very low mass halos that are below the simulation resolution threshold.

Voids are the most underdense regions of the universe. In Figure 7 we show the void size function (VSF) of the HI field, which is defined as the number density of HI voids as a function of radius. We have used the algorithm described in Banerjee & Dalal (2016) to identify voids. The VSF is an important statistic, as it contains complementary information to the one from traditional clustering observables.

We presented our network architecture in Figure 1, and we discuss its details in this section. Let us first discuss details of the input and the output of our network. The input DM box has a physical side length 2.34 h⁻¹ Mpc with a grid resolution of ∼35 h⁻¹ kpc and therefore has 64³ grid cells. As discussed in Section 2.1, the output H i box has a lower resolution (physical side length of 1.17 h⁻¹ Mpc with 8³ grid cells). A U-net typically consists of a contracting path and an expansive path of nearly equal lengths. In our case, due to lower dimensionality of the output as compared to the input, the expansive path is relatively much shorter.

The contracting path follows the typical architecture of a convolutional neural network and consists of four primary blocks. Each block consists of two successive convolutions with stride 1 and a down-sampling convolution with stride 2, each of the which is followed by batch normalization (BN) and a rectified linear unit (ReLU). For each convolutional layer, we use 3 × 3 × 3 filters. As our input DM fields are not periodic, we cannot use periodic padding similar to He et al. (2019); we instead apply a zero padding with size 1. At each down-sampling step we double the number of feature channels and reduce the number of grid points along each dimension by a factor of 2, and vice versa for the up-sampling step.

We only have one up-sampling step in our architecture. We concatenate the up-sampled map (dimensionality 8³) with two maps from the contracting path (denoted with filled boxes): one of the same dimensionality (8³) and another of a higher dimensionality (16³), but after application of a max-pooling layer. These concatenations help to train the network faster as the gradients are passed through the transverse connections. Concatenation also provides the network with various levels of granularity for the final prediction. We have also tried using a deeper neural network (which can be obtained by using stride 1 instead of stride 2 while moving between layers in Figure 1), but we did not find an improvement in our results.

As discussed earlier in Section 3.2, we needed to employ a two-phase model to tackle the problem of data sparsity. Let us now discuss in detail how we employ the two neural networks (labeled as FP (SP) for the first (second) phase hereafter). A simple way to think about the functions of the two networks is the following: FP fills the low H i mass voxels and identifies boundaries of large H i halos. SP is then used to assign an appropriate H i mass profile to the large halos.

Both FP and SP have the same U-net architecture that was shown in Figure 1. If we only use FP in our prediction, we get the results shown by a dashed line in Figure 8. From the pdf, we see that the FP is reproducing the pdf accurately for all H i voxels except the high-mass ones. We believe that this is due to a dearth of training data on the high-mass H i voxels and expect the FP to perform better with more data. To complement FP, we choose to employ an SP that is focused on predicting only the high-mass H i voxels. The training of the SP is therefore a little different from that of the FP. Equation (3) was used to generate rescaled fields to train the FP, but for SP, we use different rescaling transformations:

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{\delta }}_{\mathrm{HI}\,}^{(2)} & = & 0.03\,{\left(1+{\delta }_{\mathrm{HI}\,}\right)}^{0.5}\\ {\tilde{\delta }}_{\mathrm{DM}}^{(2)} & = & 0.0035\,{\left(1+{\delta }_{\mathrm{DM}}\right)}^{0.25}.\end{array}\end{eqnarray} \tag{ A1 }$$

This is because the high-mass tails of the rescaled DM and H i distributions are flatter for the transforms in Equation (A1) as compared to Equation (3), which makes the SP easier to train in the high-mass regime. We had used the loss function in Equation (4) for training the FP. The loss function for training the SP is obtained by substituting $\tilde{\delta }\to {\tilde{\delta }}^{(2)}$ in Equation (4), and we find that the hyperparameter values β = 0.26 and α = 2.5 give the best results for training the SP.

Let us now discuss combining the predictions of FP and SP after they are trained separately. For a given DM input field, let us denote their unscaled prediction of a particular voxel of the output H i field by δ^FP and δ^SP. We get our final result for that particular voxel by combining the predictions using weights: w^FP δ^FP+w^SP δ^SP; the weights are given by

$$\begin{eqnarray}\begin{array}{rcl}{w}^{\mathrm{FP}} & \equiv & \displaystyle \frac{1}{1+\exp \left[({\delta }^{\mathrm{SP}}-{\delta }^{\mathrm{threshold}})/100\right]}\\ {w}^{\mathrm{SP}} & \equiv & 1-{w}^{\mathrm{FP}},\end{array}\end{eqnarray} \tag{ A2 }$$

and they were chosen such that the output from SP is given more importance for high-mass voxels and vice versa. We have adopted δ^threshold = 1500 for results in this paper, but we checked that our results are not sensitive to the threshold in the range δ^threshold ∈ (1000 , 2000).

For training the network, we used the Adam optimizer (Kingma & Ba 2014) with a learning rate ranging between 10⁻⁸ and 10⁻⁴. The network takes nearly 30 epochs to train, and we choose to use a batch size of 28, as it was compatible with the memory constraints of our GPUs.

We also train the neural network to recognize polyhedral symmetries like translational and rotational invariance (He et al. 2019). We do this by generating multiple instances of a given input DM box by applying various symmetry transformations and then training the U-net on all the generated cases. This process is often referred to as data augmentation because we are generating multiple training simulations given using a single data realization. We use all transformations corresponding to the symmetry group of a cube that consists of 48 elements (octahedral group of order 24 with an additional factor of 2 to account for inversion across the origin ($\vec{r}\to -\vec{r}$)). As the cubes containing high H i densities are rare in our training set, using data augmentation greatly helps for training in the high-density regime. We also use a popular technique for tackling data sparsity called nonuniform sampling, where we input the cubes containing high H i density voxels more frequently during training. In our case, we first rank the cubes with respect to the maximum H i density inside them and then increase the frequency of occurrence of the top 1% of cubes 100-fold. However, even after using data augmentation and nonuniform sampling, we still had to use the two-phase neural network approach, as the data sparsity problem in our case is quite extreme.

We trained all of our models using New York University's High Performance Cluster using NVIDIA Tesla P-100 and P-40 GPUs.

In Section 5 we showed results from the test set H i cube. To construct this box, we used the U-net to generate small 1.17 h⁻¹ Mpc side H i subcubes and then stacked them together to produce the test 48 h⁻¹ Mpc cube. One possible concern might be that such stacking produces edge effects that will affect summary statistics at scales comparable to the length of the box (k ∼ 5 h Mpc⁻¹). We therefore check that there indeed are no artifacts on these scales and show the high-k H i power spectrum comparison in Figure 9. Furthermore, we also show a zoom-in version of the output of the test cube in Figure 10. A reason for the absence of edge artifacts could be the larger size of DM box used to predict a smaller H i box (see Figure 1), so information beyond the location of edges of the H i box is included in the U-net prediction. We also tried using the technique of annealing when concatenating the small H i subcubes to smoothen their edges, but our results did not change.

Zoom In Zoom Out Reset image size