Deep Learning for CMB Foreground Removal and Beam Deconvolution: A U-Net GAN Approach

TL;DR

This paper presents a GAN-based U-Net approach for accurate CMB map reconstruction, effectively removing foregrounds and instrumental effects using realistic simulations, achieving high fidelity even in challenging regions.

Contribution

The study introduces a novel GAN architecture tailored for CMB reconstruction that corrects for foregrounds and systematics, demonstrating superior performance over traditional methods.

Findings

Achieves less than 1% difference between input and recovered maps outside the Galactic region.

Effectively corrects for foreground contamination and instrumental systematics.

Maintains low reconstruction error within the Galactic plane for temperature and polarization maps.

Abstract

Extracting cosmological information from microwave sky observations requires accurate estimation of the underlying Cosmic Microwave Background (CMB) by removing foreground contamination, instrumental noise, and the effects of beam convolution. In this work, we develop a machine learning-based approach for CMB reconstruction using a generative adversarial network (GAN) architecture, where the generator is modeled as a U-Net-based convolutional neural network. To train the network, we generate realistic microwave sky maps by simulating Planck-like observations: scanning HEALPix-simulated skies with real Planck beam profile, actual scan patterns, and anisotropic noise consistent with Planck data. Our method achieves high-fidelity reconstruction, with the difference between the input and recovered maps being less than $1\%$ (approximately $2\mu\mathrm{K}$ for temperature and less than $0.5\mu\mathrm{K}$ for polarization) outside the Galactic region. Even within the Galactic plane, the reconstruction error stays below $2$-$3\%$ for temperature maps across most regions, and is even smaller for polarization, apart from a few isolated pixels.. Most importantly, we demonstrate, for the first time, that a GAN-based method can effectively correct for foreground contamination, the systematic effects of non-circular beams and the asymmetric Planck scan pattern for both T and E-mode skymaps. Our results demonstrate the effectiveness of our method for robust and accurate recovery of the CMB signal, even in the presence of strong astrophysical foregrounds and instrumental systematics.

Figures

Tables

Table 1. Table 1 : Planck detectors and beam models used in the simulations. The detector polarization angle ψ pol \psi_{\mathrm{pol}} , polarization leakage ϵ \epsilon , and polarization efficiency ρ = 1 − ϵ \rho=1-\epsilon are taken from the Planck RIMO. Temperature beams correspond to representative detectors at each frequency, while polarization simulations use the detectors listed below.

Frequency (GHz)	Detector (Temp)	Detector (Pol)	${\psi }_{\mathrm{pol}}$ (deg)	$ϵ$	$\rho =1-ϵ$
70	18M	18M	–	–	–
217	217-1	217-6b	-44.100	0.0236	0.9764
353	353-1	353-3b	-46.100	0.0413	0.9587
545	545-1	–	–	–	–

Table 2. Table 2 : Frequency channels with corresponding Gaussian beam FWHM (used in circular convolution) and noise standard deviation ( σ \sigma )

Frequency [GHz]	FWHM [arcmin]	Noise $\sigma$ [$\mu K$]
70	13.31	1.31
100	9.66	1.15
143	7.27	1.78
217	5.01	1.91
353	4.86	4.66
545	4.84	7.99

Table 3. Table 3 : U-Net generator architecture. C in C_{\text{in}} and C out C_{\text{out}} denote the input and output channels. Convolutions use 3 × 3 3\times 3 kernels with stride 1 ( s 1 s_{1} ) and padding 1 ( p 1 p_{1} ). Max pooling uses 2 × 2 2\times 2 kernels with stride 2 ( s 2 s_{2} ). Upsampling uses transposed convolutions with 4 × 4 4\times 4 kernels, stride 2 ( s 2 s_{2} ), and padding 1 ( p 1 p_{1} ). The input has 6 channels for temperature (3 for polarization) the output has 1 channel. Repeated channel sizes reflect skip connections. (DS = Downsampling block; s1 = stride 1; s2 = stride 2; p1 = padding 1)

