Simulating a full-sky high resolution Galactic synchrotron spectral
  index map using neural networks
We present a model for the full-sky diffuse Galactic synchrotron spectral
index with an appropriate level of spatial structure for a resolution of 56
arcmin (to match the resolution of the Haslam 408 MHz data). Observational data
at 408 MHz and 23 GHz have been used to provide spectral indices at a
resolution of 5 degrees. In this work we make use of convolutional neural
networks to provide a realistic proxy for the higher resolution information, in
place of the genuine structure. Our deep learning algorithm has been trained
using 14.4 arcmin observational data from the 1.4 GHz Parkes radio continuum
survey. We compare synchrotron emission maps constructed by extrapolating the
Haslam data using various spectral index maps, of different angular resolution,
with the Global Sky Model. We add these foreground maps to a total emission
model for a 21 cm intensity mapping experiment, then attempt to remove the
foregrounds. The different models all display different spectral or spatial
behaviour and so each provide a useful and different tool to the community for
testing component separation techniques. We find that for an experiment
operating using a cosine aperture taper beam with a primary Full Width at Half
Maximum between 1.1 and 1.6 degrees, and the principal component analysis
technique of foreground removal, there is a discernible difference between
synchrotron spectral index models with a resolution larger than 5 degrees but
that no greater resolution than 5 degrees is required.
M. O. Irfan
# Simulating a full-sky high resolution Galactic synchrotron spectral index map using neural networks

###### Abstract

We present a model for the full-sky diffuse Galactic synchrotron spectral index with an appropriate level of spatial structure for a resolution of 56 arcmin (to match the resolution of the Haslam 408 MHz data). Observational data at 408 MHz and 23 GHz have been used to provide spectral indices at a resolution of 5 degrees. In this work we make use of convolutional neural networks to provide a realistic proxy for the higher resolution information, in place of the genuine structure. Our deep learning algorithm has been trained using 14.4 arcmin observational data from the 1.4 GHz Parkes radio continuum survey. We compare synchrotron emission maps constructed by extrapolating the Haslam data using various spectral index maps, of different angular resolution, with the Global Sky Model. We add these foreground maps to a total emission model for a 21 cm intensity mapping experiment, then attempt to remove the foregrounds. The different models all display different spectral or spatial behaviour and so each provide a useful and different tool to the community for testing component separation techniques. We find that for an experiment operating using a cosine aperture taper beam with a primary Full Width at Half Maximum between 1.1 and 1.6 degrees, and the principal component analysis technique of foreground removal, there is a discernible difference between synchrotron spectral index models with a resolution larger than 5 degrees but that no greater resolution than 5 degrees is required.

keywords: Cosmology: diffuse radiation, Methods: statistical, Radio continuum: ISM

## 1 Introduction

Component separation has proven fundamental to observational cosmology; disentangling diffuse Galactic foregrounds from a cosmological signal of interest has been a central theme for Cosmic Microwave Background studies for decades (Bennett et al., 2003; Gold et al., 2009, 2011; Leach et al., 2008; Delabrouille et al., 2013; Planck Collaboration et al., 2014, 2016, 2020). More recently a plethora of low-frequency (\(<1.5\) GHz) radio cosmology experiments have started observing with the aim of measuring the redshifted 21cm hydrogen line in order to probe Cosmic Dawn (Eastwood et al., 2019; DeBoer et al., 2017), the Epoch of Reionisation (Parsons et al., 2010; Tingay et al., 2013) or the formation of large scale structure (Newburgh et al., 2014; Battye et al., 2012; Nan et al., 2011; Santos et al., 2016). At these frequencies diffuse Galactic synchrotron emission dwarfs the cosmological signal of interest; this emission is typically modelled as a power law with a spectral index which scales the temperature across frequency. The synchrotron spectral index changes both spatially and gradually across frequency. Numerous works have investigated the mitigation of simulated foreground contamination on the detection of the simulated Hi signal, e.g. Wolz et al. (2014); Shaw et al. (2015); Bigot-Sazy et al. (2015); Alonso et al. (2015); Chapman et al. (2016); Zhang et al. (2016); Mertens et al. (2018); Carucci et al. (2020); Licardo et al. (2011); Cunnington et al. (2021); Irfan & Bull (2021); Makinen et al. (2021); Yohana et al. (2021); Soares et al. (2022); Spinelli et al. (2022). It is vital that these simulated foregrounds contain an accurate level of spatial and spectral complexity to prevent a misleading simplification of the component separation problem.

Publicly available repositories of diffuse Galactic foreground models are of enormous use to the community as they provide a test-bed on which to assess the qualities and deficiencies of component separation methods. Such resources include the Global Sky Model (Zheng et al., 2017), which uses observational data to produce all-sky maps of diffuse Galactic emission between 10 MHz and 5 THz. The Global Sky Model (GSM) performs principal component analysis on empirical data sets to determine the statistically independent components within the sky maps and then interpolates this information to model the total diffuse emission temperature at any frequency within the 10 MHz to 5 THz range. Additionally, there are the Planck
Sky (Delabrouille et al., 2013) and Python Sky (Thorne et al., 2017) models which model the different physical components expected in the Galaxy due to differing emission mechanisms and use spectral index information to scale these maps across frequency.

The 408 MHz all-sky map of Haslam et al. (1982) is typically taken as a proxy for an all-sky map of synchrotron emission because synchrotron is thought to be the dominant diffuse emission at this frequency across the majority of the sky (excluding the central Galactic plane and a few specific molecular clouds). The simplest way to scale the Haslam map from 408 MHz to any other frequency (\(\nu\)) is to use the power law parametrisation:

\[T_{\nu}(p)=T_{408}(p)\,\left(\frac{\nu}{408}\right)^{\beta}\,, \tag{1}\]

and assume a single value for the spectral index (\(\beta\)) which remains constant across pixels (p). However, the synchrotron spectral index is known to vary spatially due to energy losses of the relativistic, charged particles responsible for the emission (Bennett et al., 2003). Previous works have highlighted the need to consider a variable synchrotron spectral index for the problem of 21 cm foreground removal; demonstrating that the combination of a frequency-changing beam plus a spatially-varying spectral index is a more challenging foreground removal problem than just the situation of a frequency-changing beam plus a constant spectral index (Bernardi et al., 2015; Mozdzen et al., 2016; Anstey et al., 2021).

