ForSE+: Simulating non-Gaussian CMB foregrounds at 3 arcminutes in a stochastic way based on a generative adversarial network

We present ForSE+, a Python package that produces non-Gaussian diffuse Galactic thermal dust emission maps at arcminute angular scales and that has the capacity to generate random realizations of small scales. This represents an extension of the ForSE (Foreground Scale Extender) package, which was recently proposed to simulate non-Gaussian small scales of thermal dust emission using generative adversarial networks (GANs). With the input of the large-scale polarization maps from observations, ForSE+ has been trained to produce realistic polarized small scales at $3\arcmin$ following the statistical properties, mainly the non-Gaussianity, of observed intensity small scales, which are evaluated through Minkowski functionals. Furthermore, by adding different realizations of random components to the large-scale foregrounds, we show that ForSE+ is able to generate small scales in a stochastic way. In both cases, the output small scales have a similar level of non-Gaussianity compared with real observations and correct amplitude scaling as a power law. These realistic new maps will be useful, in the future, to understand the impact of non-Gaussian foregrounds on the measurements of the cosmic microwave background (CMB) signal, particularly on the lensing reconstruction, de-lensing, and the detection of cosmological gravitational waves in CMB polarization B-modes.

1 Introduction

The cosmic microwave background (CMB) contains valuable information about the initial conditions in the early Universe and physical processes that happened during its propagation from the last scattering surface to our observation (Hu & Dodelson, 2002). In addition to the great success of observations and interpretation of the CMB total intensity signal, mainly sourced by scalar perturbations (Hinshaw et al., 2013; Planck Collaboration VI, 2020), the focus of experimental efforts has now been transferred to the observations of the CMB polarization signal. In particular, the quest for the curl or B-mode component in the CMB polarized signal has drawn the most attention nowadays, as it can be a unique probe of primordial gravitational waves (PGWs, i.e., tensor-mode cosmological perturbations) (Zaldarriaga & Seljak, 1997; Kamionkowski et al., 1997; Kamionkowski & Kovetz, 2016). These tensor perturbations are predicted by the families of inflationary theories, which describe a history of accelerated expansion in the very early Universe (Guth, 1981). Great efforts are being made with many ongoing and future telescopes, including ground-based ones, such as PolarBear (POLARBEAR Collaboration, 2022), SPT (SPT Collaboration, 2023), ACT (Madhavacheril et al., 2024), SO(Ade et al., 2019), CMB-S4 (Abazajian et al., 2016), BICEP (The BICEP/Keck Collaboration, 2022), and AliCPT (Li et al., 2018), balloon-borne ones such as SPIDER (Spider Collaboration, 2022), and telescopes on satellites, like LiteBIRD (LiteBIRD Collaboration et al., 2023). The tightest constraint on the amplitude of the PGW, in terms of the tensor-to-scalar ratio, $r\equiv A_{T}/A_{S}$, is $r<0.032$ (Tristram et al., 2022), at a $95\%$ confidence level.

While the scalar perturbations do not contribute to the CMB B-modes and only source the parity-even E-modes, the gravitational lensing effect from the forming cosmological structures along the propagation of CMB photons, with a typical deflection angle of a few arcmins, converts a small portion of E-modes into B-ones. These lensing-induced B-modes dominate over the PGWs at arcminute angular scales (Hu, 2000; Planck Collaboration VIII, 2020). Some inflationary models also predict nearly Gaussian statistics of the CMB anisotropies with small deviations known as primordial non-Gaussianity (PNG), which can be used to further constrain the vast space of inflationary theories that go beyond a simple single field slow roll scenario (Bartolo et al., 2004; Planck Collaboration IX, 2020). However, the lensing effect will also induce non-Gaussianity into the distribution of CMB photons by distorting the paths and coupling different multipole bins of power spectra.

Despite being a contaminant to PGWs and PNG, CMB lensing carries cosmological information on processes occurring during structure formation, on neutrino masses (Allison et al., 2015), and Dark Energy (Acquaviva & Baccigalupi, 2006). Effective techniques capable of reconstructing the lensing field have been developed, mainly relying on detecting the non-Gaussianity presented in the observed CMB maps (Hu & Okamoto, 2002; Maniyar et al., 2021). The latest detection of lensing potential comes from the ACT experiment, which achieves a $43\sigma$ detection of the CMB lensing signal (Qu et al., 2024).

The primary source of contamination to CMB B-modes arises from astrophysical processes causing diffuse emissions within our own Galaxy. The two main polarized components are represented by the thermal dust and synchrotron emission, dominating at high ($\gtrsim 70$ GHz) and low ($\lesssim 70$ GHz) frequencies, respectively, exceeding the CMB B-mode signal in any frequency and any position on the sky (Krachmalnicoff et al., 2016, 2018). They are caused by dust grains heated by starlight, exhibiting a quasi-thermal emission at about 18 K, and cosmic ray electrons spiraling around the lines of the Galactic magnetic field, respectively. The class of data analysis algorithms capable of extracting the cosmological B-mode anisotropies out of a multi-frequency dataset is known as component separation. Several implementations exist, exploiting the different properties that CMB and foregrounds have, such as frequency dependency and spatial distribution. Depending on the assumptions made about the foregrounds, these methods can generally be classified into blind (Yao et al., 2018; Santos et al., 2021; Zhang et al., 2019; Delabrouille et al., 2009) and parametric ones (Stompor et al., 2009; Planck Collaboration IV, 2020; BeyondPlanck Collaboration I, 2023).