Layer	$C_{in}$	$C_{out}$	Operation
DS Block 1	6	4	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=4) + BatchNorm + LeakyReLU(0.2)
DS Block 2	4	4	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=4) + BatchNorm + LeakyReLU(0.2)
DS Block 3	4	8	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=8) + BatchNorm + LeakyReLU(0.2)
DS Block 4	8	8	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=8) + BatchNorm + LeakyReLU(0.2)
DS Block 5	8	8	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=8) + BatchNorm + LeakyReLU(0.2)
DS Block 6	8	8	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=8) + BatchNorm + LeakyReLU(0.2)
DS Block 7	8	16	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=16) + BatchNorm + LeakyReLU(0.2)
DS Block 8	16	16	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=16) + BatchNorm + LeakyReLU(0.2)
DS Block 9	16	16	Conv2D(3$\times$3, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(3$\times$3, s1, p1, out=16) + BatchNorm + LeakyReLU(0.2)
US Block 1	32	16	TransConv2D(4$\times$4, s2, p1) + Dropout(0.2) + LeakyReLU(0.2)
US Block 2	32	16	TransConv2D(4$\times$4, s2, p1) + Dropout(0.2) + LeakyReLU(0.2)
US Block 3	24	8	TransConv2D(4$\times$4, s2, p1) + Dropout(0.2) + LeakyReLU(0.2)
US Block 4	16	8	TransConv2D(4$\times$4, s2, p1) + BatchNorm + LeakyReLU(0.2)
US Block 5	16	8	TransConv2D(4$\times$4, s2, p1) + BatchNorm + LeakyReLU(0.2)
US Block 6	16	8	TransConv2D(4$\times$4, s2, p1) + BatchNorm + LeakyReLU(0.2)
US Block 7	12	4	TransConv2D(4$\times$4, s2, p1) + BatchNorm + LeakyReLU(0.2)
US Block 8	8	4	TransConv2D(4$\times$4, s2, p1) + BatchNorm + LeakyReLU(0.2)
US Block 9	8	4	TransConv2D(4$\times$4, s2, p1) + BatchNorm + LeakyReLU(0.2)
Final Layer	4	1	TransConv2D(4$\times$4, s2, p1)

Table 4. Table 4 : Layer Specifications for the Discriminator Architecture

Layer	$C_{in}$	$C_{out}$	Operation
DS Block1	$C_{in}+C_{out}$	4	Conv2D(4$\times$4, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(4$\times$4, s1, p1, out=4) + LeakyReLU(0.2)
DS Block2	4	8	Conv2D(4$\times$4, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(4$\times$4, s1, p1, out=8) + BatchNorm + LeakyReLU(0.2)
DS Block3	8	16	Conv2D(4$\times$4, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D (4$\times$4, s1, p1, out=16) + BatchNorm + LeakyReLU(0.2)
DS Block4	16	32	Conv2D(4$\times$4, s1, p1, out=$n_p$) + MaxPool2D(2$\times$2, s2)
			Conv2D(4$\times$4, s1, p1, out=32) + BatchNorm + LeakyReLU(0.2)
ZeroPad1	32	32	Zero Padding(1$\times$1)
ConvLayer1	32	16	Conv2D(7$\times$7, s1, no padding) + BatchNorm + LeakyReLU(0.2)
ZeroPad2	16	16	Zero Padding(1$\times$1)
Final Layer	16	1	Conv2D(3$\times$3, s1, no padding)

Table 5. Table 5 : Reconstruction results for patches 0-5. The columns show: Input map with foregrounds, the network’s predicted CMB, the true CMB, and the difference between the predicted and true CMB.

Input	Predicted	True	Diff

Table 6. Table 6 : Reconstruction results for patches 6-11. The columns show: Input map with foregrounds, the network’s predicted CMB, the true CMB, and the difference between the predicted and true CMB.

Input	Predicted	True	Diff

Equations

$$\hat{R}(t)=R_z({\omega }_{3}t)R_x({\omega }_{2}t)R_y({\theta }_1)\left[\begin{array}{c}\cos {\theta }_2\\ \sin {\theta }_2\cos ({\omega }_{1}t)\\ \sin {\theta }_2\sin ({\omega }_{1}t)\end{array}\right],$$

$$R_z({\omega }_{3}t)=\left[\begin{array}{ccc}\cos {\omega }_{3}t& \sin {\omega }_{3}t& 0\\ -\sin {\omega }_{3}t& \cos {\omega }_{3}t& 0\\ 0& 0& 1\end{array}\right],R_x({\omega }_{2}t)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\omega }_{2}t& \sin {\omega }_{2}t\\ 0& -\sin {\omega }_{2}t& \cos {\omega }_{2}t\end{array}\right],$$

$$R_y({\theta }_2)=\left[\begin{array}{ccc}\cos {\theta }_1& 0& \sin {\theta }_1\\ 0& 1& 0\\ -\sin {\theta }_1& 0& \cos {\theta }_1\end{array}\right].$$

$$\hat{S}(t)=R_z({\omega }_{3}t)R_x({\omega }_{2}t)R_y({\theta }_1)\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right].$$

$$A_{\text{pix}}=\frac{4\pi }{N_{\text{pix}}}.$$

$$\text{Scan time}=\frac{1}{{\omega }_1}\sqrt{\frac{4\pi }{N_{\text{pix}}}}.$$

$$T_{\mathrm{obs}}(\hat{n})=\int B(\hat{n},{\hat{n}}^{\prime })T({\hat{n}}^{\prime })d{\Omega }_{{\hat{n}}^{\prime }}+T_n(\hat{n})$$

$$T_{\mathrm{obs}}^s(\hat{n})=\frac{{\sum }_{i\in \mathcal{S}}{B}^s(\hat{n},{\hat{n}}_i^{\prime })T({\hat{n}}_i^{\prime })}{{\sum }_{i\in \mathcal{S}}{B}^s(\hat{n},{\hat{n}}_i^{\prime })},$$

$$T_{\mathrm{obs}}({\hat{n}}_j)=\sum _{i\in \mathcal{S}}\left[\frac{1}{h_j}\sum _s\frac{B^s(\hat{n},{\hat{n}}_i^{\prime })}{{\sum }_{i\in \mathcal{S}}{B}^s(\hat{n},{\hat{n}}_i^{\prime })}\right]T({\hat{n}}_i^{\prime })=\sum _{i\in \mathcal{S}}\tilde{B}({\hat{n}}_j,{\hat{n}}_i^{\prime })T({\hat{n}}_i^{\prime }),$$

$$\hat{N}=\cos \phi \hat{I}+\sin \phi \hat{Z}.$$

$$x_i={\theta }_i\cos {\alpha }_i,y_i={\theta }_i\sin {\alpha }_i.$$

$$d(\hat{n})=\left[I(\hat{n}),\rho Q(\hat{n})\cos (2{\psi }_t),\rho U(\hat{n})\sin (2{\psi }_t)\right],$$

$$\vec{P}=\cos ({\psi }_{\mathrm{pol}})\vec{I}+\sin ({\psi }_{\mathrm{pol}})\vec{Z},$$

$${\hat{e}}_{\varphi }=\frac{\hat{z}\times \hat{R}}{|\hat{z}\times \hat{R}|},{\hat{e}}_{\theta }={\hat{e}}_{\varphi }\times \hat{R},$$

$${\psi }_t={\tan }^{-1}\left(\frac{\vec{P}\cdot {\hat{e}}_{\varphi }}{\vec{P}\cdot {\hat{e}}_{\theta }}\right),$$

$$B_I(\hat{n},{\hat{n}}_i^{\prime })$$

$$B_Q(\hat{n},{\hat{n}}_i^{\prime })$$

$$B_U(\hat{n},{\hat{n}}_i^{\prime })$$

$$I_{\mathrm{obs}}(\hat{n})$$

$$Q_{\mathrm{obs}}(\hat{n})$$

$$U_{\mathrm{obs}}(\hat{n})$$

$${\mathcal{L}}_{L2}=\frac{1}{N_{\mathrm{pix}}}\sum _{i=1}^{N_{\mathrm{pix}}}(y_i-G(x{)}_i{)}^2,$$

$${\mathcal{L}}_{L1}=\frac{1}{N_{\mathrm{pix}}}\sum _{i=1}^{N_{\mathrm{pix}}}|y_i-G(x{)}_i|.$$

$${\mathcal{L}}_{\mathrm{adv}}=-\frac{1}{N_{\mathrm{disc}}}\sum _{j=1}^{N_{\mathrm{disc}}}\log D(x,G(x){)}_j.$$

$$\mathcal{L}={\lambda }_2{\mathcal{L}}_{L2}+{\lambda }_1{\mathcal{L}}_{L1}+{\lambda }_{\mathrm{adv}}{\mathcal{L}}_{\mathrm{adv}},$$

$${\mathcal{L}}_{\text{Disc}}=-\frac{1}{N_{disc}}\sum _{i=1}^{N_{disc}}\left[\log D(x,y{)}_i+\log \left(1-D(x,G(x){)}_i\right)\right]$$

$$m_t$$

$${\hat{m}}_t$$

$${\beta }_1=0.9,{\beta }_2=0.999,ϵ=10^{-8}.$$

$$T(\hat{n})=\sum _{l=0}^{\infty }\sum _{m=-l}^l{a}_{lm}{Y}_{lm}(\hat{n}),$$

$$⟨a_{lm}{a}_{l^{\prime }{m}^{\prime }}^{\ast }⟩=C_l{\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},$$

$$⟨a_{lm}{a}_{l^{\prime }{m}^{\prime }}^{\ast }⟩=\sum _{LM}{A}_{ll}^{LM}{C}_{lm{l}^{\prime }{m}^{\prime }}^{LM}{C}_{N{r}^{\prime }0}^{L0}\left({\Pi }_l{\Pi }_{l^{\prime }}/{\Pi }_L\right),$$

Full text

1 Introduction

Among the various tools that we have at our disposal for understanding the early evolution and current composition of our universe, probably none has been as successful as observations of the CMB [1, 2]. This faint afterglow of the hot early phase of our universe allows us to probe the universe roughly 380,000 years after the Big Bang, offering invaluable insights into its composition, evolution, and the seeds of cosmic structure. Precise measurements of CMB temperature and polarization anisotropies enable us to constrain cosmological parameters, search for signatures of primordial gravitational waves, and test key cosmological paradigms, including inflationary models and dark matter properties. However, the observed CMB signal is heavily contaminated by astrophysical foreground emissions, including synchrotron, free-free, and thermal dust radiation, as well as instrumental effects such as the telescope’s non-ideal beam response, complex scanning strategies, instrumental noise, etc [2]. These contaminants, spanning a wide range of frequencies and spatial scales, introduce significant challenges in accurately reconstructing the intrinsic CMB signal.

The spectral characteristics of galactic foreground emissions differ markedly from those of the intrinsic CMB signal. This distinction serves as the foundation for separating foreground contamination from the CMB. Traditional component separation methods, such as Internal Linear Combination (ILC) [3, 4, 5], its extension Needlet ILC [6], and various Bayesian approaches [7, 8], have been widely employed for foreground cleaning. ILC methods minimize variance in multi-frequency maps while preserving the CMB signal. While effective to some extent, these techniques often struggle with the complex, non-Gaussian nature of foregrounds, which violates assumptions made by many statistical techniques, and require detailed prior knowledge of their statistical properties. Furthermore, the combined effects of the telescope’s beam and its complex scanning strategy introduce significant spatial distortions and noise correlations, further complicating CMB map reconstruction. Correcting for these effects using conventional deconvolution techniques can be computationally expensive and prone to systematics [9].

In recent years, machine learning (ML) has emerged as a powerful tool for astrophysical data analysis [10, 11], particularly in problems involving high-dimensional and non-linear features. Deep learning architectures, such as convolutional neural networks (CNNs) and U-Nets, have demonstrated remarkable success in image denoising and pattern recognition tasks. A graph-based Bayesian convolutional neural network based on the U-Net architecture can predict cleaned CMB with pixel-wise uncertainty estimates. The U-Net, with its unique encoder-decoder structure and skip connections, excels at capturing both local and global features, making it ideally suited for the complex task of disentangling foregrounds from the CMB [12, 13, 14, 15]. By leveraging this capability, U-Nets can effectively learn the complex mapping between contaminated and clean CMB maps, simultaneously mitigating foreground contamination and correcting for beam distortions.

This work presents a Python-based pipeline where we explore generative adversarial network (GAN)-based approach built upon a U-Net architecture for the simultaneous removal of foreground contamination and deconvolution of beam effects in CMB temperature and polarization maps. We begin by generating realistic CMB simulations incorporating the non-circular beam and complex scanning patterns of the Planck satellite. These simulations are further contaminated with realistic foreground emissions and instrumental noise to create synthetic observations that closely resemble real CMB data. The U-Net model is then trained to reconstruct the clean CMB signal from these contaminated maps. Preliminary results indicate that our approach successfully suppresses foregrounds and mitigates beam distortions while preserving the crucial statistical properties of the CMB, and retrieving the original CMB map signal with high accuracy.

To simulate realistic foregrounds, we utilize the Python Sky Model (PySM), a widely used tool for generating full-sky simulations of Galactic emissions in intensity and polarization [16]. PySM allows us to model various components, including thermal dust, synchrotron, and free-free emission, using publicly available data and phenomenological models. By incorporating PySM, our simulations capture the complexity of the foregrounds, providing a more realistic test for our U-Net-based cleaning method.

While several CMB foreground cleaning methods exist, each comes with its own limitations and systematic biases. A comparative analysis of results from multiple independent techniques, including our U-Net approach, is crucial for validating inferences and quantifying systematic uncertainties, ultimately strengthening the reliability of cosmological inferences. Our analysis demonstrates that the GAN-based model not only suppresses foreground contamination and instrumental noise, but also performs effective beam deconvolution. In particular, it corrects for the statistical isotropy (SI) violation caused by noncircular beams and the asymmetric scan pattern—systematic effects that conventional foreground-cleaning methods are fundamentally unable to remove—thereby yielding a more accurate recovery of the true CMB sky. At present we have explored our analysis only on the CMB temperature and E-mode polarization skymaps. In future, this method can be extended to B-mode polarization studies and applied to other cosmological datasets, such as HI intensity mapping etc. The computational efficiency of deep learning approaches also makes them attractive for large-scale cosmological surveys, where traditional methods may be prohibitively expensive.

This paper has been organized as follows. In Section 2, we describe the simulation pipeline used to generate the realizations of the observed sky maps. These realizations are later used to train our network. Section 3 provides details of the neural network that we use for the foreground removal. The results from our analysis are presented in Section 4. We highlight the potential of this approach to advance CMB data analysis and also discuss challenges and limitations in Section 5. The final section, Section 6, presents our conclusions.

2 Simulation Pipeline for Generating Observed CMB Sky Maps

Training neural networks demands large and diverse datasets. However, only a single realization of the CMB sky is available to us, making direct data-driven training infeasible. To address this, we develop a robust simulation pipeline. The objective is to generate synthetic CMB sky maps that closely resemble real observational data, including astrophysical foregrounds, instrumental noise, and beam effects. This provides us with a realistic and controlled environment for training the neural network and rigorously testing the performance of its cleaning and deconvolution methods.

2.1 Sky Model and Foregrounds

Our simulation pipeline for CMB map generation proceeds in several interconnected stages designed to replicate the complexities of real observational data. We begin by using CAMB [17, 18]111https://camb.info/ to compute the theoretical angular power spectrum of the CMB, from which we generate multiple realizations of the CMB sky using HEALPix [19]222https://healpix.sourceforge.io/. In addition to the temperature anisotropy spectrum $C_{\ell}^{TT}$, the polarization spectra $C_{\ell}^{EE}$ is also used to generate the corresponding Stokes $Q$ and $U$ maps. In order to simplify the validation of the polarization response in our simulations, we set the $C_{\ell}^{BB}$ and $C_{\ell}^{EB}$ spectrum to zero.

From Stokes $Q$ and $U$ maps, we construct the scalar $E$-mode polarization maps by performing a spin-2 spherical harmonic decomposition of the polarization field. This transformation is implemented using the map2alm routine from the HEALPix package, which takes as input the Stokes maps $\{I,Q,U\}$ and computes the corresponding harmonic coefficients $(a_{\ell m}^{T},a_{\ell m}^{E},a_{\ell m}^{B})$.

The $E$-mode maps are then obtained by transforming $a_{\ell m}^{E}$ back to real space using alm2map. In our simulations, the $B$-mode power spectrum is set to zero ($C_{\ell}^{BB}=0$), so only the $E$-mode component contributes to the polarization signal. Same procedure applied to foreground maps as well.

To mimic foreground contamination present in actual CMB observations, we simulate astrophysical foreground emissions using the Python Sky Model (PySM) [20, 21, 16] 333https://github.com/galsci/pysm. PySM generates realistic, frequency-dependent emission maps across six frequency bands – 70, 100, 143, 217, 353, and 545 GHz – capturing the key foreground components, namely thermal dust emission, synchrotron radiation, and free-free (bremsstrahlung) emission. These foregrounds vary both spatially and spectrally, introducing non-trivial structure that complicates CMB signal recovery. For each frequency, multiple sets of foreground maps are generated in the Stokes I, Q, and U components to account for both intensity and polarization of the emission. We combine the synthetic CMB signal with these foregrounds, and create multi-frequency observations that closely resemble the real CMB sky.

2.2 Planck Scan Pattern

The Planck satellite scans the sky using a single beam that directly measures temperature, unlike the WMAP satellite, which performed differential measurements. The beam is positioned approximately $85^{\circ}$ off the satellite’s symmetry axis, with a precession angle of around $7.5^{\circ}$. The satellite spins at a rate of $180^{\circ}$/min, and its precession traces a complete revolution over six months. The scan pattern for the Planck scan [22, 23, 24] is shown in Figure 1.

The key angular parameters and rotational velocities involved in the scan geometry (see Figure 1) are

• $\theta_{1}=7.5^{\circ}$ (precession angle)
• $\theta_{2}=85^{\circ}$ (beam tilt from the spin axis)
• $\omega_{1}=2\pi$ rad/min (spin rate)
• $\omega_{2}=2\omega_{3}=2.3908\times 10^{-5}$ rad/min (precession rate of spin axis)
• $\omega_{3}=1.1954\times 10^{-5}$ rad/min (annual orbit rate)

Due to this scanning configuration, pixels near the ecliptic poles are visited more frequently than those near the equator, resulting in non-uniform sky coverage [25, 23, 24].

Let $\hat{R}(t)$ denote the line-of-sight (LOS) direction of the beam in the ecliptic coordinate system. It is expressed as a time-dependent composite rotation:

[TABLE]

Here, $R_{y}(\theta_{1})$, $R_{x}(\omega_{2}t)$, and $R_{z}(\omega_{3}t)$ are rotation matrices representing precession, spin, and orbital motion, respectively. Note that the subscripts $x$, $y$, and $z$ refer to the axes of the local coordinate system. In the context of the $3\times 3$ rotation matrix, these subscripts indicate the axis about which the rotation is performed, i.e. they are given by

[TABLE]

The satellite’s symmetry axis, $\hat{S}(t)$, is defined by rotating a reference vector through the same transformation without the final spherical component:

[TABLE]

To estimate the scanning duration per pixel, we consider the full sky solid angle $4\pi$ steradians divided among $N_{\text{pix}}$ pixels. Each pixel subtends:

[TABLE]

Assuming uniform beam motion, the time spent on a pixel is approximately:

[TABLE]

We compute $\hat{R}(t)$ at intervals of this scan time, and then convert each $\hat{R}(t)$ vector into a corresponding HEALPix pixel index. By accumulating the number of scans per pixel over a full year, we generate a simulated hit-count map representing the Planck scanning strategy.

As discussed, this scanning geometry leads to more frequent observations near the ecliptic poles. Figure 2 shows the resulting hit-count map on a logarithmic scale from a 12-month Planck scan, displayed in galactic coordinates. This hit-count map is later used to compute the spatially varying variance of the instrumental noise.

2.3 Beam Convolution and Polarized Detector Response

2.3.1 Temperature Convolution

The observed skymap is a convolution of the intrinsic sky with the beam function.

[TABLE]

Here $T(\hat{n}^{\prime})$ is the intrinsic sky temperature. $B(\hat{n},\hat{n}^{\prime})$ denotes the beam sensitivity between the central direction $\hat{n}$ and another direction $\hat{n}^{\prime}$. $T_{n}(\hat{n})$ is the instrumental noise.

To simulate the observational response from the Planck satellite, we convolved the HEALPix maps with the real Planck beams and scan pattern. Our analysis accounts for both the asymmetric shape of the instrument beam and its orientation at each sky location.

At each pointing direction $\hat{n}$, the observed temperature $T_{\mathrm{obs}}(\hat{n})$ is computed as a weighted average over neighboring pixel temperatures, using a direction-dependent beam response function as the weighting kernel:

[TABLE]

where $\hat{n}$ is the pointing direction, $\hat{n}^{\prime}_{i}$ is the central direction of the $i$th pixel and $\mathcal{S}$ is the set of all the pixels for which $B(\hat{n},\hat{n}^{\prime}_{i})$ is non-zero. In the previous section we discussed that during the Planck scan each of the pixels get scanned multiple number of times. Here the superscript $T^{s}_{\mathrm{obs}}(\hat{n})$ indicates the observed temperature of a pixel along $\hat{n}$ direction during $s$th scan (note that $\hat{n}$ can be any direction within the pixel and does not necessarily correspond to the direction of the pixel center). As each pixel is scanned multiple times, the final observed sky temperature for a pixel centered at direction $\hat{n}_{j}$ can be obtained by averaging the temperature values measured along all the directions that fall within that pixel across different scans, i.e. $T_{\mathrm{obs}}(\hat{n}_{j})=\frac{1}{h_{j}}\sum_{s}T^{s}_{\mathrm{obs}}(\hat{n})$. Here $h_{j}$ is the hitcount of the pixel centered along $\hat{n}_{j}$, and $\hat{n}$ represents all the directions that fall on that pixel. We can rearrange the above equation to get

[TABLE]

Here $\widetilde{B}(\hat{n}_{j},\hat{n}^{\prime}_{i})$ is an effective beam function. Since different pixels are scanned a varying number of times and from different orientations, the effective beam function differs from pixel to pixel. Generating thousands of independently scanned maps for training the network would be computationally expensive and time-consuming. To address this, we precompute and store the effective beam at each pixel. This effective beam is then multiplied with maps to obtain the scanned skymaps.

For generating the effective beam, we need to know the orientation of the beam at every time step. To define the beam orientation in the local sky patch, we construct an orthonormal coordinate system at each time step using the scan and beam directions (refer to Figure 1):

• $\hat{Z}=\hat{R}\times\hat{S}$ (perpendicular to both beam and spin axis)
• $\hat{I}=\hat{R}\times\hat{Z}$ (orthogonal to $\hat{R}$ in the beam scan plane)

The $\hat{I}-\hat{Z}$ plane is a plane perpendicular to the beam pointing direction and is the plane on which the beam is located. The scan axis of the non-circular beam, $\hat{N}$, is given by:

[TABLE]

where $\varphi$ is the angle of the scan axis of the beam relative to the $\hat{I}$ axis.

To evaluate $B(\hat{n}_{j},\hat{n}_{i}^{\prime})$, we project the angular displacement between the central pixel and its neighbors into a local 2D Cartesian coordinate system defined by the scan orientation at each time step. For a given observation time, the line-of-sight (LOS) direction of the beam is represented by the unit vector $\vec{R}(t)$, and the satellite’s spin axis is denoted by $\vec{S}(t)$. These vectors define a local orthonormal frame, which is used to orient the beam response function according to the current scan geometry. In our implementation, we fix the beam orientation angle $\varphi=0$, so that the effective scan-aligned beam direction $\vec{N}$ is aligned with $\vec{I}$, simplifying the coordinate transformation.

For each central pixel, we identify neighboring pixels within a radius of 99 arcminutes using the query_disc function from the HEALPix library (After 99 arc min most of the pixels have 0 values). This typically results in approximately 2800 neighboring pixels . For each neighbor $i$, we compute the angular separation $\theta_{i}$ between the pointing direction ($\hat{n}$) and the neighbor ($\hat{n}_{i}$), as well as the angle $\alpha_{i}$ between the scan direction $\vec{N}$ and the great-circle arc connecting the two pixels. These spherical coordinates are then converted to Cartesian offsets as

[TABLE]

The coordinates $(x_{i},y_{i})$ are mapped to discrete indices on a precomputed 2D beam grid $B(x,y)$, which is centered on the beam and typically spans a $100^{\prime}\times 100^{\prime}$ region with 2 arcsecond resolution. This beam map, shown in Figure 3, derived from Planck’s measured scan beams (e.g., HFI_ScanBeam_143-1a_R2.00.fits), encodes the instrument’s angular sensitivity as a function of beam-centered displacement. The corresponding value $B(x_{i},y_{i})$ is then used for calculating the effective beam matrix, $\widetilde{B}(\hat{n}_{j},\hat{n}^{\prime}_{i})$. This effective beam is then used for convolving all the skymaps.

While the above formulation describes the temperature response of the instrument, many Planck detectors are also sensitive to linear polarization. The detector response must therefore be extended to include the Stokes parameters $Q$ and $U$.

2.3.2 Polarization Response

In addition to temperature anisotropies, several channels of the Planck High Frequency Instrument also measure the linear polarization of the microwave sky. The polarization state of the radiation field is described by the Stokes parameters $Q(\hat{n})$ and $U(\hat{n})$, which represent the amplitude and orientation of the linear polarization. Unlike temperature, which is a scalar field on the sphere, $Q$ and $U$ depend on the choice of coordinate axes in the local tangent plane.

A polarization–sensitive detector measures the projection of the incident radiation along its polarization-sensitive axis. The signal recorded by a detector pointing toward the direction $\hat{n}$ can therefore be written as

[TABLE]

where $I(\hat{n})$ is the total intensity, $Q(\hat{n})$ and $U(\hat{n})$ are the Stokes parameters, $\rho$ is the polarization efficiency of the detector, and $\psi_{t}$ is the effective polarization orientation of the detector projected onto the sky at time $t$. In the Planck instrument model, the polarization efficiency is related to the polarization leakage parameter $\epsilon$ through $\rho=1-\epsilon$. The intrinsic detector polarization angles $\psi_{\mathrm{pol}}$ and leakage parameters $\epsilon$ used in our simulations are obtained from the Planck HFI Reduced Instrument Model (RIMO).

For the LFI 70 GHz channel, the Planck beam model directly provides Stokes beam components, which describe the full polarized beam response.

The polarization orientation $\psi_{t}$ varies during the scan as the spacecraft rotates. In our implementation, following the FEBeCoP formalism, $\psi_{t}$ is obtained by projecting the detector polarization direction onto the local sky basis at each pointing direction [26]. The detector polarization direction in the local beam frame is given by:

[TABLE]

where $\vec{I}$ and $\vec{Z}$ are the local orthonormal basis vectors defined in Section 2.3.1. The local spherical basis vectors $(\hat{e}_{\theta},\hat{e}_{\phi})$ define the tangent plane at the pointing direction $\hat{R}$. These are constructed from the unit pointing vector $\hat{R}$ as

[TABLE]

where $\hat{z}$ is the unit vector along the global $z$-axis. The effective polarization angle on the sky is then computed as

[TABLE]

This formulation directly provides the physically relevant detector orientation on the sky.

The trigonometric projections $\cos(2\psi_{t})$ and $\sin(2\psi_{t})$ determine how the Stokes parameters $Q$ and $U$ project onto the detector polarization axis. To incorporate polarization into the beam convolution, the detector response is applied directly during the effective beam construction. For each neighboring pixel $i$, the beam response weight $B(\hat{n},\hat{n}^{\prime}_{i})$ is multiplied by the corresponding polarization projection factors, producing three sets of weights corresponding to the temperature and polarization components:

[TABLE]

The observed Stokes parameters are obtained by applying the corresponding effective beam kernels to each component separately:

[TABLE]

In our implementation, these three sets of weights are stored separately during the effective beam construction stage and later applied independently to the input $I$, $Q$, and $U$ maps to generate the corresponding observed components.

The resulting convolution therefore accounts simultaneously for beam asymmetry, scan geometry, and detector polarization response, producing simulated observations that closely reproduce the characteristics of Planck polarization measurements.

The polarization convolution follows the same beam–sampling procedure described in Section 2.3.1. For each central pixel, neighboring pixels within a fixed angular radius are identified using the HEALPix query_disc routine. For the polarization kernels we restrict the neighbor search to a radius of $50$ arcminutes, which captures the dominant beam response while reducing the computational cost and storage requirements associated with the polarization beam matrices.

The angular displacement between the central pixel and each neighbor is projected onto the local scan-aligned coordinate system to obtain the beam-centered offsets $(x_{i},y_{i})$. These coordinates are then mapped onto the precomputed beam grid $B(x,y)$, which encodes the instrument response. The resulting beam weight $B(\hat{n},\hat{n}_{i}^{\prime})$ is subsequently multiplied by the polarization projection factors given in Eqs. (2.15) to obtain the effective polarization kernels $B_{Q}$ and $B_{U}$.

2.4 Generating the final skymaps

In Sections 2.2 and 2.3, we describe how we convolve the real beam to generate realistic simulated skymaps. For training and testing our neural network we have generated 2 sets of simulated skymaps:

• •

Circular beam convolution: Each frequency map is smoothed using a circular Gaussian kernel. This serves as a useful baseline for comparison.

• •

Real beam convolution: Precomputed beam response maps from Planck data are used, as described in Section 2.3.

For the circular simulation, we apply a Gaussian smoothing to each frequency channel using its corresponding full-width at half-maximum (FWHM) value. This is implemented using the healpy.sphtfunc.smoothing function from the HEALPix package, which smooths a map with a symmetric Gaussian beam. The FWHM and the detector noise values used for each frequency are summarized in Table 2.

For the symmetric beam convolution maps we add uncorrelated Gaussian noise (white noise) to the each pixel. As for the symmetric Gaussian beam convolution we have not considered the effect of the scanning, the noise variance at each pixel are kept the same.

For the real beam convolution we use frequency-dependent Planck scan beam models corresponding to the detectors listed in Table 1. This allows the simulations to capture the variation of beam shapes across different observing frequencies. As discussed earlier Planck satellite follows non-uniform scan strategy, where regions near the ecliptic poles receive significantly more observations than those near the equator (see Figure 2). As a result, pixels with high hit counts exhibit lower noise levels. To account for spatially varying noise, we scale the pixelwise noise variance inversely with the hit count—i.e., the number of times a pixel was observed by the satellite.

Combined, these noise and beam effects produce simulated maps that closely match the resolution and statistical properties of real Planck data. These simulations are essential for training the neural network and testing CMB foreground cleaning methods under conditions that realistically capture observational challenges.

2.5 Input Data Preparation and Patching

The full-resolution synthetic sky maps generated by our simulation pipeline are high-dimensional and span the full celestial sphere. However, standard convolutional neural networks are designed for flat Euclidean domains and do not directly support spherical convolutions. Implementing spherical CNNs would require a major code overhaul or the use of specialized frameworks such as DeepSphere [27]. To avoid this complexity, we adopt a patching strategy that enables training on planar projections.

Each map spans six frequency channels and is partitioned into 12 HEALPix regions, each represented as a square image of size $n_{side}\times n_{side}$. This enables the neural network to operate on localized sky regions, promoting more efficient training and inference.

The pixel values are indexed and arranged using the HEALPix nested ordering scheme. After patch extraction, the data is reshaped and stored in the ring ordering convention for compatibility with standard HEALPix-based tools. Each patch retains its correspondence to the original sky location and frequency band, ensuring spatial coherence and enabling evaluation of model performance both locally and globally.

This patching strategy enables the use of conventional CNNs without spherical convolution, while still preserving the spatial and multi-frequency structure of the input data.

3 Neural Network Framework for CMB Map Reconstruction

In this work, we perform separate reconstructions for temperature and polarization. For temperature, the network takes six frequency channels—70, 100, 143, 217, 353, and 545 GHz—as input and reconstructs the CMB intensity map. For polarization, we use three frequency channels—70, 217, and 353 GHz—which have significant polarization sensitivity, and train the network to recover the scalar $E$-mode polarization map derived from the Stokes $Q$ and $U$ components (see Section 2.1).

To recover the underlying CMB signal from noisy multi-frequency observations, we adopt a generative adversarial framework consisting of a U-Net–based generator and a convolutional discriminator. The generator produces high-fidelity estimates of the CMB maps, while the discriminator guides the training by assessing the realism of these reconstructions. The architecture and training strategy are described below.

3.1 GAN Architecture

Generative Adversarial Networks (GANs), first introduced by Goodfellow et al. [28], are a class of machine learning models that learn to generate data by training two networks in competition: a generator and a discriminator. The generator aims to produce outputs that resemble real data, while the discriminator learns to distinguish between real and synthetic examples. Through adversarial training, both networks improve iteratively, resulting in increasingly realistic generated outputs. GANs have found widespread applications in various fields, including image synthesis [29], natural language processing [30], and even astrophysics, where they have been used for tasks like weak lensing map generation [31]. While GANs offer powerful generative capabilities, they can be challenging to train due to issues like mode collapse and instability.

Our approach employs this adversarial framework to guide a U-Net-based generator towards producing CMB maps that are statistically consistent with real observations. Unlike standard GAN architectures that directly generate images from random noise, our generator takes noisy, foreground-contaminated multi-frequency sky maps as input. This allows the generator to learn the complex mapping from observed data to the underlying CMB signal. The generator architecture is based on the U-Net structure, which consists of an encoder and decoder network, connected via skip connections. The encoder progressively downsamples the input, extracting increasingly abstract features. The decoder then reconstructs the high-resolution CMB map, leveraging the features learned by the encoder and the skip connections to preserve fine-grained details.

The convolutional discriminator evaluates the realism of the generator’s outputs by comparing them to ground truth CMB maps. The discriminator acts as a learned loss function, guiding the generator to produce outputs that not only minimize pixel-wise reconstruction error but also adhere to the statistical properties of real CMB data. This interaction is formalized through a minimax objective: the generator aims to create outputs that the discriminator cannot distinguish from real maps, while the discriminator strives to correctly classify real versus generated maps. Through this process, the generator gradually improves to better approximate the true data distribution.

Our architecture leverages the encoder–decoder structure of U-Net for spatially precise reconstruction, while the discriminator enforces global and perceptual realism, resulting in more accurate and physically plausible outputs.

3.2 U-Net Architecture for CMB Reconstruction

The U-Net architecture, first introduced by Ronneberger et al. in 2015, employs a symmetric encoder–decoder structure that has proven highly effective for biomedical image segmentation and, more recently, for astrophysical map reconstruction [32]. The encoder hierarchically extracts spatial features, and the decoder reconstructs the high-resolution outputs by integrating these features via skip connections. In the context of Cosmic Microwave Background (CMB) reconstruction, U-Net–based models excel at denoising and deconvolving observed maps, accurately recovering underlying temperature and polarization fields from noisy, foreground-contaminated inputs [33, 34]. Figure  4 shows the overall U-Net architecture that we have used for CMB reconstruction. This particular U-Net is designed to take 6 frequency maps (3 frequency maps for polarization) with dimensions $n_{side}=1024$ as input.

3.2.1 Encoder

The encoder is structured as a sequence of nine down-sampling blocks, each designed to progressively reduce spatial resolution while increasing feature complexity. A schematic diagram of a single down-sampling block is shown in Figure 5(a).

Each down-sampling block starts with a 2D convolutional layer with a kernel size of $3\times 3$, a stride of 1, and padding of 1, ensuring that the spatial dimensions remain unchanged. These layers operate with $n_{\text{in}}$ input channels and $n_{\text{p}}$ output channels. The parameter $n_{p}$ acts as a control for the model size and is kept constant within each block. Specifically, the input channels ($n_{in}$) start at 4 in the initial block, increase to 8 in intermediate blocks, and reach 16 in the final blocks. A max-pooling layer with a kernel size of $2\times 2$ and a stride of 2 is applied immediately after the first convolutional layer, halving the spatial dimensions from $n_{side}\times n_{side}$ to $(n_{side}/2)\times(n_{side}/2)$.

Following max pooling, a second convolutional layer, also with a kernel size of $3\times 3$, processes the $n_{\text{p}}$ input channels to generate $n_{\text{out}}$ output channels, further refining the feature representations. Batch normalization is applied after the second convolutional layer to standardize feature activations, stabilizing and accelerating training [35]. Finally, a leaky rectified linear unit (LeakyReLU) activation with a slope of 0.2 introduces nonlinearity, allowing the model to capture complex patterns in the input data [36]. The outputs of each down-sampling block are stored for skip connections in the decoding stage, ensuring that detailed spatial information from earlier stages is preserved and utilized in reconstruction [32].

3.2.2 Decoder

The decoder reconstructs the high-resolution output map from the abstract features learned by the encoder. Each layer of the decoder, referred to as an upsampling block, performs spatial upsampling while integrating feature information from the encoder through skip connections. The process begins with a transpose convolution layer, which increases the spatial dimensions of feature maps. This layer operates with a kernel size of $4\times 4$, a stride of 2, and padding of 1, effectively doubling the spatial dimensions. The transpose convolution transforms feature maps from $n_{\text{in}}$ input channels to $n_{\text{out}}$ output channels, where $n_{\text{in}}$ and $n_{\text{out}}$ decrease progressively as the feature map size increases. The architecture of the upsampling blocks are shown in Figure 4.

In the first three upsampling blocks, batch normalization is omitted to preserve delicate spatial features and avoid potential artifacts during early-stage reconstruction. Instead, dropout with a rate of 0.2 is applied to prevent overfitting by randomly deactivating a subset of neurons during training [37]. The remaining subsequent upsampling blocks incorporate batch normalization [38] to ensure training stability. Each upsampling block concludes with a ReLU activation function to introduce non-linearity and support the reconstruction of complex spatial patterns [39].

Skip connections play a pivotal role in the decoder by linking feature maps from corresponding encoder layers to decoder layers. These connections concatenate the encoder features with the upsampled decoder feature maps, ensuring the preservation of fine-grained details. This integration allows the decoder to simultaneously utilize high-level contextual information from deeper layers and low-level spatial details from earlier layers [32].

3.2.3 Final Layer

The final layer of the U-Net employs a transposed convolution to restore the spatial dimensions of the output map to match those of the input. This layer takes the upsampled feature maps from the final decoder block and produces the final cleaned multi-channel CMB map. It uses a kernel size of $4\times 4$, stride of 2, and padding of 1, ensuring that the output resolution aligns precisely with the original input.

3.3 Discriminator Architecture

The discriminator network is designed to differentiate between real and reconstructed CMB maps, acting as a crucial component of the adversarial training framework. Its architecture consists of a series of downsampling blocks, similar to those in the U-Net encoder, followed by convolutional layers that compress input feature maps into a single scalar output representing the discriminator’s assessment. The input to the discriminator is formed by concatenating the true and generated CMB maps along the channel axis. Figure 6 illustrates the overall structure, showing the sequence of downsampling operations, convolutional layers, and zero-padding stages.

The discriminator utilizes four downsampling blocks that are architecturally identical to the first four downsampling blocks described in the U-Net encoder (see Figure 5(b)). These blocks progressively reduce the spatial dimensions and increase the feature channel depth. The first dwnsampling block takes $C_{in}$ + $C_{out}$ input channels—the combined channels of the input and target maps—and subsequent blocks continue increasing the feature depth. Batch normalization is omitted in the first block to preserve low level input features but is applied in the remaining blocks to stabilize training. Each downsampling blocks concludes with a leaky rectified linear unit (LeakyReLU) activation function with a slope of 0.2, introducing non-linearity and aiding in gradient propagation.

Following the downsampling blocks, the output passed through a pair of convolutional layers, each separated by a one-pixel padding stage. The first convolutional layer has a kernel size of $7\times 7$, a stride of 1, and no padding and maps the 32 input channels into 16 output channels. After this convolution, batch normalization and LeakyReLU are applied. The final convolutional layer uses a kernel size of $3\times 3$, a stride of 1, and no padding and it maps the 16-channel feature map to a single output channel with shape of $64\times 64$. This output serves as the discriminator’s final judgment on whether diffrent portions of the input map is real or generated.

3.4 Model Workflow and the Loss functions

The overall GAN workflow is illustrated in Figure  7. The input to the generator consists of multi-frequency sky maps contaminated by astrophysical foregrounds and instrumental effects, including beam convolution and spatially varying noise. For temperature reconstruction, the input is given by $x=\{T_{\mathrm{obs}}^{\nu}\},\quad\nu\in\{70,100,143,217,353,545\}\,\mathrm{GHz},$ while for polarization reconstruction, the input consists of frequency-dependent polarization maps, $x=\{E_{\mathrm{obs}}^{\nu}\},\quad\nu\in\{70,217,353\}\,\mathrm{GHz}.$ These maps are processed through the U-Net generator, which applies a series of downsampling and upsampling operations to produce a cleaned estimate of the underlying CMB sky. The encoder compresses the input into deep feature representations while storing intermediate outputs for skip connections. The decoder then reconstructs the high-resolution map, re-integrating spatial details via those skip connections. The discriminator evaluates whether the generator’s output resembles a real CMB map given the observed sky.

The discriminator receives two types of input pairs, both formed by concatenating the multi-frequency input sky map—contaminated with astrophysical foregrounds and instrumental noise—with a CMB map: once with the true CMB map and once with the generated map, along the channel dimension. This setup enables the discriminator to learn the conditional distribution of the true CMB given the observed sky, rather than directly comparing the generated and true maps. It learns to distinguish between real samples $(T_{obs},T_{cmb})$, where $T_{cmb}$ is the true CMB, and fake samples $(T_{obs},G(T_{obs}))$, where $G(T_{obs}(\hat{n}))$ is the output of the generator. This discriminator feedback is used to improve the generator through adversarial training.

Training is driven by a combination of pixel-wise reconstruction losses and an adversarial loss. The reconstruction term comprises both mean squared error ($L_{2}$) and mean absolute error ($L_{1}$) components:

[TABLE]

[TABLE]

The $L_{2}$ term penalizes large deviations, promoting overall pixel-wise accuracy, while the $L_{1}$ term helps preserve sharp features and mitigates blurring. However, relying on pixel-wise losses alone typically limits recovery to large-scale (low $l$) structures (see Section 4.1.1).

To encourage the generator to produce outputs that are not only accurate but also perceptually realistic, an adversarial loss is added based on the discriminator’s feedback:

[TABLE]

where $D(x,G(x))$ is the discriminator’s output when given the input sky map $x$ and the generated CMB map $G(x)$. The discriminator output is a two-dimensional activation map (e.g., $64\times 64$) that provides spatially resolved realism scores across the map.

The total objective combines all terms:

[TABLE]

where $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{\mathrm{adv}}$ control the relative contributions of each term.

These weights are chosen to bring the losses to a comparable numerical scale: specifically, $\lambda_{1}=100$ for the $L_{1}$ loss, $\lambda_{2}=1$ for the $L_{2}$ loss, and $\lambda_{\mathrm{adv}}=1000$ for the adversarial loss. For example, the unscaled $L_{2}$ loss is typically in the thousands, whereas the raw adversarial loss is around 1.4. Scaling the adversarial term by 1000 ensures it contributes significantly during training.

[TABLE]

The discriminator is trained simultaneously to maximize its classification accuracy—distinguishing between real and generated samples—while the generator is trained to minimize the adversarial loss and fool the discriminator. This adversarial training loop ensures that the generator not only reduces the reconstruction error but also learns to capture the distributional characteristics of real CMB maps.

In practice, we implement the adversarial losses using the BCEWithLogitsLoss function in PyTorch. This formulation is numerically more stable than standard binary cross-entropy loss, as it combines a sigmoid activation and log loss in a single step. Since our discriminator does not apply a final sigmoid activation,BCEWithLogitsLoss is appropriate for directly comparing its raw output logits to binary targets (1 for real, 0 for fake).

Tables 3 and 4 summarize the architectures of the U-Net generator and the discriminator, respectively. These architectures are optimized for input maps with a resolution corresponding to $n_{\text{side}}=1024$, enabling high-fidelity recovery of the CMB signal from contaminated observations.

3.5 Training

The dataset comprises 1000 simulated CMB maps, which were randomly split into a training set of 990 maps and a validation set of 10 maps. To model the galactic foregrounds accurately, we utilized the Python Sky Model (PySM) to simulate various dust and synchrotron emission models. Specifically, we employed 9 different dust foreground models (comprising 2 stochastic and 7 deterministic models) and 7 different synchrotron models (comprising 1 stochastic and 6 deterministic models).

To capture the variability inherent in the stochastic foregrounds, we used different random initial seeds to generate 10 distinct maps for each stochastic model. For the deterministic models, a single map was generated per model. We explicitly allocated the s1 and s2 synchrotron models, along with the d1, d2, and d4 dust models, to serve as the validation set for the foregrounds, reserving the remaining models for the training set. As detailed in Section 2.1, the PySM simulations initially yielded Stokes $I$, $Q$, and $U$ maps, which were subsequently converted into $T$ and $E$-mode maps. These maps have a shape of (6, 1024, 1024), corresponding to 6 frequency channels at a resolution of $n_{\text{side}}=1024$ (3 frequency channels for polarization maps). To facilitate training, each map was divided into 12 patches of size (6, 1024, 1024) using HEALPix nested indexing.

Training was conducted on the Param Himalaya cluster, utilizing NVIDIA V100 GPUs with 32 GB of memory. The experiments involved using 2 GPUs on a node, allowing for the parallel training of two patches simultaneously. The GAN was implemented in PyTorch, and the Adam optimizer was employed for both the generator and discriminator. The Adam optimizer [40] updates the parameters $\theta$ at iteration $t$ according to

[TABLE]

where $g_{t}=\nabla_{\theta}\mathcal{L}(\theta_{t-1})$ is the gradient of the loss function at step $t$. The exponential decay rates, $\beta_{1}$ and $\beta_{2}$ and numerical stability constant $\epsilon$ are set to the commonly used values

[TABLE]

During training, the learning rate ($\alpha$) for the generator was initially chosen to be of order $10^{-2}$ and $10^{-3}$ for the discriminator, slowly decreasing them to around orders of $10^{-5}$ and $10^{-7}$ respectively as the loss became lower and less stable.

Multiple training runs were performed for each patch to optimize the model parameters. The average training time for one patch was approximately 40 hours when simulating circular beam and extended to 4-5 days for real beam simulations. The maps convolved with real beam are trained against HEALPix skymaps convolved with an effective gaussian circular beam. During training, the loss curves exhibited instability, generally following a pattern of initial descent followed by an increase, eventually stabilizing around 0.69 for the discriminator loss. The L1 loss and MSE loss showed exponential decay behavior.

We experimented with transferring weights from previously trained patches to initialize training for new patches. However, this approach did not yield a significant improvement in training time or loss behavior.

4 Simulations and Results

We perform three sets of simulations to assess the performance of our CMB reconstruction pipeline: one that uses a simplified circular beam for the temperature only skymap, and two that incorporate the Planck satellite’s actual scanning strategy and beam pattern for temperature and E-mode polarization skymaps respectively. The goal is to assess how well the model generalizes under increasingly realistic observational conditions. We also analyze the effect of adversarial training and the impact of model size.

4.1 Results with circular scan

In this section, we present our analysis based on the circular beam scan on CMB temperature skymap. For this convolution, no actual scanning pattern was considered; as a result, the noise variance is uniform across all pixels (see Sec. 2.4 for details). Figure 8 shows the results from our analysis. Green plot shows the reference $C_{l}$ that is used for generating the simulated CMB maps. The orange curve shows the $C_{l}$ from the outputs from our neural network. The comparison indicates that the network successfully recovers the input spectrum almost up to $l=1500$, beyond which the reconstructed power gradually drops. We show the results in both linear and logarithmic scales to allow readers to assess the agreement across both high and low multipole ranges. As mentioned earlier, $n_{p}$ is a modeling parameter, and for this analysis we set $n_{p}=8$.

4.1.1 Effect of Adversarial Training

To evaluate the impact of including a discriminator in the training process, we compared the power spectra of CMB maps reconstructed using only a generator model versus those reconstructed using both a generator and a discriminator model. Figure 8 illustrates these comparisons. The plot on the left presents the power spectra using a linear scale for the x-axis (multipole moment l), highlighting the difference at higher multipole moments. We can see that the power spectra produced without a discriminator (shown in blue) start to deviate from the original data after approximately l = 900. The inclusion of the discriminator helps maintain the power spectra closer to the original even at these higher multipole moments. The plot on the right provides the same comparison but uses a logarithmic scale for the x-axis to better visualize differences at lower l values. Here, it is evident that the model without a discriminator yields higher values than the original data for low l values. In contrast, the power spectra generated with a discriminator show a more accurate representation of the original data, particularly at these lower l values. These plots clearly demonstrate that the inclusion of a discriminator in the training process helps improve the reconstruction quality by better aligning the power spectra with the original data, especially at higher multipole moments where the generator alone fails to maintain consistency. The discriminator acts as a regularizing mechanism, preventing the power spectra from deviating significantly from the original at critical points.

4.1.2 Impact of Model Size on Reconstruction Performance

To explore the effect of model size on reconstruction performance under the circular beam assumption, we compared the power spectra of CMB maps reconstructed using different model sizes, corresponding to hyperparameter values $n_{p}$ of 2, 4, 8, and 16. Figure 9 illustrates these comparisons across various scales. The logarithmic scale plot (top-right) shows that all model sizes perform similarly well for low l values, suggesting they effectively capture the early, low-frequency features of the CMB power spectra. However, the linear scale plot (top-left) and the plot focusing on higher multipoles (bottom-left) reveal that as the model size increases, the accuracy of the power spectra improves for high l values. Specifically, models with sizes n=8 and n=16 maintain a closer alignment with the original data up to approximately l = 1700. In contrast, the smaller models diverge earlier: the $n_{p}$=2 model deviates significantly around l = 1100, and the $n_{p}$=4 model diverges around l = 1300. The bottom-right plot shows the relative error between the reference and the reconstructed $C_{l}$s, given by the absolute value of $(C_{l}^{Ref}-C_{l}^{Recons})/C_{l}^{Ref}$. The plot shows that the improvement is minimal when we change the hyperparameter $n_{p}=8$ to $n_{p}=16$. This analysis indicates that increasing model size allows the capture of more detailed features, leading to more accurate reconstructions at higher l values, better capturing high-frequency components. Nevertheless, increasing model size introduces trade-offs regarding complexity, leading to longer training times and higher computational costs. While models $n_{p}$=8 and $n_{p}$=16 perform similarly well up to l = 1700, the $n_{p}$=8 model offers a good balance between complexity and performance, making it an optimal choice in this context.

4.2 Reconstruction Results with Real Beam Convolution for temperature map

In this analysis, we convolve the foreground contaminated CMB temperature maps using Planck beam and scan pattern as described in Sec. 2 and then add the noise to generate the simulated realizations for all 6 frequency channels. Training proceeded using the discriminator-based approach with $n_{p}=8$. We have presented the results from our analysis in Figure 10. The top-left plot shows the full power spectrum covering the entire multipole range ($l=2$ to $2500$). We can see that the network can recover the spectra significantly well up-to multipole $l\sim 1000$, after which the power starts falling off. The next plot, i.e. the plot on the top-right plot shows the same power spectrum but the x-axis is set to the log scale (semilogx plot). This allows the readers to understand the behavior of the power-spectrum at low multipoles. The third figure highlights the behavior of the reconstruction specifically within the multipole range $l=200$ to $1000$. Lastly, the bottom-right figure focuses on the performance at higher multipole moments ($l>1000$), a regime where the model faces significant loss.

These figures collectively demonstrate that the incorporation of the real beam convolution into the training pipeline, coupled with the discriminator, effectively preserves the CMB power spectrum characteristics across the full range of multipoles, while also ensuring accurate reconstruction in both the intermediate and high-l regimes.

Figures 5 and 6 present a detailed visual assessment of the network’s CMB reconstruction performance for each of the 12 patches. The first 4 patches are from the upper hemisphere, the next 4 patches are from the equatorial region and the final 4 are from the lower hemisphere. Each patch is displayed across four rows: the noisy input map containing foregrounds, the predicted CMB signal from the generator network, the true underlying CMB signal, and the difference between the predicted and true CMB (Predicted - True). Here we would like to point out that there are total 6 frequency channels that are used as inputs. In the plots we display only the $70$GHz channel. For plotting the input, true and predicted maps we set the temperature range between $-200\mu K$ and $200\mu K$. For plotting the differences no such range is set. The change in the color pallet is due to difference in range. Our plots show that the foreground removal algorithm works perfectly near the polar region. The difference in the signals are less than $2\mu K$, (approximately $1\%$) of the maximum pixel temperature in the true map. Near the Galactic center, some residual signals are visible, particularly in the plots for patches 5, 6, 7, and 8. However, even close to the Galactic plane, the residual signal remains below $\sim 2$–$6,\mu\mathrm{K}$ (approximately $2\%$), except for a few isolated pixels where the difference reaches $8$–$10,\mu\mathrm{K}$.

4.3 Reconstruction Results with Real Beam Convolution for Polarization Maps

Training of the network proceeded using the same discriminator-based approach with $n_{p}=8$, targeting the reconstruction of the uncontaminated CMB polarization field. However, two key modifications were made for the polarization analysis. First, the generator architecture was adjusted to accept 3 input frequency maps, in contrast to the 6 channels used for the temperature reconstruction. Second, an additional preprocessing sanitization step was required to appropriately handle masked pixels present in the healpy maps prior to feeding the data into the network.

We present the power spectra of the reconstructed polarization maps, focusing on the $E$-mode ($EE$) spectrum, in Figure 11. The top-left panel shows the full $EE$ power spectrum covering the entire multipole range ($l=2$ to $2500$). The network successfully recovers the polarization spectra up to multipole $l\sim 1100$. Beyond this scale, the reconstructed power begins to attenuate,

To comprehensively evaluate the reconstruction quality under these realistic observational and scanning conditions, we decompose the spectra into specific multipole regimes. The top-right plot displays the same $EE$ power spectrum on a semilogarithmic scale (semilogx) to emphasize the behavior at low multipoles. The bottom-left figure highlights the intermediate regime within the multipole range $200\leq l\leq 1100$. Lastly, the bottom-right figure focuses on the high-multipole performance ($l>1000$), illustrating the regime where the reconstruction faces significant signal loss due to beam smearing and noise domination.

The raw difference in the signals are less than $0.5\mu K$, (approximately $2\%$) of the maximum pixel temperature in the true map. The residual signals Near the Galactic center, are lower compared to the temprature maps, except for a few isolated pixels where the difference reaches $1\mu\mathrm{K}$.

4.4 Exploring BipoSH coefficients

In the standard cosmological model, the Universe is assumed to be statistically isotropic (SI). The CMB temperature anisotropies can be expanded in spherical harmonics as

[TABLE]

where $a_{lm}$ are the corresponding expansion coefficients. If these anisotropies constitute a Gaussian random field—consistent with high-precision Planck measurements —their statistical properties are fully determined by the two-point correlation function. Under SI, this is equivalently expressed in harmonic space as

[TABLE]

where $\langle\cdots\rangle$ denotes an ensemble average over CMB realizations.

However, in a non-SI universe, to capture the full description of the two-point correlation function, the BipoSH spectra are required [41]. The BipoSH spectra $A_{ll}^{LM}$ are defined by the expression [42],

[TABLE]

where $\Pi_{L}:=(2L+1)^{1/2}$ and $C_{lml^{\prime}m^{\prime}}^{LM}$ are the Clebsch-Gordan coefficients. While the BipoSH spectrum $A_{ll}^{00}=C_{l}$ is the standard angular power spectrum, detecting non-vanishing power in the remaining BipoSH spectra $A_{ll}^{LM},L,M\neq 0$ forms the basic criteria for probing SI violation.

Non-circular beams, non-uniform scan patterns, and anisotropic noise introduces BipoSH coefficients in the observed CMB sky. Thus, if residual foregrounds, anisotropic noise, or uncorrected beam and scan effects remains, they can lead to isotropy violation in the reconstructed skymaps [43, 44, 45, 46, 47, 48, 49].

Figure 12 shows the $A^{20}_{ll}$ and $A^{20}_{ll-2}$ BipoSH coefficients for the maps in ecliptic coordinates from the temperature skymaps. Because the scan pattern is symmetric across the upper and lower ecliptic hemispheres, these two coefficients are expected to capture the dominant contributions in the scanned sky maps. In the top plots we have the BipoSH coefficients from all the maps after beam convolution (before adding any noise and foreground). The reference maps, generated with HEALPix, are statistically isotropic and therefore yield BipoSH coefficients consistent with zero. In the plot we can see that all the data points are consistent with [math]. In contrast, the scanned and beam convolved maps display clear nonzero values, reflecting isotropy violations induced by the noncircular beam and scan strategy. The amplitude of the BipoSH coefficients is determined by the beam shape [46]. To emphasize the BipoSH coefficients from both the reference and reconstructed maps in the bottom panel, we show the BipoSH coefficients for the reference (green) and reconstructed (orange) maps, together with those corresponding to the smallest beam, i.e., the 545 GHz channel.

The addition of anisotropic noise and Galactic foregrounds should, in principle, further enhance such violations. However, the reconstructed maps show no such signatures, with the BipoSH coefficients returning to values consistent with zero. This demonstrates that the neural reconstruction not only suppresses residual foregrounds and anisotropic noise but also effectively removes the distortions introduced by the noncircular beam and scanning pattern—something conventional methods are unable to achieve. This makes the neural reconstruction approach significantly more powerful than existing techniques.

In Figure 13, we present similar plots the E-mode polarization maps at different frequencies. We observe that the BipoSH signatures induced by the beam shapes and the scan strategy are completely removed. However, the standard component-separation algorithm fails to eliminate the SI signal introduced by the beam convolution and the scan pattern.

Various previous studies have also noted that splitting maps into patches and training a neural network on each patch can introduce edge effects, potentially leading to anomalous signals in the final skymap [50]. Any such effect would leave an imprint in the BipoSH coefficients. As the BipoSH coefficients from the final maps are consistent with zero, we conclude that no statistically significant signal due to patch-wise reconstruction is present in the final maps.

5 Discussion

The results presented provide important insights into the performance of the U-Net architecture, particularly when augmented with a discriminator and subjected to realistic observational conditions. Key findings emerge regarding the discriminator’s role, the impact of model size, and the challenges posed by real beam convolution.

Our comparisons clearly demonstrate the benefit of including a discriminator during training (Section 4.1.1). It significantly enhances reconstruction quality by improving the consistency of the reconstructed power spectra with the original data. Figure 8 shows noticeable improvements at low multipole (up to $l\approx 100$) as well as at high multipoles ($l>900$) compared to the generator-only model.

Concerning model size (Section 4.1.2), initial experiments with Gaussian beam convolution suggested that performance gains diminish beyond $n=8$, with models like $n=8$ and $n=16$ maintaining good accuracy up to $l\approx 1700$ (Figure 9). The analysis incorporating a real beam convolution (Section 4.2) reveal a more stringent performance limit under these observational conditions. Using the $n=8$ model, chosen based on the earlier trade-off analysis, we found that fidelity with the ground truth power spectrum is maintained reasonably well only up to $l\approx 1000$ (Figure 10). Beyond this multipole, the reconstructed spectrum shows increasing deviation, indicating that the model struggles to accurately recover the CMB signal at smaller angular scales when faced with beam smoothing effects.

It is important to note that for circular beam convolutions, we use different FWHM values for different frequency channels, as listed in Table 2. The noise variance is uniform across the sky. In contrast, for the real scan, the noise variance varies spatially, depending on the hit count distribution—further increasing the complexity of accurate reconstruction.

One limitation of this analysis is that the generation of realistic sky maps is computationally intensive, particularly for larger model sizes or higher-resolution maps. In this work, we use maps with a resolution of $n_{\text{side}}=1024$, whereas the native Planck resolution is $n_{\text{side}}=2048$. Generating maps at the full Planck resolution would substantially increase both computation time and storage requirements, primarily because of the beam convolution and map synthesis steps. In addition, training the neural network on data at this resolution would further intensify the computational load. For the polarization analysis, we have restricted ourselves to E-mode polarization and to three frequency bands, constrained by available computational resources. Including all frequency maps and both polarization modes would enable us to investigate polarization leakage and assess whether our neural network–based method can effectively mitigate it. However, such comprehensive studies are computationally expensive and have therefore been postponed to future work.

Future improvements to the neural network should focus on enhancing performance beyond $l\approx 1000$ under realistic beam conditions. While implementing model quantization or pruning could address computational efficiency, and advanced training algorithms might reduce runtime, the core challenge lies in improving high-$l$ fidelity. This could involve enhancing the model architecture (e.g., incorporating attention mechanisms, exploring different convolutional approaches suited for deconvolution tasks), refining the loss function to be more sensitive to high-$l$ power spectrum features, or developing multi-stage refinement processes.

These findings underscore the crucial need to balance reconstruction quality with computational feasibility and, critically, to develop deep learning models robust to realistic observational effects like beam convolution. Achieving accurate recovery of cosmological information encoded at smaller angular scales ($l>1000$) remains a key area for future work in applying neural networks to CMB analysis.

6 Conclusion

In this work, we demonstrate that a U-Net–based generator combined with a discriminator network offers a powerful framework for reconstructing CMB maps from simulated data. Several recent studies have applied machine learning and neural networks to foreground cleaning in temperature and polarization skymaps [13, 51, 52, 53]. However, these works do not take into account the effects of non-circular beams and the instrumental scan pattern. Others use simulated skymaps from Planck [54], which inherently include the beam effect, although little attention has been given to understanding its impact. The instrumental effects can leave imprints in statistical measures such as the Bipolar Spherical Harmonic (BipoSH) coefficients [55, 46]. In contrast, our analysis fully incorporates both the non-circular beam effects and the scan strategy, resulting in a more sophisticated and realistic treatment of the data. For the first time, we demonstrate that a neural network–based algorithm can not only clean the foregrounds but also effectively remove the distortions introduced by the noncircular beam and scan pattern. Moreover, most researchers rely on simple U-Net–based architectures for foreground cleaning. In contrast, our approach incorporates a discriminator network alongside the generator, which enhances training efficiency and accelerates convergence.

To reduce computational cost, we adopt a patching strategy in which the full-sky map is divided into smaller, manageable segments. This approach significantly reduces the memory and time requirements compared to training on the entire sky at once. As map resolution increases, full-sky convolutional training becomes increasingly time-intensive and impractical. However, the patching method helps mitigate this issue. For example, a high-resolution map with $n_{side}=2048$ can be partitioned into 48 patches of size $1024\times 1024$, allowing parallel training of each patch if sufficient computational resources are available. By forgoing full-sky spherical convolutional techniques like DeepSphere [27] in favor of this patchwise approach, we achieve a practical compromise between computational feasibility and foreground cleaning accuracy. Importantly, our analysis shows that the neural network developed in this work does not introduce any statistically significant edge effects, countering concerns raised in previous studies and reinforcing the robustness of the patching strategy [50].

Previous studies on neural-network–based CMB foreground cleaning have primarily relied on the mean squared error (MSE) or $L_{2}$ loss function alone. Although combining multiple loss functions is a common practice in machine learning, it has not been explored in this context. In our analysis, we incorporate a combination of loss functions and find that including an $L_{1}$ term alongside the $L_{2}$ loss leads to more stable training convergence.

While applying neural networks for foreground removal on the CMB sky may not always offer significant speed advantages over traditional statistical methods, particularly given that there is only one CMB sky, the approach holds substantial promise in other domains. In areas such as HI intensity mapping, where large volumes of data across multiple frequency channels are involved, neural networks can drastically accelerate data processing. For instance, tasks like RFI removal, calibration, and solar contamination mitigation in Tianlai observations currently demand considerable computational resources, often taking a full day to process a single day’s data using hundreds of processors [56, 57, 58]. Training a neural network for such tasks is computationally expensive upfront, but once trained, the inference time is minimal, enabling fast, large-scale data cleaning.

Beyond cosmology, similar architectures can be adapted for a variety of applications, including weather data analysis and other time-series or imaging datasets where traditional methods are limited by processing time [59]. Thus, the proposed neural network framework is not only valuable for CMB analysis but also offers a scalable, efficient solution for broader scientific data challenges.

Data and Code Availability

The codes for generating simulated sky maps by convolving the real beam and Planck scan pattern, can be found in https://github.com/SJaiswal711/cmbnn. The codes for training and testing the U-Net GAN model is available at https://github.com/Obasho10/cmbnn.

Acknowledgments

Computational resources for this study were provided by the Param Himalaya Supercomputing Facility at IIT Mandi. S.J. acknowledges financial support from Dr. Nirmalya Kajuri’s SERB-funded project (Project No. IITM/SERB/NK/395) and Dr. Vikram Khaire’s SETI-funded project (Project No. PHY2425003SETIKHAI) during the course of this work. S.D. thanks Aditya Rotti and Nidhi Pant for letting us use their code for calculating BipoSH spectrum.