Miville-Deschenes et al. (2008) used the 408 MHz all-sky map together with the _WMAP_ 23 GHz polarisation map to determine an all-sky model for the synchrotron spectral index between 0.408 and 23 GHz. The 23 GHz map of polarisation intensity was specifically used as at GHz frequencies synchrotron emission is no longer the dominant diffuse Galactic emission in intensity, but it is the only emission expected to be present in polarised intensity. To increase the signal-to-noise ratio of the polarised data the maps were smoothed to 5 degrees, enabling the production of an all-sky synchrotron spectral index map at 5 degree resolution. This map is used in both the Planck Sky Model and Python Sky Model to scale the 408 MHz Haslam map across different frequencies. While it is generally accepted that the synchrotron spectral index changes spatially, the relationship between spectral index steepness and Galactic latitude is still not fully understood. Some spectral index maps, such as the 45 to 408 MHz Guzman et al. (2011) map, show a range of values with the shallowest indices occurring within the Galactic plane. Dickinson et al. (2009) display six different spectral index maps from the literature, all with very different spatial features, some of which display a trend in spectral index steepness with Galactic latitude while some do not, and discuss how the visible spatial features depend on the variation assumptions made for the spectral index calculation. Guzman et al. (2011) attribute the presence of shallow spectral indices across the Galactic plane to free-free emission (spectral index \(\sim-2.1\)) due to the increased amount of warm ionised hydrogen within the Galactic plane. Miville-Deschenes et al. (2008) actually present three distinct maps of the synchrotron spectral index. Two are made from the Haslam 0.408 and _WMAP_ 23 GHz intensity maps; at 23 GHz in intensity several diffuse Galactic emissions are present across the sky. While model 1 models the free-free contribution at 23 GHz, model 2 models both the free-free and anomalous microwave emission contributions. Both of these models also show shallower indices running across the plane. The third model, however, is made using the Haslam 408 MHz intensity map and the _WMAP_ 23 GHz polarised intensity map. As synchrotron emission is believed to be the only non-negligible emission present in polarised intensity at 23 GHz, this is the only spectral index model that doesn't rely on a free-free or anomalous microwave emission model; this is also the only model out of the three which does not display shallow spectral indices within the Galactic plane. It is this spectral index model which is used in the Planck Sky Model and the Python Sky Model, therefore we follow suite and also use this spectral index map as our basis.

Previous works have already demonstrated the need for a spatially complex synchrotron spectral index model and provided such models; the highest resolution model available being at 5 degrees. We aim to determine if this level of spatial accuracy is enough or if a spectral index map with the same resolution as the Haslam data themselves would further complicate the foreground removal process. Inspired by the recent success of Krachmalnicoff and Puglisi (2021) who used Convolutional Neural Networks (CNNs) to learn high resolution features of thermal dust models from low resolution input models, we aim to simulate high resolution information for the Miville-Deschenes et al. (2008) spectral index map. We use a CNN trained on a spectral index map constructed using both Haslam and Parkes (CHIPASS) (Calabretta et al., 2014) observational data. CNNs are a deep learning technique often employed for the task of image segmentation (assigning a label to each pixel of an image). We do not attempt a physically motivated model of synchrotron emission, as in Waelkens et al. (2009) and Fauvet et al. (2011) where the Galactic magnetic field itself is modelled. We simply attempt to construct a plausible model representative of the synchrotron spectral index for use in the testing of Hi data reduction pipelines. Previous models of the synchrotron spectral index and emission maps have used Gaussian realisations to provide additional spatial resolution; for instance Remazeilles et al. (2015) use Gaussian realisations to provide a higher resolution estimate of the 408 MHz map. Component separation methods, however, behave differently when attempting to clean Gaussian or non-Gaussian structure; Spinelli et al. (2022) show that a variety of different techniques all struggle to deal with non-Gaussian foregrounds as viewed through a non-Gaussian (Airy) beam, while the same techniques do a better job of approaching the Hi power level when only Gaussian foregrounds are considered. Therefore we aim to determine if the additional complexity of non-Gaussian, 56 arcmin resolution spatial structure will provide a more accurately challenging test-bed for such techniques.

In this work we create a new all-sky spectral index template and assess how high resolution, non-Gaussian structure impacts the ability of a blind component separation method, Principal Component Analysis (PCA), to clean an emission cube of diffuse Galactic synchrotron emission in an attempt to measure the Hi auto-correlation power spectrum. Numerous foreground cleaning methods, both blind and parametric, are available for use and each method has different advantages and disadvantages when faced with spatial and spectral structure. It is not the aim of this work to provide a complete review of all existing component separation methods, as of such we choose to select one mainstream (i.e. often applied to intensity mapping data) technique to test here and make our spectral index map publicly available for the community to test alongside all other component separation techniques. We use our high resolution spectral index map, alongside the Haslam data at 408 MHz to form diffuse synchrotron emission templates at various frequencies, which can then be compared to existing synchrotron emission models. In this work we compare four models for synchrotron emission: 1) the GSM model, 2) the Haslam data scaled using the 5 degree Miville-Deschenes et al. (2008) spectral index map, 3) the Haslam data scaled using a version of the 5 degree Miville-Deschenes et al. (2008) spectral index map, which has had higher resolution angular detail up to 56 arcmin added to it using Gaussian realisations and 4)
the Haslam data scaled using the 5 degree Miville-Deschenes et al. (2008) spectral index map, which has had higher resolution angular detail up to 56 arcmin added to it using our trained CNN.

We choose to only consider one foreground emission: diffuse Galactic synchrotron emission, so as to assess the effect our different models have on cleaning. In reality, low level diffuse free-free emission, extragalactic points sources, residual radio frequency interference (RFI) and possibly even residual ground emission pick-up may also be present in observational data. The simulations include Hi emission, white noise and diffuse synchrotron emission and each frequency channel is convolved with a frequency-dependant beam. We use the online Hi  simulations repository FastBox1 to provide the test-bed set-up and so these simulations are specifically focused on foreground removal for a single-dish intensity mapping experiment using the MeerKAT dishes; an experiment such as MeerKLASS (Santos et al., 2016). However, the spectral index map that we have produced can be used by any experiment (single-dish or interferometric) alongside the Haslam data to get an estimate of diffuse synchrotron emission.