Unlike the CMB, diffuse foregrounds generically exhibit a large degree of non-Gaussianity. This comes primarily from large-scale observations (Coulton & Spergel, 2019; Allys et al., 2019) and is expected to also be true at smaller scales, since dust grains are highly locally distributed and the magnetic field in the diffuse interstellar medium is highly turbulent. The foreground non-Gaussianity introduces coupling between different modes in the angular power spectra, which induces non-diagonal terms in the covariance matrices of the observed signal. On the other hand, covariance matrices are usually treated as diagonal under the assumption that all the sky components behave like Gaussian random fields. Additionally, approaches relying on power spectra will inherently ignore the non-Gaussianity, as the two-point functions are insufficient in capturing its presence. Therefore, the existence of non-Gaussianity in the Galactic foregrounds will potentially induce systematic errors in primordial B-mode detection.

Abril-Cabezas et al. (2024) claim that the constraints on $r$ are robust against the non-Gaussianity in the polarized dust at the power spectra level, but they only focus on $30\leq\ell\leq 300$, the angular scales at which the PGWs cause a maximum B-mode contribution known as recombination bump. At smaller scales where lensing dominates, non-Gaussian foreground will possibly cause bias to lensing reconstruction (Beck et al., 2020) and de-lensing. Moreover, non-Gaussianity also affects PNG measurements (Cabella et al., 2010).

The main reason for neglecting non-Gaussianity is that observations of foregrounds at arcminute scales are lacking, and so is our capability of generating realizations of simulations of their patterns. Current observations of polarized foregrounds are limited to degree scales over a large sky fraction and can only reach a higher resolution at sub-degree or arcminute scales for portions of the sky (Planck Collaboration XI, 2020; Bernardi et al., 2004; Remazeilles et al., 2015).

The way to simulate foregrounds relies on data and phenomenology. Foreground observations at a given frequency are used as templates and extrapolated to other frequency bands by making assumptions about their spectral energy distributions (SEDs). These templates are mostly based on the observations by the Planck (Planck Collaboration IV, 2020) and Wilkinson Microwave Anisotropy Probe (WMAP) satellite observations (Remazeilles et al., 2015; Bennett et al., 2013), which characterize the signal down to a limited angular resolution (around $1^{\circ}$). There exist several widely used packages that simulate Galactic foreground maps, such as the Python Sky Model (PySM, Thorne et al. (2017)) and Planck Sky Model (Delabrouille et al., 2013). The latest version of PySM3 package (Zonca et al., 2021; PanEx Collaboration, 2024) uses foreground templates from the latest available observations and extrapolates them to arcminute angular scales. These packages usually fit a power law to the observed large-scale power spectra of foregrounds and extrapolate to small scales with Gaussian realizations of them. Although commonly used, these packages are not able to produce foreground small scales with significant non-Gaussianity comparable to the observed large scales.

Clark & Hensley (2019) provides a data-driven framework to construct three-dimensional Stokes parameter maps of thermal dust emission in position–position–velocity space using only neutral hydrogen data on the basis that HI is strongly correlated with the dust in the diffuse interstellar medium (Lenz et al., 2017). By integrating over the velocity space, they obtained the polarized dust emission maps over the full sky at a resolution of $16.2\arcmin$, corresponding to that of HI data, which are in good agreement with the Planck observed 353 GHz dust maps in terms of several physical properties such as the polarization fraction and power spectra, although they do not have a thorough discussion about the non-Gaussianity contained in their maps.

Hervías-Caimapo & Huffenberger (2022) uses a large number of filaments in the distribution of the thermal dust grains together with the large-scale template from data to reproduce the main features of the Planck 353 GHz map. These include the power spectrum slopes of intensity and polarization maps, the ratios between EE, BB, and TE power spectra, and the level of non-Gaussianity in the total intensity map, which can be controlled by the density of filaments. When focusing on scales of arcminutes or tens of arcminutes, $\ell=300-1200$, however, the Planck 353 GHz total intensity map has more non-Gaussianity with respect to the generated small scales from the model based on filaments.

Diffuse foreground models are also generated by exploiting magnetohydrodynamic (MHD) simulations. They model the physical processes of the interstellar medium, such as heating and cooling of the gas and the interaction between dust grains and the turbulent magnetic field. The thermal dust emission can be obtained by integrating the emission of dust grains along the line of sight (Padoan et al., 2001; Planck Collaboration Int. XX, 2015; Kim et al., 2019). The simulated maps are thus non-Gaussian because of the MHD processes and can reach small scales according to the resolution of the MHD simulation itself. However, the MHD simulations can only reproduce the statistical properties of the Galaxy, and thus fail to generate the specific morphology of Galactic foreground emission. Another important limiting factor of MHD simulations is that they are computationally expensive, especially in achieving a high resolution.

An additional technique consists of simulating the foreground by exploiting innovative algorithms capable of modeling the emission pattern at a high resolution on the basis of the observed one at moderated and low angular scales. The ForSE model introduced by Krachmalnicoff & Puglisi (2021) (hereafter KP2021) is able to generate small scales, up to $12\arcmin$, starting from the low-resolution polarized dust observations. It utilizes a generative adversarial network (GAN), which is trained to inject small-scale features with statistical properties like the ones observed at a high resolution in the intensity maps.

In this study, we further extend the ForSE algorithm to: i) produce dust polarization maps at $3\arcmin$ resolution, which is vital for checking the stability of lensing reconstruction and de-lensing algorithms in CMB observations, ii) stochastically generate multiple realizations of small-scale features, which are most important in estimating the uncertainty associated with foreground variations.

Additional uses of deep generative models for dust simulations exist in the literature. For example, Aylor et al. (2021) use GANs to generate simulated total intensity maps from the observed Planck generalized needlet internal linear combination (GNILC) map at 353GHz. Other common generative models, such as variational autoencoders (Thorne et al., 2021) or diffusion models (Heurtel-Depeiges et al., 2023), are also used.

The outline of the paper is as follows. In Section 2, we summarize the first version of the ForSE algorithm and present the extension and methodology used to achieve the new version operating at arcminute scale. In Section 3, we describe the pre-processing, training, and post-processing of the new model and also present our validation procedure. In Section 4, full-sky maps are presented and compared with maps from the latest version of PySM3 package. Conclusions are given in Section 5.

