Deep learning for intensity mapping observations: Component extraction

We generate a number of mock intensity maps for training and testing in the following manner. First, we populate a cubic volume of $280h^{-1}~{}\rm Mpc$ with dark matter haloes using the publicly available pinocchio code (Monaco et al., 2013). We set the minimum halo mass of the catalog to be $3\times 10^{10}h^{-1}~{}\rm M_{\odot}$. We have tested and confirmed that the map properties such as the total line intensity are not significantly affected by this choice of the minimum halo mass.

We derive the halo mass-luminosity relation using the outputs of the cosmological hydrodynamics simulation Illustris-TNG (Nelson et al., 2019). We use the TNG300-1 dataset which has a simulated volume of $V_{\rm box}=(302.6~{}\rm cMpc)^{3}$. We compute the line luminosity from a simulated galaxy as

where $A_{\rm line}$ accounts for attenuation by dust. We adopt $A_{\rm H\alpha}=1.0$ mag and $A_{\rm[O{\sc iii}]}=1.35$ mag. We use the photoionization simulation code cloudy (Ferland et al., 2017) to compute the coefficient $C_{\rm line}(Z)$ as a function of the mean metallicity of the galaxy. The cloudy computation is done in the same manner as in Moriwaki et al. (2018) except that we adopt typical values of the electron density $n=100~{}\rm cm^{-3}$ and the ionization parameter $U=0.01$. We compute the luminosity of a simulated halo by summing up the luminosities of its member galaxies. The halo mass-luminosity relation, i.e., the mean $\overline{L}_{i}$ and the variance $\sigma_{i}$ within each halo mass bin, is then obtained from the Illustris output. To generate an emissivity field from a pinocchio halo catalogue, we assume a Gaussian distribution with a mean $\overline{L}_{i}$ and a variance $\sigma_{i}$ and assign luminosities to the haloes in $i$-th mass bin.

We perform the above procedure for H$\rm\alpha$ and [Oiii] indenpendently at the respective redshift. We generate two-dimensional intensity maps by projecting the three-dimensional emissivity fields along one direction. The total area subtended by a map is $(3.4~{}\rm deg)^{2}$, and we assume a spectral resolution $R=40$ that corresponds to the expected resolution of SPHEREx. We find that the relative contribution from the [Oiii] emission (at $z=2.0$) is $\sim 60\%$ of the H$\rm\alpha$ map (at $=1.3$), which is consistent with other theoretical studies (Fonseca et al., 2017; Silva et al., 2018). To generate a large number of training data set, we use 300 different halo catalogues.

For each realization, 100 maps with an area of $(1.7~{}\rm deg)^{2}$ are generated by projecting along random direction. We obtain 30,000 training data in total. In this way, we obtain training maps with various mean intensities. Each map has $256\times 256$ pixels, corresponding to a pixel size of $(0.4~{}\rm arcmin)^{2}$. For the test data set, we produce another 1,000 halo catalogs and generate 1000 independent maps. We smooth the training and test maps with a Gaussian beam with $\sigma=1.2~{}\rm arcmin$ before giving them to the networks.¹¹1We have performed the same analysis in the present paper to the ”raw” data without smoothing. We have found the performance of the cGANs is somewhat degraded when the networks are trained with the unsmoothed data populated with a number of discrete sources.

We construct cGANs based on pix2pix by Isola et al. (2016).²²2 https://github.com/yenchenlin/pix2pix-tensorflow We train the networks so that they reconstruct both H$\rm\alpha$ and [Oiii] images from an observed image. This kinds of one-to-many image translation networks are studied by, for instance, Lee, Yang & Oh (2018) for separating transparent and reflection scenes.

We have two pairs of adversarial convolutional networks called generator and discriminator. They are denoted by $(G_{1},D_{1})$ for H$\rm\alpha$ map and by $(G_{2},D_{2})$ for [Oiii] map. The generators try to reconstruct H$\rm\alpha$ and [Oiii] maps from an observed map $X_{\rm obs}$, whereas the discriminators try to distinguish the true maps $X_{\rm true,i}$ and the reconstructed maps $G_{i}(X_{\rm obs})$. In other words, for an input $(X_{\rm obs},X)$ with $X$ being either $X_{\rm true,i}$ or $G_{i}(X_{\rm obs})$, the discriminator returns a probability that $X$ is $X_{\rm true,i}$. Here, $X_{{\rm true},i}(i=1,2)$ denote the true maps of H$\rm\alpha$ and [Oiii], respectively. The generator consists of 8 convolution layers followed by 8 de-convolution layers and the discriminator consists of 4 convolution layers. Two generators $G_{1}$ and $G_{2}$ share the first 8 convolution layers. The kernel size of the convolutions is $5\times 5$. In each layer, batch normalization,³³3 During test phase, we set is_training = False in batch normalization to use fixed normalization parameters. dropout, and skip connection are also performed (see Isola et al., 2016, for more details).

During the training phase, the performance of the generators and the discriminators are evaluated by a linear combination of the cross-entropy losses and the mean L1 norms:

where