Footnote 1: [https://github.com/philbull/FastBox](https://github.com/philbull/FastBox)

The paper is laid out as follows: Sect. 2 is split between the description of the construction of our high resolution spectral index map using CNNs and the description of our simulation test set-up i.e. each of the simulated components and the chosen component separation method. We also assess the success of our CNN in building a realistic level of high resolution structure within our spectral index map in Sect. 2. In Sect. 3 we go on to use spherically-averaged auto-correlation power spectra to assess the impact of the different foreground maps on the cleaning ability of PCA. Our conclusions are presented in Sect. 4.

## 2 Method

### Constructing a high resolution spectral index map

Spectral index maps can be formed from data taken at two different frequencies (\(\nu_{A}\) and \(\nu_{B}\)); the map produced gives the average spectral index per pixel required to scale the emission from \(\nu_{A}\) to \(\nu_{B}\) (as long as said emission can be modelled as a power law). The spectral index is an average over frequency if and only if it is believed that the spectral index changes across frequency. It is important to select two sets of observational data which can be believed to only contain the emission of interest. For example, in this work we want to create a map of the synchrotron spectral index and so we need two data sets observed across either frequencies or regions of the sky we believe to be dominated by synchrotron emission. The Haslam 408 MHz data is the standard proxy for diffuse synchrotron emission. To create a spectral index map using these data we required another synchrotron dominated observational data set of the same, or higher resolution. Figure 1 shows both full and partial radio continuum maps, publicly available from the _WMAP_ legacy archive 2 as filled circles and two possible sources of future radio maps, from the MeerKLASS (Santos et al., 2016) and Bingo (Bartye et al., 2012) Hi intensity mapping experiments as dotted lines. As the MeerKLASS data span a range of resolutions they are denoted using a dotted rectangle. The only two maps available with higher resolution than the Haslam data (which can also be seen on Figure 1) are both at 1.4 GHz. We choose to use the CHIPASS radio continuum survey, which is available at 14 arcmin resolution.

Footnote 2: [https://lambda.gsfc.nasa.govl](https://lambda.gsfc.nasa.govl)

We only consider publicly available maps under 1.5 MHz, as at higher frequencies free-free emission is no longer a negligible component across the full sky. In Figure 2 we plot the typical decrease in synchrotron and free-free emission temperatures across frequency for both low and high Galactic latitudes. The _Planck_ full focal plane simulation data were used to provide the emission amplitudes 3; we used the mean synchrotron and free-free emission temperatures within a 5 degree squared region centred at Galactic latitudes 75\({}^{\circ}\) (to represent high latitudes) and 10\({}^{\circ}\) (to represent low latitudes). Both emissions are modelled as a power law with a synchrotron spectral index of -2.7 (Irfan et al., 2022) and a free-free spectral index of -2.1 (Bennett et al., 1992).

Footnote 3: [https://pla.esac.esa.int/maps](https://pla.esac.esa.int/maps)

#### 2.1.1 The training and testing data

To train our network we used high resolution (14 arcmin) data from the CHIPASS experiment (Calabretta et al., 2014) which cover the Southern sky at declinations < 25\({}^{\circ}\). Even if the data were full sky, it would be risky to try and construct a full sky synchrotron spectral index between 0.408 and 1400 MHz using the Haslam and CHIPASS data. At frequencies over 1 GHz free-free emission is no longer believed to be negligible across the majority of the sky; free-free emission can contribute up to 50 per cent of the total emission close to the Galactic plane (Platania et al., 1998). Therefore the full data were not used to train our network, instead we focused on the North Polar Spur (NPS) region which is believed to be a strong synchrotron emission feature (Haslam et al., 1964). This region is highlighted using red dotted lines in Figure 3 which shows the full CHIPASS data in Galactic coordinates. We restricted our training spectral index maps to come from this red dotted region only.

Figure 1: Publicly available, full and partial sky observational data under 1.5 GHz, as collated by the _WMAP_ legacy archive. The survey central frequencies and resolutions are shown and each survey is labelled by observatory name; some surveys are made up from several telescope observations and in those cases only one observatory name has been listed. Full survey details can be found in the relevant papers: DRAO 10 MHz (Caswell, 1976), DRAO 22 MHz (Roger et al., 1999), MU (Guzman et al., 2011), Parkes 85 MHz (Landecker and Wielebinski, 1970), Parkes 150 MHz (Landecker and Wielebinski, 1970), EDA2 (Kriele et al., 2022), Jodrell (Haslam et al., 1982), Dwingeloo (Berkhuijsen, 1972), Parkes 1.4 GHz (Calabretta et al., 2014), Stokert (Reich and Reich, 1986). Two future radio surveys, using the Bingo and MeerKLASS telescopes, are also plotted using dotted lines. The grey strip highlights surveys with resolutions higher than or equal to 56 arcmin.
The 1.4 GHz CHIPASS map was smoothed (assuming Gaussian beams) to 56 arcmin and used alongside the 56 arcmin, 408 MHz Haslam map to produce a spectral index map:

\[\beta=\ln\left(\frac{T_{\rm{\nu_{1}}}}{T_{\rm{\nu_{0}}}}\right)/\ln\left(\frac{ \nu_{1}}{\nu_{0}}\right), \tag{2}\]

where \(\nu_{0}=408\) MHz and \(\nu_{1}=1.4\) GHz. The destriped, reprocessed version of the Haslam map (Remazeilles et al., 2015) was used. As this map is available at HEALPix(Gorski et al., 2005) N\({}_{\rm side}\) 512, we downgraded the N\({}_{\rm side}\) 1024 CHIPASS map to N\({}_{\rm side}\) 512.

For the synchrotron spectral index map to have the correct mean both the Haslam map and the CHIPASS data must have the correct zero-levels. As the Haslam data are so widely used as a proxy for synchrotron emission, considerable work has already been done to determine the map zero-level. Wehus et al. (2017) used linear regression between multiple data sets to fit a zero-level of 8.9 K to the Haslam data. We adopt that value and subtract it from the Haslam data. To find the zero-level of the CHIPASS data we then used the same temperature-temperature linear regression technique used in Wehus et al. (2017). The linear regression between the CHIPASS and Haslam data within our selected NPS region is given as:

\[T_{\rm{1400}}(p)=m\times T_{\rm{408}}(p)+c, \tag{3}\]

where

\[c=m\times c_{\rm{408}}+c_{\rm{1400}}, \tag{4}\]

where \(T_{\rm{1400}}\) is the temperature per pixel within the NPS region in the 1400 MHz map, \(T_{\rm{1400}}\) is the equivalent at 408 MHz, \(m\) is the gradient fitted from the linear regression and the fitted offset (\(c\)) is a combination of the offsets in both maps. By taking the Haslam offset (\(c_{\rm{408}}\)) as 8.9 K, we could then calculate the CHIPASS offset. Figure 4 shows the linear regression within the NPS region between the Haslam and CHIPASS data; the Haslam data have already had the map zero-level of 8.9 K removed. We found a fitted zero-level of 3.21 K for the CHIPASS data; both the Haslam and CHIPASS maps were then used with their respective zero-levels subtracted to form a spectral index map. We show the spectral indices for our NPS region in Figure 5. In Sect. A we explore the uncertainty on our fitted zero-level of 3.21 K by performing linear regression across regions of different Galactic latitude and examine how this uncertainty propagates to uncertainties within our final spectral index map. We find a 1\(\sigma\) deviation of 0.1 K on the zero-level, which effects both the mean level of the spectral indices determined, by 3 per cent, as well as their spatial variations, by 11 per cent at the resolution of 56 arcmin.

The spectral index map produced was then smoothed to a resolution of 5 degrees, to provide both a low and high resolution perspective of the same data, and the data within our NPS region was cut up into numerous smaller sized spectral index maps. To increase the ease with which the network learnt the features across the numerous training spectral index maps, each separate map was normalised to contain spectral indices spanning from -1 to 1. Normalising the spectral index maps set the mean-level for all the maps to zero. This was not a problem however, as our interest is in the spatial variations around this mean; specifically the relationship between these variations at 56 arcmin and at 5 degree resolution. We selected 39 overlapping patches within this region, each spanning 7.3 degrees and having dimensions \(64\times 64\) pixels. A smaller number of patches of larger dimensions could have been chosen, but we opted for smaller images to reduce the number of features the network had to learn and to increase the number of training examples. Additionally, to maximise the number of training spectral index maps, we rotated each map three times by 90 degrees giving a total of \(39\times 4\) training maps. We formed these maps at both 56 arcmin and 5 degrees resulting in 156 pairs of maps.

Figure 4: Linear regression between the Haslam and CHIPASS data within the NPS region. The fitted offset is stated on the plot.

Figure 3: The CHIPASS 1.4 GHz map. The region used to create our synchrotron spectral index training data is highlighted in red.

Figure 2: Typical temperature variation across frequency for both synchrotron and free-free emission at high (\(b_{\rm{hi}}\)) and low (\(b_{\rm{ho}}\)) Galactic latitudes. The typical emission amplitude values were taken from the _Planck_ full focal plane simulations and power law models were assumed for both with a spectral index of -2.7 for synchrotron and -2.1 for free-free emission.
#### 2.1.2 The network

CNNs can be built up using a variety of architectures. Krachmalnicoff and Puglisi (2021) use generative adversarial neural networks (GANs) to learn the 12 arcmin features from pairs of 12 and 80 arcmin thermal dust images taken from the GNNLC all-sky thermal dust model (Planck Collaboration et al., 2016). Initially we began by training our CNN using a GAN architecture but found this structure failed to converge for our training data, possibly due to the small number of training information available as we could only use a small, synchrotron dominated fraction of the full sky. In the future, when experiments like MeerKLASS and Bingo publicly release large-area maps of the MHz sky, the GAN architecture can be revisited. For the current level of data availability however, we found the U-Net architecture optimal for our goals (Ronneberger et al., 2015) and made use of the Keras python library. U-Net CNNs use a 'U'-shaped symmetric structure of convolutional layers followed by deconvolution layers.

We have 156 pairs of synchrotron spectral index maps each of dimension \(64\times 64\) pixels. Table 1 shows the details of each layer in the U-Net network used in this work. The first layer has 8 filters expanding the original map size from \(64\times 64\) to \(64\times 64\times 8\), for the following layers the filters are doubled using a stride of 2. We picked a kernel size of 8, ensuring that the kernel size is divisible by the stride size to reduce the 'checkerboard' effect in the final images 4. Following the example of the generator 5 used in Krachmalnicoff and Puglisi (2021) we used the LeakyReLU activation function with a slope of 0.2, batch normalisation (to reinstate a zero mean and a variance of 1) after each convolution and a tanh activation for the final layer.

Footnote 4: [https://distill.pub/2016/deconv-checkerboard/](https://distill.pub/2016/deconv-checkerboard/)

Footnote 5: [https://github.com/ai4cmb/ForSE](https://github.com/ai4cmb/ForSE)

We held back 25 per cent of the spectral index map pairs for testing and used the other 75 per cent for training. To train our network we minimised the Mean Squared Error (MSE) loss function using the Adam optimiser for gradient descent with an initial learning rate of 0.0001. Use of the ReduceLRonPlateau option enabled the network to reduce the learning rate by a factor of 0.1 after seven iterations (epochs) with zero improvement in the loss function. The network was set to train in batches of size six and was allowed to stop whenever the loss function ceased to decrease after ten iterations. The testing spectral index maps were not used to train the network, but instead to evaluate it (see Sect. 2.1.4).

#### 2.1.3 Processing full-sky data

After the network had been trained the goal was to take the full-sky spectral index map of Miville-Deschenes et al. (2008) at 5 degree resolution and use it as the input for the network to generate a 56 arcmin resolution version. We obtained the Miville-Deschenes et al. (2008) map from the _Planck_ Full Focal Plane simulations (Planck Collaboration et al., 2016), using Equation 2 where \(T_{\mathrm{vp}_{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{ \mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{\mathrm{ \mathrm{    }}}}}}}}}}}}}}}}\) are the simulated synchrotron temperature maps at 353 and 217 GHz. The _Planck_ Full Focal Plane synchrotron simulations are available at N\({}_{\mathrm{side}}\) 2048, so we downgraded the synchrotron spectral index map to N\({}_{\mathrm{side}}\) 512. The \(12\times 512\times 512\) pixels were then projected into