2 The ForSE+ model

In this section, we briefly review the basic structure and assumptions of the ForSE algorithm presented in KP2021. We then introduce our new version of the code, ForSE+, which allows one to generate maps of the thermal dust emission with non-Gaussian structures at $3\arcmin$ angular resolution, in both a deterministic and a stochastic way.

2.1 Review of the ForSE model: From $80\arcmin$ to $12\arcmin$

As has already been mentioned, the ForSE model, introduced and validated in KP2021, is based on GANs and allows one to produce non-Gaussian full-sky maps of polarized thermal dust emission at an angular resolution of $12\arcmin$ from low-resolution Planck observations at $80\arcmin$.

Generative adversarial networks are a particular family of networks (Goodfellow et al., 2014) whose characteristic feature is to be composed of two sub-networks called Generator ($G$) and Discriminator ($D$), which are trained to compete against each other. In practice, the goal of $G$ is to produce new images that are compared by $D$ with a set of real images that are the training set. Once the training of the two sub-networks is done in an adversarial way, $G$ is able to produce images that have the same statistical properties as those belonging to the training set, in such a way that mock and real images are no longer distinguishable by $D$.

In the ForSE implementation, the input to $G$ are images at a low resolution ($80\arcmin$) of the thermal dust emission observed by the Planck satellite and processed through the GNILC method (Planck Collaboration IV, 2020). The output of $G$ are small-scale features at $12\arcmin$, which are then compared by $D$ with real observations in total intensity at the same angular resolution. We note that different training is performed for total intensity and Stokes $Q$ and $U$ maps, but always having as the training set images of small scales $\tilde{m}_{12^{\prime}}^{I,20^{\circ}}$ in total intensity. KP2021 therefore assumed that thermal dust statistical properties of small scales in polarization are the same as for Stokes $I$ maps. We also rely on that assumption in this work.

In KP2021 and in this work, the following definition of “small scales” is used. Let $M_{TOT}$ be a foreground map at a given angular resolution, which can be seen as the sum of two maps containing, respectively, only large- and small-scale structures:

$$M_{TOT}=M_{LS}+M_{SS}.$$ (1)

Assuming that the map encoding small-scale features, $M_{SS}$, is modulated by the large scales, $M_{LS}$ – that is, $M_{SS}=M_{LS}\cdot m_{ss}$ – we have

$$M_{TOT}=M_{LS}+M_{LS}\cdot m_{ss}=M_{LS}\cdot\tilde{m}_{ss},$$ (2)

where $\tilde{m}_{ss}=m_{ss}+1$ represents the small-scale map generated by network $G$, having as an input $M_{LS}$.

Although there exist neural networks designed to work on the sphere (Krachmalnicoff & Tomasi, 2019), our GAN deals with flat two-dimensional images. For this reason, we had to project the input maps onto flat patches, and project output patches back onto the sphere, after the application of the trained GAN model. We used the same projection strategy as that described in the appendix of KP2021, which divides the GNILC thermal dust Stokes $Q$ and $U$ maps at 80$\arcmin$ with $N_{side}=2048$ into 174 square patches that have $320\times 320$ pixels and a physical side length of $20^{\circ}$. We note that this projection on flat patches and reprojection on the sphere can introduce distortions in the final full-sky map. In order to mitigate this effect, the flat patches overlap each other. We estimate that the final level of distortion induced by our procedure is always less than 7% of the signal (around 2% on average).

From now on we will use $\tilde{m}_{z^{\prime}}^{X,y^{\circ}}$ and $M^{X,y^{\circ}}_{z^{\prime}}$ to refer to the small-scale (output of the network) and large-scale (input of the network) patches, respectively, where $X=I/Q/U$ defines the Stokes parameters, $y^{\circ}$ the physical dimension of the patch in degrees, and $z\arcmin$ its angular resolution in arcminutes.

The GAN architecture, training procedure, input maps, and output validation of the ForSE model are fully described in KP2021.

Figure 1: Structures in ForSE (KP2021) and ForSE+ (this work). ForSE is designed to achieve a $12\arcmin$ resolution from input large scales at $80\arcmin$. We have three kinds of newly trained models here: ForSE+S12 and ForSE+S3 to produce stochastic maps at $12\arcmin$ and $3\arcmin$, and ForSE+D to generate a deterministic map at $3\arcmin$. The output, $M_{12\arcmin}^{Q/U,20^{\circ}}$, from ForSE and ForSE+S12 are the input to ForSE+D and ForSE+S3 to get deterministic and stochastic small scales at $3\arcmin$, respectively.

2.2 ForSE+: Producing non-Gaussian dust maps at $3\arcmin$

In this work, we implement ForSE+, an updated version of the ForSE code and its training procedure, with two objectives:

1. To allow the generation of full-sky polarization maps with non-Gaussian features at an enhanced resolution of $3\arcmin$.

2. Starting from the same low-resolution maps, to generate multiple realizations of stochastic small scales that still have the correct non-Gaussian statistical properties.

In the sections below, we explain the assumptions and methodology used to achieve these two goals. Figure 1 sketches the input and output of the new ForSE+ models and in Table 1 we summarize the new models trained in this paper.

Model	Description (Output)	Input
ForSE	KP2021 version to simulate deterministic thermal dust emission maps at $12\arcmin$	Planck GNILC maps at $80\arcmin$
ForSE+S12	new model to simulate stochastic thermal dust emission at $12\arcmin$	Planck GNILC maps at $80\arcmin$ plus random component
ForSE+D	new model to simulate deterministic thermal dust emission maps at $3\arcmin$	Maps at $20\arcmin$, smoothed from ForSE $12\arcmin$ maps
ForSE+S3	new model to simulate stochastic thermal dust emission maps at $3\arcmin$	Maps at $20\arcmin$ smoothed from ForSE+S12 plus random component