where $N_{\rm pix}=(256)^{2}$. In each training set, the generators (discriminators) are updated to decrease (increase) the loss function $\mathcal{L}$ averaged over a mini batch.⁴⁴4 Mini batch is a randomly selected set of training data, $\{X_{{\rm obs},i},X_{{\rm true},i}\}_{i=0}^{n_{\rm b}}$, where $n_{\rm b}$ is batch size. In training phase, the networks pass through all the training data without duplication. When we set the number of epochs $n_{\rm e}>1$, this passing through is repeated for $n_{\rm e}$ times. For $n_{\rm d}$ training data, updates are performed for $n_{\rm d}n_{\rm e}/n_{\rm b}$ times in total. We adopt $\lambda_{1}=\lambda_{2}=\lambda_{\rm tot}=100$ and a batch size of 4. The networks are trained for 8 epochs. We use the Adam optimizer (Kingma & Ba, 2014) with learning rate 0.0002, and decay rate parameters $\beta_{1}=0.5$ and $\beta_{2}=0.999$ for updating the parameters.

We study the performance of our networks with 1000 test data. Fig. 1 shows an example of true and reconstructed maps. In our fiducial case of $\lambda_{\rm obs}=1.5\mu$m, the contribution from H$\rm\alpha$ is larger than [Oiii], and then outstanding structures in the observed map mostly originate from the H$\rm\alpha$ emission at $z=1.3$. It is thus remarkable that not only the H$\rm\alpha$ distribution but also the weaker [Oiii] intensity is reproduced well.

It is important to study whether statistical quantities are also reproduced accurately. We first examine the peaks in our intensity maps. We select as ”peaks” local maxima with heights greater than $3\sigma$. We find 24089 (18859) and 24800 (17631) peaks in the true and the reconstructed H$\rm\alpha$ ([Oiii]) maps over our 1000 test data sets. Among them, 18262 (5095) peaks are matched correctly. This means that 76% (27%) of the true peaks are reproduced, and 74% (29%) of the reconstructed peaks are true.

If our purpose is to study individual peaks or other individual structures, we may require much higher accuracy for reconstructed structures. Previous studies developed cGANs that also learn the reliability of reconstructed maps (Lee, Yang & Oh, 2018; Kendall & Gal, 2017). In principle, we can use these methods to quantify the reliability of the outputs. Another promising idea is to combine multiple networks. To test this idea, we use 5 networks that have an identical set of convolutional layers but are trained with different sets of data. Intensity maps reconstructed by the 5 trained networks are similar, but not exactly the same. We find that it is generally difficult to reproduce the true intensity in portions where these networks commonly fail. For H$\rm\alpha$ ([Oiii]) maps, the number of peaks detected by all the 5 networks is 14332 (895). Among them, 13018 (539) peaks are true, which means a 91% (60%) confidence level for our peak detection. We note that if we take the average or the median of the reconstructed maps by multiple networks on a pixel-by-pixel basis, dark structures in void regions and small-scale structures are smoothed out.

Summary statistics such as the mean intensity and the power spectrum are primary tools to study the distribution and the properties of the emission-line galaxies. These can then be used for galaxy population studies or for cosmological parameter inference. In this section, we examine how well the mean intensity and the power spectrum are reconstructed. We take medians of the reconstructed statistics by 5 different networks and compare them with true ones.

Fig. 2 shows the correspondence of the true and the reconstructed mean intensities. The mean intensities are widely distributed because of the cosmological variance of the underlying density field. We see clear correlations between the true and the reconstructed mean intensities. Fig. 2 shows that the mean H$\rm\alpha$ ([Oiii]) intensity of each $(1.7~{}\rm deg)^{2}$ map can be estimated with $\sim 10$% ($\sim 20$%) accuracy. We note that the residual uncertainty is comparable to those resulting from the luminosity function estimates by recent galaxy surveys (Sobral et al., 2013; Khostovan et al., 2015). Planned LIM observations would have a much larger observational area. For instance, SPHEREx (Doré et al., 2016) will perform a deep survey over $200~{}\rm deg^{2}$ and thus the estimated mean intensities would have a much smaller statistical uncertainty.

We test if our networks generate accurate images (intensity maps) if the input observed map is significantly different from the training data. To this end, we input intensity maps with the mean differing as much as 20 %. Some of these samples have mean intensities below or above the range plotted in Fig. 2. We find that the networks reconstruct the H$\rm\alpha$ and [Oiii] intensities with accuracy similar to those shown in Fig. 2 when both maps are scaled with the same factor. However, we find systematic offset when only one map is scaled more than 10% while the other being unchanged. In order to reconstruct these ”outliers” accurately, we need to consider a wide variety of training data, and/or to combine multiple networks trained with maps in different mean intensity ranges.

Another important statistic is two-dimensional power spectrum. Fig. 3 shows the variation of the true (shaded regions) and the reconstructed (error bars) power spectra. For reference, we also show the power spectra of unsmoothed true maps.⁵⁵5 Power spectrum of the unsmoothed map $P(k)$ can be recovered from that of the smoothed map $P_{\rm sm}(k)$ by $P(k)=\exp(k^{2}\sigma^{2})P_{\rm sm}(k)$, where $\sigma$ is the smoothing scale of the Gaussian beam. Clearly, our networks learn the clustering of galaxies even though we do not explicitly teach that galaxies at different redshifts have different clustering amplitudes. The top panel of Fig. 3 shows the difference between the true and the reconstructed power spectra normalized by the square root of the variance of the true power spectra $\sigma_{\rm true}$. We note that the variance of the training data is also $\sigma_{\rm true}$. For H$\rm\alpha$ map, the difference is typically less than $\sigma_{\rm true}$ at large scales; our network is able to recover the power spectrum of H$\rm\alpha$ at $z=1.3$ with an accuracy of $\sim 10\%$ from a confused map.