\begin{table}
\begin{tabular}{c c c} \hline \hline
**Operation** & **Dimensions** & **Hyperparameters** \\ \hline \hline Input & \(64\times 64\times 1\) & \\ Convolution 2D & \(64\times 64\times 8\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Convolution 2D & \(32\times 32\times 16\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Convolution 2D & \(16\times 16\times 32\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Convolution 2D & \(8\times 8\times 64\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Convolution 2D & \(4\times 4\times 128\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Convolution 2D & \(2\times 2\times 256\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Up sampling 2D & & \\ Convolution 2D & \(16\times 16\times 32\) & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\ Leaky ReLU & & \\ Leaky ReLU & & \(\alpha=0.2\) \\ Batch normalisation & Momentum = 0.5 \\
768 spectral index maps of dimensions \(64\times 64\) pixels using healpy functions (Zonca et al., 2019).

Once the 768 generated spectral index maps were projected back onto the sphere, we enlisted the technique of Cycle-spinning (Coifman and Donoho, 1995) to remove any border effects caused as a result of the projections between the 3D sphere and the 2D image plane. Cycle-spinning involves rotating the full-sky, 5 degree spectral index map, breaking up the map into 2D maps which are then used by the network to estimate high resolution 2D maps, re-projecting these maps to form a high resolution all-sky map and then performing the inverse rotation. We perform 12 rotations in total, plus the original map: Axis (X,Y): \((0^{\circ},0^{\circ})\), \((45^{\circ},0^{\circ})\), \((0^{\circ},90^{\circ})\), \((-45^{\circ},0^{\circ})\), \((0^{\circ},-90^{\circ})\), \((45^{\circ},90^{\circ})\), \((45^{\circ},-90^{\circ})\), \((-45^{\circ},90^{\circ})\), \((-45^{\circ},90^{\circ})\), \((-45^{\circ},90^{\circ})\), \((-45^{\circ},-90^{\circ})\), \((90^{\circ},90^{\circ})\), \((90^{\circ},-90^{\circ})\), \((-90^{\circ},90^{\circ})\) and \((-90^{\circ},-90^{\circ})\). These thirteen maps were then averaged together. Lastly we smoothed the maps from the pixel resolution (6.87 arcmin for an N\({}_{\rm side}\) 512 map) to 56 arcmin.

#### 2.1.4 Assessing the synchrotron spectral index map

An example of a training spectral index map pair given to the network is shown in Figure 6. The top left panel shows a normalised \(64\times 64\) pixel spectral index map at 5 degree resolution, while the top right shows the same but for the 56 arcmin spectral index map. In the lower panel we have the map generated by the trained network having been given the 5 degree map as input. The network can be seen to have learnt the key features in the high resolution map. There are clearly pixel effects in the generated map, such as a faint checkerboard effect and the odd spurious (not true to the high resolution map) pixel. However, as we are not attempting to create a spectral index map at the pixel resolution (6.87 arcmin for an N\({}_{\rm side}\) 512 map) these spurious effects get removed by the Cycle spinning and smoothing processes described in Sect. 2.1.3.

To see how well the network can generate high resolution spectral index maps from a low resolution map not used in the training we use a map from the test subset. The top left panel of Figure 7 shows a low resolution map from the test subset, and the accompanying high resolution map in the top right panel. The network generated map is shown in the lower panel. It can be seen that the reproduction of the small scale structure is far less faithful to the true high resolution map than in the case of the training data shown in Figure 6. However, on visual inspection the level of spatial structure in the generated map seems appropriately detailed. The histogram distributions of map pixels can be used as a method to assess image complexity. In Figure 8 we show the histograms for the same maps displayed in Figure 7. The low resolution training map has one strong peak but other than that has a very flat histogram distribution, while both the 56 arcmin map and the network generated map show considerably more structure. The generated test maps are not correct, in that they are not identical to the true 56 arcmin test data, but they contain high resolution spatial structure and they remain faithful to the large scale structure in the map therefore they can be used to provide a high resolution spectral index map.

The original 5 degree resolution, all-sky spectral index map of Miville-Deschenes et al. (2008) is displayed in the top panel of Figure 9, whilst our 56 arcmin version is presented in the lower panel. Our 56 arcmin map is publicly available 6. Figure 10 selects

Figure 8: Histogram distribution of the normalised spectral indices in a test high and low resolution maps and the network generated high resolution model. These distributions are for the same data show in Figure 7

Figure 6: _Top left:_ One of the normalised 5 degree spectral index training maps, _top right:_ the equivalent 56 arcmin version of the map, _bottom:_ the equivalent map generated by the trained network.

Figure 7: _Top left:_ One of the normalised 5 degree spectral index test maps, _top right:_ the equivalent 56 arcmin version of the map, _bottom:_ the equivalent map generated by the trained network.
a smaller region of the all-sky maps to clearly demonstrate the additional detail in the high resolution version of the Miville-Deschenes et al. (2008) spectral index map. The healpy anafast library was used to calculate the angular power spectra for both the original and our new spectral index map, shown in Figure 11. The angular power spectrum for the original map drops off at \(\ell\sim 40\) / \(\theta\sim 4.5^{\circ}\) while the CNN generated spectral map power only starts to drop off at \(\ell\sim 150\) / \(\theta\sim 1.2^{\circ}\). We can see that the high resolution spectral index map suffers from a slight power loss (between 5 and 15 percent) from the original map between \(\sim 6\) and 18 degrees, as our CNN fails to perfectly reconstruct high-resolution images different to those that it has been trained with. We also add the angular power spectrum of a synchrotron spectral index map generated by adding high resolution structure to the 5 degree map using a Gaussian realisation with Gaussian structure up to 56 arcmin. The high resolution Gaussian structure was added using a power-law in \(\ell\): \(A\times(36/\ell)^{\beta}\), where the amplitude was set by the 5 degree spectral index map power at \(\ell=36\) and \(\beta=2.4\) following the parametrisation for synchrotron emission detailed in Santos and Cooray (2006). This will be pertinent for the next section, where we discuss the full simulation and four options for synchrotron emission models.

### The simulation

Having constructed a high resolution spectral index map, we wished to investigate whether or not such a map would impact the results of a foreground cleaning simulation. We used the Hi simulation suite FastBox to set up two data cubes of dimensions \(128\times 128\times 128\). Both cubes cover a 3600 square degree region with a pixel resolution of \(0.47^{\circ}\) (\(60^{\circ}\) across 128 pixels). The first data cube spans the frequency range 1220 to 1363 MHz over 128 pixels, giving a frequency resolution of 1.1 MHz and the seconds cube spans from 1084 to 1226 MHz, also with a frequency resolution of 1.1 MHz. The total data cubes consist of Hi emission plus diffuse synchrotron emission and are convolved with the MeerKAT beam, instrumental white noise is then added into the simulations after the beam convolution.

The FastBox Hi signal is simulated as:

\[\Delta T_{b}(\vec{x},z)=\overline{T}_{b}(z)\,b_{\rm HI}(z)\,\delta_{m}(\vec{x },z), \tag{5}\]

details of the mean brightness temperature, HI bias and fractional HI density used can be found in Irfan and Bull (2021). Log-normal transformations of the Gaussian field are applied to ensure a physical density distribution and the effect of Redshift Space Distortions are included by shifting each 3D pixel of the transformed field.

For the foreground contribution we only simulate diffuse synchrotron emission, as we wish to explore different synchrotron models. All our foreground models cover the 3600 degree square region of \(25^{\circ}<\) Galactic latitude (\(b\)) \(<85^{\circ}\) and \(220^{\circ}<\) Galactic longitude (\(l\)) \(<280^{\circ}\). The first model is the Global Sky Model; the GSM is in fact a model of the total diffuse Galactic emission at the user-selected frequency. As the simulated data are at high Galactic latitude and MHz frequencies the assumption is that synchrotron emission will be the dominant emission, but this assumption may not necessarily hold true within the GSM map which will contain some fractional level of extrapolated free-free, anomalous microwave and thermal dust emission. Models two, three and four are all models of pure synchrotron emission assuming a power law model form and using the Haslam 56 arcmin data as the emission amplitude at 408 MHz. These three models only differ in the resolution of the spectral index map used. The 5 degree spectral index map of Miville-Deschenes et al. (2008) is the base for all three spectral index maps and provides all the spectral index information for model two. Model three uses the new spectral index map produced by this work, which has high resolution information between 56 arcmin and 5 degrees provided by our CNN and model four has the high resolution information

Figure 11: Angular power spectra of the 5 degree, CNN and Gaussian generated synchrotron spectral index map.

Figure 10: A \(20^{\circ}\) by \(20^{\circ}\) zoom-in of the all-sky spectral index map at (_left:_) 5 degree and (_right:_) 56 arcmin resolution.

Figure 9: Full sky synchrotron spectral index between 408 MHz and 23 GHz at 5 degrees (_top_) and 56 arcmin (_bottom_) resolution.
simulated using a Gaussian realisation. Figure 12 shows the three different spectral indices per pixel for our 3600 square degree region. Common large-scale structure is clear in all maps but the 5 degree map is clearly missing the high resolution structure seen in the other two. The Gaussian high resolution map looks to be simply more noisy than the CNN map, as opposed to containing high resolution structure. The four foreground emission models are summarised in Table 2.

By introducing high resolution spatial structure we aim to determine whether said structure poses a problem for component separation methods after the sky signal has been convolved with a complex beam structure, i.e. a frequency changing beam with sidelobes which cannot be modelled at each frequency as Gaussian. Through the FastBox set-up we make use of the L-Band, Stokes I katbeam model 7 which models the frequency changing beam as a cosine aperture taper (Mauch et al., 2020). A 1D slice through the 2D beam pattern as a function of frequency is show in Figure 13. The approximate Gaussian FWHM for the MeerKAT beam for the frequency ranges under investigation in this work is 1.54\({}^{\circ}\) to 1.13\({}^{\circ}\).

Footnote 7: [https://github.com/ska-sa/Katbeam](https://github.com/ska-sa/Katbeam)

Only the combination of the synchrotron emission plus the Hi signal is convolved with the beam; the final constituent of the total emission model: instrumental noise, is added after beam convolution. We assume Gaussian instrumental noise with a standard deviation calculated from the radiometer equation:

\[\sigma_{\rm rms}=\frac{T_{\rm sys}}{\sqrt{N_{d}\,t_{\rm res}\,\delta\nu}}, \tag{6}\]

where

\[T_{\rm sys}=T_{r}+60\left(\frac{\nu}{300}\right)^{-2.5}, \tag{7}\]

where \(T_{r}=16\) K (based on typical MeerKAT receiver temperatures (Wang et al., 2021)), \(N_{d}\) is the number of available dishes which we set to 64 to match the MeerKAT array, \(\delta\nu\) is the frequency resolution and \(t_{\rm res}\) is the observational time per pixel. Figure 14 shows the total emission (Hi plus synchrotron emission convolved with the beam and then added to the instrumental noise) at 1273 MHz for the four different synchrotron emission models. The three models which use a scaling of the Haslam data for the foregrounds are indistinguishable by eye, whereas the total emission model which uses the GSM to provide the synchrotron emission template is quite distinct. All the models share common identifiable features however, such as regions of compact emission.

In Figure 15 the spectral form for a single pixel in each of the four total emission models is shown. The pixel temperature is multiplied by frequency squared in order to highlight any deviations from a simple power law model in the spectral form. The GSM total emission cube can be seen to display a completely different spectral form to the other three models, which is unsurprising as it is the only model not formed using a power law parameterisation. The other three models were formed from power laws but convolution with the MeerKAT beam has resulted in the averaging together of neighbouring pixels, which destroys the simple power law spectral form over frequency. The more complex the beam, the more complex these spectral perturbations and in-turn the harder it becomes to remove foreground structure when using a component separation technique that replies on spectral smoothness.

\begin{table}
\begin{tabular}{c c c} \hline \hline
**Model name** & **Emission** & **Resolution** \\ \hline GSM & total diffuse Galactic & 56 arcmin \\
5\({}^{\circ}\) & diffuse Galactic synchrotron & 56 arcmin/ 5\({}^{\circ}\) \\ CNN & diffuse Galactic synchrotron & 56 arcmin/ 56 arcmin \\ Gaussian & diffuse Galactic synchrotron & 56 arcmin/ 56 arcmin \\ \hline \hline \end{tabular}
\end{table}
Table 2: The four foreground emission models investigated in this work and their resolutions. The last three models require both amplitude and spectral index templates which is why they have two resolutions stated for each; the first for the amplitude, the second for the spectral index.

Figure 12: The three spectral index maps investigated in this work: the 5\({}^{\circ}\) model is on top, the CNN model is in the middle and the model with high resolution information filled in using a Gaussian realisation is on the bottom row.

Figure 13: A 1D slice through the katbeam beam model as a function of frequency.
### The foreground clean

Numerous methods of component separation can be used to attempt to clean total emission maps from foregrounds. Our aim is to determine if an individual technique behaves differently for our four sets of total emission cubes. Therefore we select a single method of foreground cleaning and explore the successes and failures of that technique for our four simulation sets. The method of choice for this work is Principal Component Analysis, which is a blind technique meaning that it requires no parametric information about the foreground emissions. PCA relies on the assumption that the foreground emission is smooth over frequency, whilst the cosmological signal is mainly stochastic over frequency. The technique works by performing an eigen decomposition of the data frequency-frequency covariance matrix;

\[C(\nu,\nu^{\prime})=\frac{1}{N_{\rm pix}}\sum_{j}^{N_{\rm pix}}\delta T_{j}(\nu )\delta T_{j}(\nu^{\prime}). \tag{8}\]

The largest eigenmodes are then labelled as foregrounds and subtracted away from the total data. The number of these large eigenmodes to remove is the only input parameter required by the PCA algorithm from the user.

## 3 Results

To asses the impact of our four different total emission models on the success of a PCA clean to the full data cube we enlist the use of spherically averaged auto-correlation power spectra. In place of plotting the power for each cube, we plot the ratio between the cleaned cube power spectrum and the Hi plus white noise power spectrum. An ideal foreground removal technique would remove all the synchrotron emission leaving only the Hi signal (after convolution with the beam) plus the instrumental noise thus giving a power spectra ratio of 1 on our plots. However, our focus is not on which foreground model results in the ratio closest to one but rather, if there are significant differences between the different total emission models.

In Figure 16 we summarise the current understanding of synchrotron spectral index modelling for the problem of component separation. The auto-correlation power spectra shown are all for PCA-cleaned total emission cubes, the number of PCA modes removed from the total emission cubes is always 1. The synchrotron emission models for the four spectra were provided by the Haslam data scaled to different frequencies using four different spectral index maps. The first, and simplest spectral index map is simply a constant value at each pixel. We use the mean value of the 5 degree Miville-Deschenes et al. (2008) spectral index map, -2.93, as our constant spectral index value. The second map is a Gaussian distribution of values with a mean of -2.93 and the same standard deviation as the 5 degree spectral index map, 0.06, with a spatial resolution of 5 degrees. The third map is the actual 5 degree spectral index map and the fourth map is the third map smoothed to a resolution of 10 degrees. As previously identified in the literature, using a single, constant value for the synchrotron spectral index across all

Figure 14: The total emission models at 1273 MHz for _top to bottom_: the GSM model, the 5\({}^{\circ}\) model, the CNN model and the Gaussian model.

Figure 15: The spectral form of a single pixel in each of the four total emission models.
pixels presents a simplified version of the true problem to the component separation algorithm. The fact that the three 5 and 10 degree power spectra display the largest power excesses over the Hi  plus white noise power indicates that the foreground removal is particularly challenging for a spatially complex synchrotron spectral index. The Gaussian 5 degree spectral index map, however, can be seen to place an excess of power across different (smaller) angular scales to the true 5 degree spectral index map and therefore is not an accurate representation of the foreground cleaning problem at medium angular scales. Whereas, the power spectra ratios with the largest values at large scales (small \(k\) values), as opposed to medium scales, are the two where a non-Gaussian, spatially changing spectral index has been used. As there is a significant difference between the cleaned power spectra for the true 10 and 5 degree spectral index map this plot motivates the goal of this paper, which is to determine if a new spectral index map with non-Gaussian spatial structure and a resolution of 56 arcmin is required to further bolster Hi  component separation test-beds.

In the next plot we return to the four emission models set-up in Sect. 2.2: namely the GSM, the 5 degree spectral index model, the 56 arcmin CNN spectral index model and the 56 arcmin Gaussian spectral index model. Figure 17 presents the auto-correlation power spectra plots for the four total emission models across the two frequency ranges investigated. The top plot is for the 1220.0 to 1362.6 MHz frequency range whilst the bottom two are for the 1084.0 to 1225.6 MHz frequency range. At lower frequencies the foreground emission is brighter and so a larger number of modes may be needed to achieve a suitably effective clean, which is why we include both the 1 and 2 mode removal results for the 1084.0 to 1225.6 MHz range. The GSM total emission cube was only ever cleaned using 2 mode removal as 1 mode left visible foreground structure in the cleaned maps. It is also clear from the plots in Figure 17 that the GSM model results in a completely different cleaned result compared to the other three models which are formed from a power law extrapolation of the Haslam data. The GSM is a total diffuse emission model and therefore should be used as such, as opposed to as a proxy for purely synchrotron emission. It has a more complex spectral structure than a simple power law, which is an advantage to those who wish to test their component separation technique of choice on such a data set. However, for those who intend to add their own compact sources, free-free emission, residual RFI etc., using the GSM as a proxy for solely synchrotron emission may make for an unrealistically hard to clean total emission cube.

The three models which are based on a power law extrapolation of the Haslam data can be seen to perform similarly in cleaning from Figure 17. In fact for the 1084.0 to 1225.6 MHz frequency range when two foreground modes are removed and the total emission cube is over-cleaned, i.e. Hi signal itself is removed causing the power spectra ratio to drop below one, the three power law models are essentially identical. For the 1 mode clean, for both frequency ranges investigated, we can see a difference in the three power law models at the large scale, low \(k\) end of the power spectra. At the largest scale plotted the 5 degree model and the Gaussian model behave similarly, while the CNN model gives different results. For the first frequency range the cleaned maps are between 3 and 3.5 times higher than the ideal measured signal (Hi  plus white noise), the difference between 3 and 3.5 is noteworthy, although within the error budget. However for the second range, all three methods are more than 100 times higher in power than the ideal signal and

Figure 16: Spherically averaged auto-correlation power spectra. The power spectra ratio between cleaned total emission cubes and data cubes containing only Hi emission convolved with the beam plus instrumental noise are plotted. In all four synchrotron emission models the Haslam data have been scaled using a spectral index map. The four different plots are for four different synchrotron spectral index maps.

Figure 17: Spherically averaged auto-correlation power spectra. The power spectra ratio between cleaned total emission cubes and data cubes containing only Hi  emission convolved with the beam plus instrumental noise are plotted. The top plot is for total emission cubes covering the 1220.0 to 1362.6 MHz frequency while the bottom two plots are for the 1084.0 to 1225.6 MHz. The difference between the middle and the bottom plot is the number of foreground modes removed for the PCA clean.
so the difference between the three methods becomes insignificant compared to the overall inability to recover the Hi signal. The conclusion to be taken from this is that it is worthwhile to pursue a more spatially complex foreground model, such as the CNN model presented in this work, if the component separation method under testing is performing well. If on the other hand, a new method is being honed which is orders of magnitude away from recovering the Hi signal, there is no need to pursue a more complex test-bed set-up than using the Haslam map extrapolated to different frequencies with the Miville-Deschenes et al. (2008) spectral index map. Specifically, for the MeerKLASS experimental set-up and the foreground cleaning method of PCA there is no significant difference between any of our 56 arcmin and the 5 degree spectral index models i.e. spectral index spatial complexity beyond a 5 degree resolution is not required for this particular test set-up.

## 4 Conclusions

Prior to this work it was known that assuming a spatially constant synchrotron spectral index would be an unrealistic simplification of the component separation problem for 21 cm intensity mapping experiments with frequency-dependant beams. It had also been shown that Gaussian foreground realisations are easier to separate from the Hi signal than non-Gaussian foregrounds when using blind foreground removal techniques. Motivated by the knowledge that increasingly spatially complex synchrotron spectral index models have been required to test the limitations of component separation techniques we created a 56 arcmin spectral index model to act as a proxy until future observational data allows us to measure the all-sky synchrotron spectral index at high angular resolution.

We have presented a model of the diffuse Galactic synchrotron emission spectral index between 0.408 and 23 GHz created from the 5 degree map of Miville-Deschenes et al. (2008) using a CNN trained on low and high resolution spectral index maps from CHIPASS 1.4 GHz observational data. Our map contains spatial structure up to 56 arcmin resolution. The intent is for this map to be used alongside the 56 arcmin 408 MHz map to create more realistic models for diffuse synchrotron emission across frequency for use as part of any simulation suite designed to test component separation techniques. It can be seen from the test subset of the data used to train the CNN that the small-scale structure generated is not an accurate representation of the true (empirically measured) small-scale structure, this coupled with the mid-range resolution power loss seen in the angular power spectrum means that our spectral index map cannot be used for science analysis of diffuse Galactic synchrotron emission. However, we believe that the CNN model presented in this work offers a useful contribution to simulation test-beds designed to probe the advantages and disadvantages of different component separation techniques.

Deep convolutional generative adversarial networks (DC-GANs) have been implemented for the simulation of thermal dust emission maps (Krachmalnicoff and Puglisi, 2021; Aylor et al., 2021) but for our input resolution (5 degrees) and relatively small number of training maps (156) we found it problematic to get a GAN network to converge and so found instead, the U-NET CNN architecture to be the optimum set-up for our purposes. The availability of high resolution CHIPASS data enabled the calculation of spectral indices between 408 and 1400 MHz at 56 arcmin which could then be smoothed so our network could be trained using pairs of 56 arcmin and 5 degree spectral index maps. Ideally we would have used publicly available high resolution MHz data, as at 1.4 GHz free-free emission is no longer negligible in certain regions of the sky. To combat the lack of available MHz data, we used the CHIPASS data but only within the North Polar Spur region; an area known to be dominated by diffuse synchrotron emission. In the future, however, we will be able to redo this analysis using MeerKAT or Bingo data to train our CNN.

To emphasize the use of a spatially complex foreground emission model, we set up four different total emission models and tested the ability of PCA to clean away the foregrounds, leaving the Hi plus instrumental noise. Our four emission models each contained the same Hi signal, the same level of instrumental Gaussian noise and the foreground and sky were convolved by the frequency-dependant katbeam beam model for each total emission model. The only difference between the four models were the foreground contributions. Emission model one had the synchrotron emission provided by the Global Sky Model, emission models two to four were all power law models of synchrotron emission made by extrapolating the Haslam map over frequency using a spectral index map. For model two the spectral index map used was the 5 degree map of Miville-Deschenes et al. (2008), for model three our high resolution CNN spectral index map was used and for model four the high resolution spatial information for the 5 degree spectral index map was provided using a Gaussian realisation. None of these four models are perfect, the GSM model is for total emission and therefore has a far more complex spectral form than a simple power law and the other three models assume that only synchrotron emission is present in the Haslam data. Having a selection of possible models with different characteristics, however, is useful for testing component separation techniques as that allows for an investigation into how each technique responds to different contaminants. Using the GSM allows for the testing of a technique in the presence of a complex spectral structure caused by multiple emission sources, whilst using a power law model with a spatially varying spectral index simulates the complex spectral structure caused by the interaction between a foreground emission and the telescope beam pattern. For our experimental setup (cosine aperture taper beam with a FWHM between 1.1 and 1.6 degrees and foregrounds cleaned using PCA) and an upper resolution bound of 5 degrees we have shown that increasing the resolution of the spectral index spatial structure changes the cleaned power spectrum and that this change occurs over different \(k\) scales depending on whether the spatial structure is Gaussian or non-Gaussian. At resolutions equal to or greater than 5 degrees, however, there is no significant difference between using the Miville-Deschenes et al. (2008) 5 degree spectral index map, a Gaussian 56 arcmin spectral index map or a non-Gaussian 56 arcmin spectral index map. This case-study is of particular pertinence for the MeerKLASS (Santos et al., 2016) component separation effort, while the 56 arcmin spectral index map made publicly available by this work is relevant to any other 21 cm experiment testing component separation techniques on their own unique experimental set-up.

In this work we have chosen to use the Miville-Deschenes et al. (2008) all-sky map as our estimate for the synchrotron spectral index at all resolutions larger than 5 degrees. As of such our CNN spectral index map is completely tied to this per-pixel spectral index estimate. We believe this is an appropriate choice given the prevalence of this spectral index map within both the CMB and intensity mapping communities as it is used within both the Python Sky Model and Planck Sky Model. However, it must be noted that any limitations or inaccuracies associated with the Miville-Deschenes et al. (2008) 5 degree spectral index map are shared by the CNN 56 arcmin spectral index map presented in this paper.
An additional complexity for the synchrotron spectral index is the very likely possibility that it changes, not only across pixels but also over frequency. Our CNN spectral index map, like the 5 degree spectral index map, provides the average (across frequency) spectral index per pixel between 0.408 and 23 GHz. The per-pixel spectral index value can, however, be scaled across frequency using a curvature model from the literature (Kogut (2012) for example), if required. We leave it to the community to scale the spectral index map as desired and to extend the tests shown in this work to include any other component separation technique and other non-Gaussian beam models.

## Data Availability

The 56 arcmin map of the simulated synchrotron spectral index is publicly available here: [https://github.com/melisirfan/synchrotron_emission](https://github.com/melisirfan/synchrotron_emission) and the Jupyter notebooks used in this analysis have been added to the Fastbox repository.

## Acknowledgements

M.I acknowledges support from the South African Radio Astronomy Observatory and National Research Foundation (Grant No. 84156) and would like to thank Mosima Nagip for her valuable insight on CNN architecture and Phil Bull and Mario Santos for the useful discussions.