Table 1: Summary of three newly trained models in this paper and the first version of the model proposed in KP2021.

In order to incorporate these extensions, we used the same GAN architecture as for the ForSE model, re-implemented using the new Tensorflow framework (version 2.6.0), and we performed additional fine-tuning steps of some hyper-parameters of the networks to improve our results.

2.2.1 Scale invariance assumption

The first goal of our implementation of ForSE+ is to generate a map of polarized dust emission with non-Gaussian features at arcminute angular scales. As we anticipated, this is crucial to understand the possible impact of non-Gaussianity in the extraction of the CMB $B$-mode lensing signals, which peak at these scales. In order to achieve this, we need to find a way to overcome the current limitation in the observational data at our disposal. As a matter of fact, in order to train our GANs we can only rely on a training set composed of 350 patches, with dimensions of $20^{\circ}\times 20^{\circ}$ and $320\times 320$ pixels, at an angular resolution of $12\arcmin$, taken from the total intensity GNILC Planck map at 353 GHz, described in KP2021. No other observations of thermal dust emission at a higher angular resolution in a portion of the sky large enough to be used as training set are available.

The idea to circumvent this restriction and still be able to reach a resolution higher than $12\arcmin$, even in the absence of a proper training set, is to make a scale-invariance assumption, applying our GAN model in an iterative way. In practice, our set of training squared patches with dimensions of $20^{\circ}\times 20^{\circ}$ and a resolution of $12\arcmin$ can equally be treated as having dimensions of $5^{\circ}\times 5^{\circ}$ and resolution of $3\arcmin$, since the network does not have a sense of physical units and is only sensitive to the ratio between the dimension of the patch and its angular resolution (i.e., $\tilde{m}^{I,20^{\circ}}_{12^{\prime}}\equiv\tilde{m}^{I,5^{\circ}}_{3^{\prime}}$). In this way, the same dataset used to train the first version of ForSE in KP2021 can be used to train a new GAN model. This model takes the output of the first iteration of the code as input and generates non-Gaussian scales at $3\arcmin$. As has been mentioned, this implies that we are assuming scale invariance for the statistical properties of dust emission; that is, scales at $12\arcmin$ have the same properties as those at $3\arcmin$. This assumption is justified by the fact that current observations of the dust polarization power spectra shows that they can, at the first order, be approximated as a power law as a function of the angular scales (Planck Collaboration XI, 2020).

Additional pre- and post-processing of the input patches (including upsampling, smoothing, and the sub-division of patches) is needed to train the GAN in the correct way and to be able to restore full-sky maps, as is described in Section 3.2.1.

2.2.2 Stochasticity

Our second goal is to produce different realizations of the non-Gaussian small-scale structures. This is important to estimate the variance of the signal we are simulating as well as the correlation among different angular scales. The way we achieved this was by simply adding a random component to the large-scale maps that are the input of our GANs. We then trained new models on these signal + random component maps, always using the 350 $\tilde{m}_{12^{\prime}}^{I,20^{\circ}}$ patches as the training set.

The random component that we considered was simply generated as a random realization from a uniform distribution in the range $[-1,1]$. Since our input maps were always re-scaled to have pixel values ranging in the same interval to be compatible with the input of our neural networks, we had a signal-to-noise ratio ($\rm SNR)$ of $\sim 1$. We refer to this as “stochastic training;” in other words, a random component was added to the input signal, ForSE+S, as opposed to the deterministic case (ForSE+D), in which we did not add any random component to the input maps.

3 Pre-processing, training, and post-processing of ForSE+

In this section, we describe the training details of ForSE+S and ForSE+D, including the pre-processing and post-processing steps, and present the results on flat-squared patches, including maps and power spectra, before reprojecting them to full-sky maps.

3.1 ForSE+S to $12\arcmin$

We first describe how we generate stochastic maps with non-Gaussian features at a resolution of $12\arcmin$. We also applied a similar procedure to construct maps at $3\arcmin$, as is described in Section 3.2. For clarity, we will call ForSE+S12 the model that generates stochastic maps at $12\arcmin$ and ForSE+S3 the one that goes up to $3\arcmin$.

As was mentioned above, we injected stochasticity into our generative process by simply adding a random component to the low-resolution maps that are the input of our GAN. By doing so, ForSE+S is able to generate non-Gaussian small-scale features that still partially depend on the real observed large-scale structures but will have a different morphology as we change the realization of the random component in the input dataset.

3.1.1 Training and post-processing

As in KP2021, we trained two models for $Q$ and $U$ maps separately. The inputs to the generator, $G$, were the 174 $M^{Q/U,20^{\circ}}_{80\arcmin}$ patches that together cover the full sky, with dust signal plus an additional random component. The training set, which encodes the target statistical distribution of small scales, was the 350 $\tilde{m}^{I,20^{\circ}}_{12\arcmin}$ maps. The weights of the neural works were updated for $2\times 10^{5}$ epochs and saved every 500 steps. Since we do not want to generalize the usage of the trained model – the trained model only needs to predict the output from the training data – it is not a problem if there is an over-fitting during training. Therefore, there was no consideration of a separate validation set during the training process.

Following KP2021, we used the Minkowski functionals (MFs), as implemented in (Mantz et al., 2008), as the criteria with which to choose the best epoch. For two-dimensional fields, three kinds of MFs can be built: $\mathcal{V}_{0}$, $\mathcal{V}_{1}$, and $\mathcal{V}_{2}$. They fully describe the statistical properties of the field and represent the covered area ($\mathcal{V}_{0}$), the boundary length ($\mathcal{V}_{1}$), and the number difference of connected domains and holes of the image’s feature ($\mathcal{V}_{2}$), as a function of the threshold, $\rho$ (Hadwiger, 1957).

