Foreground Cleaning for HI Intensity Mapping using GNILC Generalized Needlet Internal Linear Combination

Foreground Cleaning for HI Intensity Mapping using GNILC Generalized Needlet Internal Linear Combination Topollo Naketsana Astrophysics, University of the Western Cape, UWC Supervisors: Prof Mario Santos & Dr Karin Fornazier ∗ Abstract Intensity mapping is a powerful observational tool used to probe the large-scale structure of the Universe. It does this by measuring the combined light from many faint sources, helping us understand how matter is distributed. The 21 cm emission line of neutral hydrogen (HI) provides a powerful probe of the underlying matter distribution through intensity mapping. However, this cosmological signal is buried under bright Galactic and extragalactic foregrounds which are several orders of magnitude stronger than the HI signal. In this study, we implement the Generalized Needlet Internal Linear Combination (GNILC) method for foreground cleaning in HI intensity mapping using simulated maps and later the MeerKAT single-dish data within the redshift range z≈0.4–1.4. Simulations were constructed to include the cosmological HI signal and dominant astrophysical contaminants, forming a multifrequency data cube analyzed in needlet (wavelet) space. The GNILC algorithm employs a localized covariance analysis and constrained principal component separation to recover the HI signal while minimizing residual foreground contamination. Our results show a high correlation (rp= 0.956) between the recovered and true HI signal, with a root mean square residual of 0.26K, demonstrating the method’s strong ability to separate the faint HI signal of interest from the foregrounds (excluding point sources for the first simulation). These findings validate GNILC as a promising foreground removal approach for future large-scale structure surveys such as MeerKLASS, with further improvements needed. Keywords: HI 21 cm, intensity mapping, GNILC, MeerKAT, foreground cleaning, largescale structure, cosmology 1 Introduction The large-scale structure of the Universe carries the fingerprints of how matter has evolved over cosmic time. One of the most effective ways to trace this structure is by studying the 21 cm emission line from neutral hydrogen (HI), which was the most abundant gas during the Epoch of Reionization (EoR). This line comes from a spin-flip transition in the hydrogen atom and allows us to map where matter is distributed across different redshifts. By observing this signal, we can learn about how galaxies formed, how dark matter and dark energy influence the Universe, and measure important features like the Baryon Acoustic Oscillations (BAO) that trace the expansion history of the cosmos. The MeerKAT radio telescope, and in particular its MeerKLASS survey, provides an opportunity to study the HI signal in the redshift range z≈0.4–1.4. In this project, we focus on using MeerKAT in single-dish mode, which is especially sensitive to the large-scale. A major challenge in our study is that the HI signal is extremely faint compared to other astrophysical sources. Foregrounds like Galactic synchrotron emission, free–free emission, and extragalactic point sources are several orders of magnitude brighter. This ∗Email: topollonak[email protected] 1 makes foreground cleaning an important step to take in cosmology. To address this problem, we use the Generalized Needlet Internal Linear Combination (GNILC) method. GNILC is a component separation technique that works in needlet (wavelet) space and uses local covariance analysis to separate the faint HI signal from the much brighter foregrounds. In this project, we apply GNILC to simulated MeerKAT data to test how well it can recover the HI signal. The main goal is to understand how accurately GNILC performs in idealized conditions and what improvements are needed when additional components, such as point sources, are included. 2 Simulations and Experimental Setup To implement the Generalized Needlet Internal Linear Combination (GNILC) method in separating the HI cosmological signal from astrophysical foregrounds, we created a set of simulated sky maps that reproduce the main emission components observed by the MeerKAT telescope in single-dish mode. We used most of the python packages for preprocessing and healpy for spherical map handling and visualization. Each component of the total observed sky emission is represented as: xi(p)=si(p)+ni(p) (1) where xi(p) is the total observed brightness temperature at frequency channel iand pixel p, si(p) is the true HI cosmological signal, and ni(p) represents all astrophysical foregrounds and instrumental noise combined. 2.1 Simulated Components The simulation includes three main astrophysical components: 1. HI 21 cm signal: A Gaussian random realization generated using an input power spectrum consistent with ΛCDM cosmology. 2. Galactic synchrotron emission: Modeled as a spatially correlated field with a spectral index βsynch ≈ −2.8, scaled across frequency. 3. Extragalactic point sources: Simulated as randomly distributed sources with flux density following a power-law S∝ν−α, where α= 0.7. These maps are combined into a multi-frequency data cube with dimensions: Nch ×Npix = (252,786432) corresponding to 252 frequency channels and a HEALPix resolution of nside = 256. This corresponds to an angular resolution of approximately 1.55◦, suitable for probing multipoles in the range: from 30 to 300. The observing band corresponds to the Ultra High Frequency (UHF) range, covering redshifts z= 0.4−1.4. 2 Figure 1: Simulated foreground components at nside = 256.Ma: cosmological HI 21 cm signal. Top left HI 21 cm signal, top right is the Galactic synchrotron and extragalactic point sources. Each map is shown in brightness temperature units (mK) and corresponds to a single frequency slice 0.6 GHz in the simulated data cube. 2.2 Data Preparation and Channel Compression Since the GNILC method scales computationally with the number of frequency channels, we reduced the dimensionality from Nchan = 252 to Nch = 20. By averaging adjacent channels within narrow redshift bins. This compression preserves the large-scale cosmological information while reducing covariance noise in the analysis. 3 Component Separation using GNILC The Generalized Needlet Internal Linear Combination (GNILC) method is a foreground cleaning technique developed to separate the faint cosmological HI 21 cm signal from much brighter astrophysical emissions. It extends the standard Internal Linear Combination (ILC) approach by operating in needlet space - a localized wavelet domain on the sphere, which allows both frequency and spatial adaptivity (Olivari, 2016). 3.1 GNILC METHOD We model the observed data vector at each pixel pand frequency channel ias: xi(p)=si(p)+ni(p),(2) where si(p) represents the true HI signal, and ni(p) contains all foreground components and instrumental noise. For all channels, this can be written in matrix form: x(p)=s(p)+n(p),(3) where x,s, and nare Nch ×1 column vectors. 3 The total covariance matrix of the observed data is then given by: Rx(p)=⟨x(p)xT(p)⟩=RHI(p)+Rn(p),(4) where RHI is the covariance of the cosmological HI component and Rnrepresents the combined covariance of all foregrounds and noise. 3.2 Needlet Decomposition The GNILC algorithm first decomposes the observed sky maps into multiple needlet scales to capture spatially varying foreground properties. This provides localization in both pixel and multipole (ℓ) space. At each scale j, the data are filtered using needlet functions hj(ℓ), producing needlet coefficients: xj(p) = X ℓm hj(ℓ)xℓm Yℓm(p),(5) where Yℓm(p) are the spherical harmonics and hj(ℓ) defines the band-pass filter centered around a given multipole range. 3.3 Signal Subspace Identification To identify the dimensionality of the HI signal subspace, GNILC computes the local covariance of the data Rx(p) and compares it with the prior HI covariance Rprior HI (p). The covariance is then written as: ˜ Rx= (Rprior HI )−1/2Rx(Rprior HI )−1/2.(6) An eigen-decomposition of ˜ Rxallows separation of the subspaces corresponding to the HI signal and the dominant foregrounds. If λiare the eigenvalues of ˜ Rx, the subspace with λi≈1 corresponds to the HI component, while larger eigenvalues correspond to foreground contamination, this is called constrained PCA method. 3.4 Internal Linear Combination (ILC) After determining the dimension of the HI subspace, GNILC reconstructs the HI component using an Internal Linear Combination: ˆs(p)=W(p)x(p),(7) where W(p) is the weight matrix that satisfies two conditions: It provides unit response to the HI signal: W(p)a= 1, where ais the frequency mixing vector of the HI signal. It minimizes the total variance of the residual foregrounds: minimize ⟨ˆs2(p)⟩. ß The analytical solution for the optimal weight vector can be calculated using the Lagrangian multiplier method: W(p) = aTR−1 x(p) aTR−1 x(p)a.(8) 3.5 GNILC Workflow Summary The complete GNILC pipeline applied in this work follows these steps: 1. Input multi-frequency maps (xi(p)) containing HI + foregrounds. 4 2. Decompose the maps into needlet scales using hj(ℓ) filters. 3. Compute local covariance matrices Rx(p) in needlet space. 4. Compare Rx(p) with the prior Rprior HI (p) to determine the HI subspace via eigen-analysis. 5. Apply the ILC weights W(p) to reconstruct the HI component ˆs(p). 6. Recombine all needlet scales to produce the final recovered HI map. This method adapts to local variations in the signal-to-noise ratio, allowing better separation in regions dominated by strong Galactic foregrounds. Its performance is later assessed through the correlation coefficients, power spectra, and residual analyses described in the next section. 4 Results The GNILC algorithm was applied to the simulated data cube described in the previous section to test its ability to recover the faint HI 21 cm signal from strong astrophysical foregrounds. The performance was evaluated by comparing the recovered HI map with the true HI simulation and analyzing their correlation, residuals, and angular power spectra. Reults consists of two parts, the first one is when we run GNILC using only Synchrotron as the foreground and HI map (excluding the bright extragalactic point sources) and the second part is when we include the point sources. 1. HI Signal Recovery from Galactic synchrotron Figure 2: Shows the comparison between the true HI map and the recovered HI map obtained after GNILC cleaning when only Galactic synchrotron emission is included as foreground. The method successfully reconstructed the HI signal. Interpretation: The left map is the reconstructed HI signal with masked off parts of the sky by GNILC - thats where the galactic synchrotron is a strongest. The the right, the angular power spectrum of the input map, the true signal and the recovered signal. When additional components such as extragalactic point sources (radio galaxies, quasars, and AGNs) are included, the quality of recovery decreases, particularly in regions of strong contamination. 5 Figure 3: Recovered HI maps in the presence of both Galactic synchrotron and extragalactic point-source contamination. Left: GNILC reconstruction HI. Right: corresponding Angular power spectrum. These results shows that GNILC was able to recover the 96% of the signal in the absence of Point sources and couldn’t recover the signal for the case of Point sources. 4.1 Validation Metrics The recovered and true HI maps were compared in both pixel and harmonic space. The pixelspace correlation coefficient rpis defined as: rp=Cov(Trec, Ttrue) σrec σtrue .(9) For the simulation without point sources, we obtained: rp= 0.956, which indicates a strong correlation between the recovered and true HI signal. To further evaluate performance at different angular scales, we computed the multipole correlation coefficient rℓbetween the recovered and true angular power spectra Cℓ. Figure 4: Correlation coefficient rℓas a function of multipole ℓbetween recovered and true HI maps, showing strong agreement across most scales. Finally, the Root Mean Square (RMS) residual between the recovered and true temperature maps was computed as: RMS(∆T) = 1 NsX i (Trec −Ttrue)2,(10) which gives: RMS(∆T) = 0.26 K. 6 This low RMS value confirms that GNILC effectively minimizes foreground contamination while preserving the cosmological HI signal. Overall, these results demonstrate that GNILC successfully recovers the HI signal with high accuracy in the absence of point sources, validating the approach for future intensity mapping analyses using MeerKAT and similar experiments. 5 Discussion The GNILC method performed well in recovering the simulated HI signal from dominant foregrounds. From the results, a strong pixel-space correlation of rp= 0.956 and a low residual RMS of 0.26 K were achieved, showing that most of the large-scale HI structure was successfully reconstructed. The method was especially effective when only Galactic synchrotron emission was included, as shown in Figure 2. However, when point sources were added, the quality of recovery decreased (Figure 3). These sources introduce small-scale non-Gaussian structures that are difficult for GNILC to model, leading to higher residuals in the reconstructed maps. This indicates that further preprocessing, such as source masking, may be needed for real observations. Overall, the results show that GNILC can efficiently separate diffuse foregrounds and recover the HI 21 cm signal under idealized conditions, providing a solid step toward applying this method to MeerKAT single-dish data. 6 Conclusions We successfully tested the GNILC algorithm for foreground cleaning in simulated HI intensity mapping data. The recovered maps showed a strong correlation with the true HI signal (rp= 0.956) and a low residual RMS of 0.26 K, confirming that GNILC can efficiently separate diffuse foregrounds such as Galactic synchrotron emission. The method performs best under idealized conditions but struggles when bright point sources are included, indicating the need for improved preprocessing and masking strategies. Identified key parameters and derived quantities used in this analysis. These values correspond to the MeerKAT UHF band and define the scales probed by the GNILC-recovered HI maps. Parameter Formula Value z z = 1420/νobs −1 0.4−1.4 θFWHM θ= 1.22 λ/D ∼1.55◦ ℓ ℓ ≈π/θ 30 −300 SNR SNR = σHI/σn3.6891 Table 1: Main observational and derived parameters for the simulated MeerKAT single-dish HI data used in this study. In summary, GNILC shows strong potential for use in future HI intensity mapping surveys, such as MeerKLASS, by effectively recovering the cosmological signal at large scales. Further development, including point-source mitigation and beam-frequency modeling, will improve its performance on real observational data. 7 Acknowledgements Department of Physics and Astronomy at the University of the Western Cape. Thanks to my supervisors, Prof. Mario Santos and Dr. Karin Fornazier, for their invaluable support and guidance throughout this research. Their expertise and mentorship were essential to the development of this work. Acknowledge the MeerKLASS collaboration for their research framework and data references that guided this study, as well as the broader MeerKAT science community for their ongoing efforts in advancing HI intensity mapping in South Africa. Finally, I thank my peers and colleagues within the department for their encouragement and discussions that helped shape the direction of this project. References [1] Olivari, L. C. (2016). Generalized Needlet Internal Linear Combination (GNILC): Foreground cleaning method for HI intensity mapping. [2] Olivari, L. C. (2018). Intensity Mapping: A new approach to probe the large-scale structure of the Universe. [3] Santos, M. G., Cluver, M., Hilton, M., Jarvis, M., Jozsa, G. I. G., Leeuw, L., Smirnov, O., Taylor, R., Abdalla, F., & Moodley, K. (2017). MeerKLASS: The MeerKAT Large Area Synoptic Survey. [4] Dai, X., & Ma, Y. (2025). Expanded GNILC for Multi-frequency CMB/HI Separation. [5] De Caro, B., Carucci, I. P., & Fornazier, K. (2025). eGNILC approach for HI Intensity Mapping. 8