During the GAN training process, the goal is to produce small-scale feature maps with statistical properties as close as possible to the ones of the training set. We quantified the level of agreement by calculating the overlapping fractions between the distributions of MFs of the maps in the training set ($\tilde{m}^{I,20^{\circ}}_{12\arcmin}$) and those of the generated ones ($\tilde{m}^{Q/U,20^{\circ}}_{12\arcmin}$). The distributions of MFs are indicated by the $\pm 1\sigma$ variation among the training set or generated maps and the overlapping fractions were computed as the ratio between the intersection area and the total area spanned by the two distributions. In practice, we calculated the MFs overlapping for each saved epoch of $G$ and selected as the best model the one with the highest score. In doing so, we computed the MFs for the output maps, $\tilde{m}^{Q/U,20^{\circ}}_{12\arcmin}$, by applying $G$ to maps with the realization of a random component different from the one used for training.

The best models are obtained after 5500 epochs and 6000 epochs for $Q$ and $U$, respectively. Their MFs are shown in Figure 2 and compared with the ones from the training set. The overlap among the distributions is at a level of $50\%-60\%$, comparable with the one obtained in Krachmalnicoff & Puglisi (2023) (that includes corrections to KP2021). In comparison, Gaussian small scales have MFs with obviously different shapes from these two sets of maps, as is illustrated in Fig.7 of KP2021.

Figure 2: MFs of small scales at $12\arcmin$ produced by ForSE+S, compared with the ones from the intensity maps in the training set. The overlapping fractions (OFs) of each pair of MFs are also shown. The dashed line represents the mean over the set of different patches and the shaded area is their standard deviation. The distribution of $Q$ maps is shown in the upper panel, and $U$ in the lower panel.

We also computed the overlapping fraction of the MFs by considering 100 different realizations of the small-scale maps (obtained by changing the random component in the input). The mean (standard variation) values of the overlapping fraction for $\mathcal{V}_{0}$, $\mathcal{V}_{1}$, and $\mathcal{V}_{2}$ are 59.1%(1.5%), 62.9%(0.6%), and 55.8%(0.7%), respectively, for Stokes $Q$ maps, and 58.8%(2%), 56.3% (1.3%), and 45.5%(1.1%), respectively, for $U$ maps. These numbers show the robustness of ForSE+S in generating stochastic small scales with non-Gaussian high-order statistics.

Since in the training procedure both the input maps and the training sets are normalized in the range of [-1,1], the output maps also have pixel values in this range, and therefore need to be normalized to restore physical units. We achieved this by following the procedure of KP2021, hence ensuring that, for each patch and for both $Q$ and $U$, the amplitude of the power spectrum of the produced small scales matches the extrapolation at higher multipoles of the power spectrum of the low-resolution input maps at $80\arcmin$.

3.1.2 Results of ForSE+S at $12\arcmin$

Figure 3 shows $M_{12\arcmin}^{Q/U,20^{\circ}}$ patches with two different realizations of small scales at $12\arcmin$ from ForSE+S12, after the normalization mentioned above, in the second and third columns, compared with the deterministic results from ForSE in the first column. Both the differences at small scales and consistency at large scales between the outputs from ForSE and ForSE+S12 show the capability of the trained model to produce small-scale features in the map space.

Figure 3: Maps of $20^{\circ}\times 20^{\circ}$ patches at $12\arcmin$. From left to right are the deterministic map from ForSE and two stochastic realizations from ForSE+S12. Throughout this paper, all the maps are shown in Galactic coordinates.

Figure 4: EE and BB power spectra of the squared patches at $12\arcmin$ shown in Figure 3. Dash-dotted black lines are the power spectra for GNILC $80\arcmin$ patches and green lines are for the deterministic $12\arcmin$ patches. Gray lines show the power spectra of all the 100 realizations from ForSE+S12 and red lines are the means.

To further validate these maps, we calculated the second-order statistics – the power spectrum – using the Namaster package (Alonso et al., 2019). In Figure 4, we show the EE and BB power spectra from the $QU$ squared patches of low-resolution maps, the output from ForSE, and 100 realizations of ForSE+S12 in black, green, and gray, respectively. The mean values of power spectra from 100 stochastic realizations are also shown in red. The output maps from ForSE+S12 are consistent with the deterministic output maps in terms of the power spectra, as was expected.

3.2 ForSE+D and ForSE+S3 to $3\arcmin$

We describe now the procedures that we followed to generate maps at the resolution of $3\arcmin$, by using our GAN model iteratively both in the case of ForSE+D and for ForSE+S3. The whole procedures, including several pre- and post-processing steps, are sketched in Figure 5 and described in the following (see also Foschi (2021)).

Figure 5: Pipeline to prepare the training data for ForSE+D (upper row, from left to right) and post-processing to transform the output patches from ForSE+D into the full-sky map at $3\arcmin$ in the HEALPix format (lower row, from right to left). 174 in the upper right corner of each tile means the number of images. The first number beneath the images indicates the number of pixels on each side of the image and the second expression has the same subscript and superscript as that described in Section 2.1. We note that after the “Dividing” and before the “Composing” steps, each original image was divided into 49 sub-patches, with a side length of $5^{\circ}$. At the end, all the flat patches were reprojected onto a HEALPix map with $N_{side}$ equal to 4096.

3.2.1 Pre-processing for the training of ForSE+D

As was mentioned above, we reached the resolution of $3\arcmin$ by assuming scale invariance in the statistics of the thermal dust emission. We exploited the output of the first GAN model as the input to a second generative step by using the same training set and considering it to be composed of 350 patches with dimensions of $5^{\circ}\times 5^{\circ}$ at a resolution of $3^{\prime}$, as was explained in Section 2.2.1. Therefore, since the first iteration of ForSE allows one to go from maps with dimensions of $20^{\circ}\times 20^{\circ}$ and a resolution of $80^{\prime}$ to maps at a resolution of $12^{\prime}$, the second iteration can produce $5^{\circ}\times 5^{\circ}$ maps at a resolution of $3^{\prime}$. In order to preserve the proportions among all the relevant quantities (i.e., patch dimension, resolution of input, and resolution of output), the input patches for the second iteration should have dimensions of $5^{\circ}\times 5^{\circ}$ and a resolution of $20^{\prime}$. We can obtain those patches by smoothing and subdividing the output of ForSE.

In practice, we pre-processed each of the 174 $M^{Q/U,20^{\circ}}_{12\arcmin}$ maps with $320\times 320$ pixels obtained from the first iteration in the following way:

1. We upsampled each map with $320\times 320$ pixels to $1280\times 1280$ by repeating each pixel four times.

2. We smoothed these maps with a Gaussian kernel in order to reach a resolution of $20^{\prime}$.

3. We subdivided each of these large squared patches into $5^{\circ}\times 5^{\circ}$ patches, each one containing $320\times 320$ pixels. We divided each large patch into 49 small ones, with a large overlap among them, made of 160 pixels on each side. This overlap is needed in order to avoid border effects when the composition and reprojection on the sphere is performed.

4. At the end of this procedure, we have a set of $174\times 49=8526$ patches for $Q$ and $U$, which is the total amount of patches covering the full sky and represents the $M^{Q/U,5^{\circ}}_{20\arcmin}$ that will be used as the input to the second GAN iteration.

We applied the same pre-processing to the output of ForSE and ForSE+S12 in order to produce maps at $3^{\prime}$ in both the deterministic case and the stochastic one. In the stochastic case, we added an additional random component, as is described in Section 2.2.2.

3.2.2 Training and post-processing of ForSE+D

Once the pre-processing steps described above are performed, the obtained 8526 patches can be used to train a new GAN model. However, the time cost is basically linear with the number of patches. We note that the 8526 patches have repeating pixels among them (see the steps above to get 8526 patches); thus, they are actually not independent. In order to save computational time, we fed as input to the generator, $G$, only a subset of 696 patches, randomly selected from the total 8526 patches. Once the network was trained, we applied it to the remaining patches.

In Figure 6 we show the MFs of the generated small scales at $3\arcmin$. The overlap with the target distribution is at a level of 70-80% in the deterministic case. In the stochastic case (not shown), on the other hand, the overlap ranges between 60% and 70%, showing therefore that we are also able to generate non-Gaussian small-scale features at this higher resolution.

Figure 6: MFs overlapping between deterministic 3$\arcmin$ small scales and intensity small scales.

The output of the GAN models are patches, $\tilde{m}^{Q/U,5^{\circ}}_{3\arcmin}$, containing the small-scale features that, as in the first iteration, have pixel values in the range of $[-1,1]$. We therefore normalized them in physical units before multiplying them with the large-scale maps at the resolution of $20^{\prime}$, obtaining 8526 $M^{Q/U,5^{\circ}}_{3\arcmin}$ maps.

We recall that each of these patches comes from a subdivision of the larger maps with physical dimensions of $20^{\circ}\times 20^{\circ}$ (as is described in Section 3.2.1) into 49 sub-patches with a large overlap among each other. Therefore, before reprojecting them on the sphere to obtain the full-sky maps, we needed to recombine them. In order to avoid border effects, we did this by using $\cos^{2}$ apodization window functions, as is shown in Figure 7: each sub-patch was weighted with the corresponding window function, then they were summed together to form the 174 $M^{Q/U,20^{\circ}}_{3\arcmin}$ maps that could then be reprojected on the sphere.

Figure 7: Apodized window function of each sub-patch at different positions to compose 49 sub-patches with a side length of $5^{\circ}$ into a patch with a length of $20^{\circ}$.

3.2.3 Results of ForSE+D and ForSE+S3 at $3\arcmin$

Figure 8: $20^{\circ}\times 20^{\circ}$ patches of three different resolutions. From left to right are GNILC maps at $80\arcmin$, ForSE maps at $12\arcmin$, and ForSE+D maps at $3\arcmin$.

Figure 9: EE and BB power spectra of the squared patches at $3\arcmin$ shown in Figure 8. Dash-dotted black lines are the power spectra for GNILC $80\arcmin$ patches and purple lines are for the deterministic $12\arcmin$ patches. Green lines depict the behavior of ForSE+D maps, while gray lines show the power spectra of all the 100 realizations from ForSE+S3 and red lines are the means.

Figure 10: Upper panel: Two stochastic realizations of $Q$ maps, $\tilde{m}_{3\arcmin}^{Q,5^{\circ}}$, from ForSE+S3. Lower panel: $M_{3\arcmin}^{Q,5^{\circ}}$ maps of the upper panel after the normalization step. This patch is at the position of the dashed red box in the upper right plot of Figure 8.

In Figure 8, we show the comparisons of a selected patch at $80\arcmin$, $12\arcmin$ from ForSE and maps at $3\arcmin$ from ForSE+D, from left to right. It is clear that ForSE+D can inject small-scale structures, following the modulation of the large-scale emission.

The power spectra, computed on the same patch, are shown in Figure 9. As can be seen, the amplitude at low $\ell$ matches the one from the low-resolution GNILC maps, and the generated small scales have power at higher multipoles, with no breaks in the power spectrum that follows a power law, as was expected.

ForSE+S3 was utilized to generate stochastic small scales at $3\arcmin$, as is mentioned in 2.2.2. Two realizations of small scales for the $Q$ map in the range of [-1, 1] with a side length of $5^{\circ}$ and centered on $[288^{\circ},-61.5^{\circ}]$ (i.e., at the position of the dashed red box in the upper right plot of Figure 8) are shown in the upper panel of Figure 10, and in the lower panel we show the normalized maps combined with the large scales. The differences between the two realizations should be noted and both of them trace the large-scale features of the input maps. We also show the power spectra of 100 realizations in Figure 9, where we also plot the mean, validating the results of stochasticity on the power spectrum level.

4 Validations of full-sky maps from ForSE+

In this section, we show the maps and power spectra after reprojecting the flat patches back to the sphere, following the algorithm in the appendix of KP2021. Before showing the results, we introduce a further step to adjust the E-to-B ratio in the simulated maps to match the observations at large scales.

4.1 Fine-tuning of the E-to-B ratio to match observations

Figure 11: E-to-B ratios of non-tuned and tuned maps in green and red lines, for different sky masks.

We recall that, due to the limitation of observational data in polarization, we used the statistical properties of the total intensity small scales to be the ground truth of both $Q$ and $U$ during the training process. This implies that the output high-resolution map will have the same power for $E$ and $B$ modes. However, high-frequency observations of the Planck satellite have shown the existence of an asymmetry between the measured power of thermal dust emission in the $E$ and $B$ modes, with $A_{BB}/A_{EE}\sim 0.5$ over the multipole range of $40\leq\ell\leq 600$ (Planck Collaboration XI, 2020). Therefore, in order to match the real observations, we applied a fine-tuning to the full-sky maps, obtained after all the steps mentioned above.

We first transformed $QU$ maps into the harmonic space to obtain the $a_{\ell m}^{E}$ and $a_{\ell m}^{B}$ coefficients. Then we applied a factor of $\sqrt{0.5}$ to the $a_{\ell m}^{B}$ to decrease its power, since $a_{\ell m}^{E}$ and $a_{\ell m}^{B}$ have the same variance as was just mentioned. Finally, we transformed the tuned $a_{\ell m}s$ back to the map space to get maps that match the observed ratio of 0.5 at the power spectra level. We note that the large scales in the output maps of ForSE/ForSE+ remain to be the observed ones, so the tuning was applied only to the small scales; that is, $a_{\ell m}s$ belonging to $\ell>135$ and the transition between large and small scales was smoothed with a sigmoid function to avoid discontinuities in the final power spectra.

The E-to-B ratios of the power spectra calculated for different fractions of sky are shown in Figure 11. The red lines show the ratios after the tuning steps, which were indeed corrected to $\sim 2$ for the injected small scales.

4.2 Full-sky maps

We are now ready to present the full-sky maps and specifically, we show the ForSE+D maps at $3\arcmin$, at Nside = 4096, in Figure 12, compared with the input GNILC ones at $80\arcmin$. The difference between the two highlights the presence of small-scale features in the ForSE+D maps. When zooming into a patch that consists of two $M^{Q/U,20^{\circ}}_{3\arcmin}$ patches we also find no clear edge effects. These results show that border effects from reprojection are negligible.

In Figure 13, the power spectra of the maps on the sphere at different resolutions are shown for different sky fractions. As was expected, our maps at $3\arcmin$ can further extrapolate the power up to $\ell\sim 3600$, which corresponds to the drop scale of $3\arcmin$. We also note the absence of discontinuities at the transition from large scales to small scales ($\ell>200$), implying that the small scales are produced with a similar pattern as the one at large scales, attributed to the normalization step to rescale the generated small scales in the range of [-1, 1] to the correct amplitude, as was mentioned in Section 3.2.2. In order to make a comparison with the latest d9 dust model from PySM3 (PanEx Collaboration, 2024), we also show the power spectra of d9 maps at 353GHz in blue. Although the latter was produced with entirely different methods, the power spectra of ForSE+D maps at $3\arcmin$ are impressively close to those of d9 maps, even for different sky fractions.

Figure 12: Top panel: Full-sky polarization maps (left: Stokes $Q$, right: Stokes $U$) for the GNILC template at a $80\arcmin$ angular resolution. These maps are the input to our algorithm. Middle panel: Maps with small-scale features, up to $3\arcmin$, generated by ForSE+D. Bottom panel: The difference between the two maps. The residuals mostly encode smaller angular scales, as was expected.

Figure 13: EE and BB Power spectra of GNILC $80\arcmin$, deterministic $3\arcmin$ in gray and red, with $f_{sky}$ = 1, 0.8, and 0.4 in solid, dash-dotted, and dashed lines, respectively. In order to make a comparison, we also show the power spectra of PySM3 d9 dust map at 353GHz in blue.

We also calculated the power spectra of 100 realizations of full-sky maps at $3\arcmin$, generated with ForSE+S3, with a $80\%$ sky mask, up to $\ell_{max}=4096$ and with a bin width of $\Delta\ell=160$. The covariance matrix from these 100 realizations is shown in Figure 14. The correlation among multipoles at small scales is further evidence that the small scales at $\ell>800$ were synthesized in a non-Gaussian way. In fact, if the small-scale features were produced with Gaussian properties, we would not observe any non-diagonal correlation, which is verified from our experiment to calculate the covariance matrix of a purely Gaussian field. We devote the next section to further deepening the non-Gaussian properties of our maps.

Figure 14: Covariance matrices of maps at $3\arcmin$ on the sphere with an $80\%$ sky mask, shown in absolute values. Calculated for power spectra of an $\ell$-bin range of [40, 4096] with a bin width of 160. Non-diagonal terms are normalized with the diagonal values, i.e., $R_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}C_{jj}}}$.

4.3 Non-Gaussianity measured on the sphere

In previous sections, we considered the MFs in order to characterize non-Gaussianity for flat patches. Here, we expand the analysis to the full-sky maps, by exploiting the algorithms described in Grewal et al. (2022), in order to calculate MFs for spherical maps in the HEALPix format. We focus on the small scales injected by ForSE and ForSE+D and here we used the multipole range [200, 2048] to band-pass filter the raw maps and applied a $80\%$ sky mask to ignore the Galactic plane.

In Figure 15 we show the results of the MFs statistics applied to the sphere, confirming our previous claims. Differences in the ForSE+D map are visible with respect to the case of Gaussian maps in dashed gray for all three kinds of MFs. The Gaussian maps were generated from a random realization from the power spectra of the ForSE+D $Q$ map. We also show the results of the latest PySM3 d9 dust maps in blue, which exhibits a deviation from Gaussianity, although being similar to the modulated Gaussian maps in orange, whose non-Gaussianity is supposed to derive only from the modulation of large scales. We further checked that the results are robust for the $40\%$ sky mask and for $\ell_{min}$ = 500 and 1000.

Figure 15: The three kinds of MFs for the ForSE+D $Q$ full-sky map at $3\arcmin$ with an $80\%$ sky mask (red). The raw map is band-pass filtered to keep small scales only, corresponding to the multipole range [200, 2048]. The dashed gray lines show the MFs of the Gaussian maps. Results for PySM3 d9 and modulated Gaussian maps are also shown in blue and orange, respectively.

5 Discussion

In this paper, we extend the ability of the ForSE model proposed in KP2021 based on the GAN technique to simulate non-Gaussian stochastic small scales of polarized thermal dust emission up to $3\arcmin$. We have trained three new models, which are summarized in Table 1, together with the one trained in KP2021.

Based on the results obtained in KP2021, our first test was to bring stochasticity into our models by adding a random component into the input and to train a new network called ForSE+S12 so that it can generate different maps with different seeds. Minkowski functionals were used to quantify the level of non-Gaussianity in the maps and in Figure 2 we demonstrate that the polarized thermal dust small scales have a similar level of non-Gaussianity to that in the intensity small scales. Different realizations of maps at $12\arcmin$ for a specific patch and the corresponding power spectra are shown in Figure 4, which indeed have expected variations at small scales and the correct amplitude scaling with multipoles for power spectra.

We then considered the case of a $3\arcmin$ angular scale, where lensing $B$-modes have the strongest signal, and we studied the deterministic case first, ForSE+D. By relying on the scale invariance assumption discussed in Section 2.2.1 and the pre- and post-processing steps outlined in Figure 5, we trained the model to generate maps at $3\arcmin$ out of those at $12\arcmin$, which are the output from ForSE. MFs distributions in Figure 6 represent a verification of the training process. Maps and power spectra of a selected patch of maps at $3\arcmin$ are shown in Figure 8, where the $3\arcmin$ power is present as expected. ForSE+S3 was finally trained in order to generate stochastic small scales at $3\arcmin$. Two realizations of small scales are shown in Figure 10, with their power spectra shown in the right panel of Figure 8, showing the expected multipole scaling.

After obtaining the small scales on flat patches, we reprojected them onto the celestial sphere in order to get full-sky maps in the HEALPix format, with $N_{side}=2048$ for maps at $12\arcmin$ and 4096 for maps at $3\arcmin$. We show the Mollview projection of the deterministic maps at $3\arcmin$ in Figure 12 and compare them with the observed Planck GNILC maps, validating both the effectiveness of the injected small scales and the reprojection process. The power spectra of deterministic maps at $80\arcmin,12\arcmin$, and $3\arcmin$ with $100\%,80\%$, and $40\%$ sky masks are plotted in Figure 13, indicating that the small scales generated preserve the same anisotropies as the low-resolution observations.

We further obtain the covariance matrix of power spectra up to $\ell_{max}\sim 4000$ from 100 realizations of maps at $3\arcmin$ in Figure 14, which has a strong correlation between different multipoles at small scales, and thus highlights the non-Gaussianity in our maps. We note that for a Gaussian field, the variances along the diagonal line are smaller than the variances we obtained here and the off-diagonal correlation is zero. We repeated the calculation of the covariance matrix for the PySM3 d11 model, which is also designed to generate multiple realizations of small scales of thermal dust emission (PanEx Collaboration, 2024), and find that off-diagonal elements are also close to zero, meaning that the level of non-Gaussianity of the injected small scales on the full sky is small.

Finally, we measured the level of non-Gaussianity using MFs in the simulated full-sky maps in Figure 15, and compared it with the small scales simulated within the Gaussian assumption. The difference between the MFs obtained from the two kinds of maps is another clear indication of the non-Gaussian small-scale component generated with ForSE+.

Once the steps above are achieved, we can use the synthetic dust Stokes $Q$ and $U$ maps at 353GHz at $3\arcmin$ resolution as a template and appropriately scale across different frequencies by taking the observed dust SEDs. The series of maps can then serve as a simulation suite to obtain the scatter of component separation against variations in the foreground realization. Moreover, high-resolution maps are essential in order to validate the lensing reconstruction methods, tested so far on foreground models where the arcminute structure is either absent or Gaussian (Lonappan et al., 2023), or based on models of the Galactic magnetic fields in order to reproduce a non-Gaussian pattern (Beck et al., 2020). In a forthcoming work, we will investigate if any variation of the results shown in Beck et al. (2020) occurs when considering the foreground maps exhibiting non-Gaussianity generated through GANs in this paper.

With the new version of PySM3 under development, we capitalize on integrating our maps into a new model of the latest PySM packages for public use. Maps will be shared upon request to the authors. The code to generate maps has also been made publicly available. Finally, we report here future improvements in the ForSE algorithm. First, more observations are needed. In fact, the assumptions made in Section 2.2 are somewhat a compromise due to the lack of observation at the required resolution. The models will get more reliable as more observations become available. Second, by considering the network architecture, the loss function accounting for the non-Gaussianity of the targets produced by the GAN may be considered, as in the current implementation this process is not automatized. A way to quantify the level of non-Gaussianity with a formalism that is differentiable with respect to the input pixels would be desirable, as then we could construct the loss function of the network in order to include the information of non-Gaussianity, which will effectively guide the generator. The last point we want to mention is that the operation of adding noise is to some extent like the training process of diffusion models in deep learning, which is more natural than what is done in this work, so if we turn to use the network of diffusion model, it may have a better performance.