A Constrained NILC method for CMB B mode observations

The objective of the ILC method is to accurately isolate the specific emission of interest, such as CMB, from a set of multi-frequency maps. This technique assumes that the emission of interest remains consistent across different frequencies, allowing for the construction of a single cleaned map that predominantly represents the desired signal.

Suppose that all the frequency maps are calibrated to CMB units and the same resolution, the observation can be modeled as:

$\displaystyle\mathbf{y}(p)=\mathbf{1}s(p)+\mathbf{f}(p)+\mathbf{n}(p).$

(2.1)

Here, $\mathbf{y}(p)$, $\mathbf{1}$ and $\mathbf{n}(p)$ are $n_{\nu}$ dimensional vector, where $n_{\nu}$ denotes the number of frequency channels, $p$ can be any domain like pixel, harmonic, needlet, etc. The terms $\mathbf{f}$ and $\mathbf{n}$ describe the contribution from galactic diffuse foreground contamination and instrumental noise respectively. To minimize the impact of the “junk” term $\mathbf{f}$+$\mathbf{n}$, the ILC provides an estimator $\hat{s}_{\rm ILC}$ of $s$ by forming the linear combination $\hat{s}_{\rm ILC}(p)=\mathbf{w}^{T}\mathbf{y}(p)$, which has minimum variance and also preserves the CMB signal with $\mathbf{w}^{T}\mathbf{1}=1$. Such the weights are given by [10]

$\displaystyle\mathbf{w}=\frac{\hat{\mathbf{C}}^{-1}\mathbf{1}}{\mathbf{1}^{T}\hat{\mathbf{C}}^{-1}\mathbf{1}},$

(2.2)

where $\hat{\mathbf{C}}$ is the empirical covariance matrix of the observations.

Suppose the number of frequency channel is larger than the combined number of foreground components and the CMB signal, and in an ideal scenario with no noise, the ILC estimator is unbiased. However, in the real world, bias always exists and the bias of ILC may come from different situations as follows:

• •

  The non-stationary of foreground and noise within the selected $p$ domain introduces a challenge as it leads to varying weights, which can introduce bias in the ILC method. To address this issue, various approaches have been explored to ensure the stationary of the field $\mathbf{f}$ and noise $\mathbf{n}$. One approach is to carefully choose appropriate $p$ domains by 1) decomposing the sky maps into several regions[10], 2) harmonic domain[12], 3) needlet domain[31]. Additionally, it is beneficial to scan the sky as uniformly as possible during observation.

• •

  The residual of “junk” term can be defined as $r=\mathbf{1}^{T}\hat{\mathbf{C}}^{-1}(\mathbf{f}+\mathbf{n})/(\mathbf{1}^{T}\hat{\mathbf{C}}^{-1}\mathbf{1})$. When $\mathbf{f}$ and $\mathbf{n}$ are random field with the zero means, $\hat{s}_{\rm ILC}$ is unbiased, but when there exits non-zero higher order terms of auto-term and cross-term of $\mathbf{r}$, for example $\left<r^{2}\right>+2\left<sr\right>$ will introduce bias[13]:

  a) The auto-term bias can the written as following with the assumption that there is no correlation between foreground and noise:

  $\displaystyle\left<r^{2}\right>=\frac{\mathbf{1}^{T}\hat{\mathbf{C}}^{-1}(\left<\mathbf{f}\mathbf{f}^{T}\right>+\left<\mathbf{n}\mathbf{n}^{T}\right>)\hat{\mathbf{C}}^{-1}\mathbf{1}}{(\mathbf{1}^{T}\hat{\mathbf{C}}^{-1}\mathbf{1})^{2}},$

  (2.3)

  both the bias from foreground and noise increase with the noise level. However, we can easily remove noise bias by employing noise-only simulations with the knowledge of statistical property of noise which is described in Section 4.2. Unfortunately, the same approach cannot be applied to the foreground. For low SNR $B$-mode physics, such bias is typically non-negligible. It is possible to obtain an unbiased ILC estimator by replacing the $\mathbf{C}$ with $\mathbf{C}-\mathbf{N}$, where $\mathbf{N}$ is the noise covariance, but this will leads to a substantial increase in the variance of CMB power[32].

  b) The cross-term bias mainly arise from the difference between the true underlying covariance matrix $\mathbf{C}$ and the estimated empirical covariance matrix $\hat{\mathbf{C}}$. The bias values approximated amounts to $2(1-n_{\nu})\sigma^{2}_{s}/N_{p}$, which means that the ILC estimator underestimates the CMB power by roughly $2(1-n_{\nu})C_{\ell}/N_{p}$, more details refer to [13].

In this paper, we focus on reducing the ILC bias with the case of low SNR. To minimize this bias, more additional information are needed to distinguish the Galactic foreground signal from the instrumental noise. In other words, extra constraints, which depend on the additional information, can be imposed on the weights.

Suppose the emission law of all the Galactic foreground and the CMB signal are known, thus the observation can be modeled as

$\displaystyle\mathbf{y}(p)=\mathbf{As}(p)+\mathbf{n}(p).$

(2.4)

Here we also assume that all maps have been re-smoothed to the same resolution. The mixing matrix $\mathbf{A}$ are known from the prior. $\mathbf{s}(p)$ is a column vector with $N_{\rm comp}$ elements denote the templates for all the components (CMB + foreground). And thus the term $\mathbf{n}(p)$ includes the instrumental noise only. The goal of constrained ILC is to reduce the impact of residual foreground signal while preserving the CMB signal. Thus the constraint becomes $\mathbf{A}^{T}\mathbf{w}=\mathbf{e}$ with $e_{\rm CMB}=1$, $e_{\rm non-CMB}=0$. The weights for minimizing the variance of the resulting map with the constraints can be solved using Lagrange multipliers[11]. The solution is

$\displaystyle\mathbf{w}=\hat{\mathbf{C}}^{-1}\mathbf{A}(\mathbf{A}^{T}\hat{\mathbf{C}}^{-1}\mathbf{A})^{-1}\mathbf{e}.$

(2.5)

The foreground bias in this issue will arise bias of $\mathbf{A}$ in two main aspects: misunderstanding of the number of foreground components and the bias associated with the individual elements of $\mathbf{A}$.

From our previous work, we have shown that the NILC on $B$ maps is very promising on $r$ constraint, so here we choose to estimate the weights in needlet domain. The following will describe the constrained NILC in detail, the multi-frequency data is in the form of needlet coefficients which given by

$\displaystyle\bm{\beta}_{jk}=\sum_{\ell m}\mathbf{y}_{\ell m}h_{\ell}^{(j)}Y_{\ell m}\left(\hat{\xi}_{k}^{(j)}\right)$

(2.6)

where $\mathbf{y}_{\ell m}$ is the spherical harmonic coefficients of the observational $\mathbf{y}(p)$ which can be arbitrary scalar or pseudo-scalar map. $\hat{\xi}^{(j)}_{k}$ is a set of grid points on sphere. In practice, they are exactly the center of the pixel $k$ in Healpix pixelization scheme with different $\texttt{NSIDE}(j)$. $h_{\ell}^{(j)}$ is the spectral windows satisfying $\sum_{j}(h_{\ell}^{(j)})^{2}=1$. The NILC estimation on needlet space is

$\displaystyle\hat{\beta}^{\rm NILC}_{jk}=\mathbf{w}^{T}_{jk}\bm{\beta}_{jk}.$

(2.7)

The final estimation $\hat{s}_{\rm NILC}$ in real space is

$\displaystyle\hat{s}_{\rm NILC}(p)=\sum_{j}\sum_{\ell m}\hat{\beta}_{j,\ell m}^{\rm NILC}h_{\ell}^{(j)}Y_{\ell m}(p)$

(2.8)

The weights in (2.7) are calculated from the equation (2.5) with an empirical covariance $[\hat{\mathbf{C}}]_{jk}$ which are estimated by the formula

$\displaystyle\hat{\mathbf{C}}_{jk}=\frac{1}{N_{\mathcal{D}}}\sum_{k^{\prime}}\omega_{j}(k,k^{\prime})\bm{\beta}_{jk^{\prime}}\bm{\beta}_{jk^{\prime}}^{T}.$

(2.9)

It means that the covariance matrix obtained by computing the mean of the needlet coefficients. In practice, the weights in (2.9) is a 2D Gaussian function. Due to the locality of the NILC method, it is convenient to extend equation (2.5) by using different mixing matrices $\mathbf{A}$ to compute weights at different positions. The estimation of $\mathbf{A}$ matrices from the multi-frequency maps will be shown in section 4.

We adopted a sky model that is as close as possible to the real observation. It contains CMB signal, diffuse foreground emission composed of synchrotron, thermal dust, free-free and spinning dust. The polarized foreground maps are simulated by using the Planck Sky Model[33] (PSM). Here are the details of the sky model:

1. 1.

   CMB: the input CMB maps are Gaussian realizations from a particular power spectrum obtained from the Boltzmann code CLASS²²2http://lesgourg.github.io/class_public/class.html[34], using Planck 2018 best-fit cosmological parameters [35] with the tensor scalar ratio $r=0$. The CMB dipole was not considered in the simulations.

   The effect of CMB distorted by weak gravitational lensing was also considered in our simulation. Weak lensing distorts a part of $E$-mode signal to $B$-mode and it is not a linear effect and therefore it is non-Gaussian. The precise way to simulate such lensing distortion effects is to rearrange the pixel on the map level. For this part, the package LensPix³³3https://cosmologist.info/lenspix/ is used to simulate it. the algorithm is presented in the reference[36], which distorts the primordial signal given a realization of lensing potential map from $C_{\ell}^{\phi\phi}$.

2. 2.

   Synchrotron: Galactic synchrotron radiation is produced by charged relativistic cosmic rays that are accelerated by the Galactic magnetic field. At frequencies below 80 GHz, synchrotron radiation is the primary polarized foreground. In our simulations, we did not account for the effects of Faraday rotation, and instead modeled the synchrotron emission across frequencies using template maps based on specific emission laws. To simplify matters, we modeled the frequency dependence of synchrotron radiation in the CMB frequency range as a power-law function (in Rayleigh-Jeans units):

   $\displaystyle\begin{bmatrix}T_{\rm sync}(\nu)\\ Q_{\rm sync}(\nu)\\ U_{\rm sync}(\nu)\end{bmatrix}\propto\begin{bmatrix}T_{template}\\ Q_{template}\\ U_{template}\end{bmatrix}\times\nu^{\beta_{s}}.$

   (3.1)

   The spectral index $\beta_{s}$ typically has a value of $-3$. When generating temperature maps, we utilized the 408 MHz radio continuum full-sky map of Haslam et al. [37, 38] as our template. For polarization maps, we used the SMICA component separation of Planck PR3 (Planck Public Data Release 3) [39] polarization maps as templates, with the same synchrotron spectral index.

3. 3.

   Thermal dust: At higher frequencies (above 80 GHz), the thermal emission from heated dust grains becomes the dominant foreground signal. In the PSM, we modeled the emission of thermal dust as modified black-bodies [40],

   $\displaystyle I_{\nu}\propto\nu^{\beta_{\rm d}}B_{\nu}(T_{\rm d}),$

   (3.2)

   where $T_{\rm d}$ is the temperature of the dust grains, $A_{\rm d}$ is the amplitude and $\beta_{\rm d}$ is the spectral index. $B_{\nu}(T_{\rm d})$ is the Planck black-body function at temperature $T_{\rm d}$ at frequency $\nu$. For polarization, the Stokes $Q$ and $U$ parameters (in ${\rm MJy}/{\rm sr}$ unit) are correlate with the intensity,

   $\displaystyle Q_{\nu}(\hat{\mathbf{n}})$

   $\displaystyle=f(\hat{\mathbf{n}})I_{\nu}\cos(2\gamma_{d}(\hat{\mathbf{n}})),$

   $\displaystyle U_{\nu}(\hat{\mathbf{n}})$

   $\displaystyle=f(\hat{\mathbf{n}})I_{\nu}\sin(2\gamma_{d}(\hat{\mathbf{n}})),$

   (3.3)

   where $f(\hat{\mathbf{n}})$ and $\gamma_{\rm d}(\hat{\mathbf{n}})$ denote polarization fraction and polarization angle at different position $\hat{\mathbf{n}}$. We use the dust template in Planck PR3 at $353$ GHz in intensity and polarization as our template, which includes the template maps of the dust grains $T_{\rm d}$, the amplitude for all $T,Q,U$ at 353 GHz, and the spectral indices.

4. 4.

   free-free: Free-free emission is produced by electron-ion scattering in interstellar plasma and is generally fainter than synchrotron or thermal dust emission. It is modeled as a power-law with a spectral index $\beta_{\rm f}$ of approximately $-2.14$ [41] and is intrinsically unpolarized. Free-free emission only contributes to temperature measurements, and we used the free-free emission map from the Commander method [16] in Planck 2015 data PR2 [41] as our template.

5. 5.

   Spinning dust: Small rotating dust grains can produce microwave emission if they have a non-zero electric dipole moment, which can also be modeled as a power-law within the CMB observation frequency bands. We used the template spinning dust map from Commander in Planck PR2 as our input template. Spinning dust emission is weakly polarized, and in our simulation, the polarization fraction of spinning dust was set to 0.005.

In this work, we consider a ground-based observational experiment that can only cover a partial sky area. The scanning strategy is optimized for detecting primordial gravitational waves on Northern hemisphere, in other words, the scanning strategy was designed to perform a deep scan of a small sky region. The hit-map corresponding to the scanning strategy is shown in Figure 1. The deep scanning sky patch is centered at RA $10$h$04$m, Dec. $61.49^{\circ}$. The effective area of the sky-patch is about $7.5\%$ of the full-sky.

Considering the atmospheric window, we have set up two observation frequency bands of 95 and 150 GHz for the ground-based experiments. In order to achieve better foreground removal results for ILC, we also simulated some space-based observed sky maps and taking them as the inputs. The maps including the K and V bands of WMAP, as well as all bands of Planck High Frequency Instrument (HFI). The beam and noise information of them was obtained from the LAMBDA website⁴⁴4https://lambda.gsfc.nasa.gov.

The frequency bands and the related beam resolution are listed in Table 1. All the maps are in the HealPix pixelization scheme at $\texttt{NSIDE}=512$. The noise in the sky patch is not homogeneous. For simulated sky maps of ground-based observations, the noise variance on each pixel depends on the hit-map, which shown in Figure 1. The variance is given by

$\displaystyle\sigma^{2}_{p}=\frac{{\rm NET}^{2}}{H_{p}\Omega_{p}},$

(3.4)

where $H_{p}$ denotes the value of hit-map at pixel $p$ (in ${\rm hour}/{\rm deg^{2}}$ unit) and $\Omega_{p}$ is the solid angle of pixel $p$. In the case of simulating space-based observed sky maps, the noise variance map can be obtained from the LAMBDA website.

To perform a constrained NILC on the simulated observational sky maps, there is still a key issue to be solved, which is how to obtain the mixing matrix $\mathbf{A}$ in equation (2.5). Generally speaking, the term $\mathbf{A}$ encodes the emission law of the foreground. Based on the current understanding of the foreground model, only synchrotron radiation and thermal dust emission can contribute polarization signal. From the synchrotron and dust model equations (3.1) and (3.2), it appears that we need to know these three parameters $\beta_{\rm s},\beta_{\rm d}$ and $T_{\rm d}$ in advance.

In our implementation, the parameter $T_{\rm d}$ is fixed to be $19.6$ K according to the results of Planck[43]. As for $\beta_{\rm s}$ and $\beta_{\rm d}$, We will fit them from the observational intensity and polarization sky maps. The fitting is based on three assumptions: 1. The foreground emission is dominated by large-scale signals. 2. The spectral indices of any foreground components modeled by power-law change very little on small scales. 3. The frequency dependence of foreground emission is the same for both intensity and polarization signals. We fit the values of these two parameters at each pixel using the $8\ (frequency\ channel)\times 3\ (T,Q,U)$ observed sky map. The model at each pixel $p$ is

	$\displaystyle T_{\nu}(p)$	$\displaystyle=A_{\rm s}(p)\frac{g_{\rm RJ}(\nu)}{g_{\rm RJ}(\nu_{\rm s})}\left(\frac{\nu}{\nu_{\rm s}}\right)^{\beta_{\rm s}}+A_{\rm f}(p)\frac{g_{\rm RJ}(\nu)}{g_{\rm RJ}(\nu_{\rm s})}\left(\frac{\nu}{\nu_{\rm s}}\right)^{-2.14}+A_{\rm sd}(p)\frac{g_{\rm Jy}(\nu)}{g_{\rm Jy}(\nu_{\rm s})}\frac{f_{\rm sd}(\nu)}{f_{\rm sd}(\nu_{\rm s})}$
		$\displaystyle+A_{\rm d}(p)\frac{g_{\rm Jy}(\nu)}{g_{\rm Jy}(\nu_{\rm d})}\frac{B_{\nu}(T_{\rm d})}{B_{\nu_{\rm d}}(T_{\rm d})}\left(\frac{\nu}{\nu_{\rm d}}\right)^{\beta_{\rm d}(p)}+T_{\rm CMB}(p),$		(4.1)
	$\displaystyle Q_{\nu}(p)$	$\displaystyle=Q_{\rm s}(p)\frac{g_{\rm RJ}(\nu)}{g_{\rm RJ}(\nu_{\rm s})}\left(\frac{\nu}{\nu_{\rm s}}\right)^{\beta_{\rm s}}+Q_{\rm d}(p)\frac{g_{\rm Jy}(\nu)}{g_{\rm Jy}(\nu_{\rm d})}\frac{B(\nu,T_{\rm d})}{B(\nu_{\rm d},T_{\rm d})}\left(\frac{\nu}{\nu_{\rm d}}\right)^{\beta_{\rm d}(p)}+Q_{\rm CMB}(p),$		(4.2)
	$\displaystyle U_{\nu}(p)$	$\displaystyle=U_{\rm s}(p)\frac{g_{\rm RJ}(\nu)}{g_{\rm RJ}(\nu_{\rm s})}\left(\frac{\nu}{\nu_{\rm s}}\right)^{\beta_{\rm s}}+U_{\rm d}(p)\frac{g_{\rm Jy}(\nu)}{g_{\rm Jy}(\nu_{\rm d})}\frac{B(\nu,T_{\rm d})}{B(\nu_{\rm d},T_{\rm d})}\left(\frac{\nu}{\nu_{\rm d}}\right)^{\beta_{\rm d}(p)}+U_{\rm CMB}(p).$		(4.3)

Here $\nu_{\rm s}\equiv 23$ GHz and $\nu_{\rm d}\equiv 353$ GHz. The factors $g_{\rm RJ}(\nu)$ and $g_{\rm Jy}(\nu)$ denote the conversion factors from antenna temperature $\mu\rm{K}_{\rm RJ}$ or flux density Jy/sr to CMB unit $\mu{\rm K}_{\rm CMB}$ at frequency $\nu$. We considered the spinning dust emission are unpolarized in our fitting model. $f_{sd}$ is the external template calculated from SpDust2 [44]. Totally there are $12$ free parameters in the model include the amplitudes and the spectral indices of the foreground, and the CMB signal.

In order to achieve better fitting results, all the observational maps are re-smoothed to the same resolution $52.8$ arcmin and down graded the NSIDE to $8$. This step can average out the noise and CMB signal on each pixel, thereby improving the fitting results. After smoothing and downgrading, the spectral index for each pixel is the average of the true spectral index values within the pixel area. It can be easily figure out by perform Taylor expansion

$$s_{\rm FG}(p)=\int{\rm d}\bm{n}A(\bm{n})\nu^{\beta(\bm{n})}\simeq\bar{A}\nu^{\bar{\beta}}+\bar{A}(\ln\nu)^{2}\int(\beta(\bm{n})-\bar{\beta})^{2}+\dots,$$

(4.4)

with $\bar{A}\equiv\int{\rm d}\bm{n}W(\bm{n})A(\bm{n})$ and $\bar{\beta}\equiv\int{\rm d}\bm{n}W(\bm{n})\beta(\bm{n})$, where $W(\bm{n})$ is the weight that encodes the information for smoothing and downgrading. Based on our aforementioned assumptions, we only consider the zeroth-order term.

The likelihood for each pixel, denoted as $\mathcal{L}(p)$, is given by

$\displaystyle-2\log\mathcal{L}(p)=\sum_{\nu}\frac{\left[\hat{T}_{\nu}(p)-T_{\nu}(p)\right]^{2}}{\sigma^{2}_{T}(\nu,p)}+\frac{\left[\hat{Q}_{\nu}(p)-Q_{\nu}(p)\right]^{2}}{\sigma^{2}_{P}(\nu,p)}+\frac{\left[\hat{U}_{\nu}(p)-U_{\nu}(p)\right]^{2}}{\sigma^{2}_{P}(\nu,p)}+const.$

(4.5)

Here $\sigma^{2}_{T}$ and $\sigma^{2}_{P}$ represent the variances of instrumental noise of the temperature map and polarization map, respectively. The variances are recalculated at $\texttt{NSIDE}=8$. For each pixel, we sampled the likelihood (4.5) by using Markov Chain Monte Carlo (MCMC). For efficiency purpose, we use the python package NumPyro⁵⁵5https://github.com/pyro-ppl/numpyro[45, 46] to sample the likelihood. When running MCMC, the priors for all parameters are uniform distributed. For $\beta_{\rm s}$ and $\beta_{\rm d}$, we use the ranges $\beta_{\rm s}\in[-4,-2]$ and $\beta_{\rm d}\in[1,3]$. For every pixel $p$, we take the best-fit value sampled from MCMC as the fitting results $\beta_{\rm s}(p)$ and $\beta_{\rm d}(p)$. Thus, $\beta_{\rm s}$ and $\beta_{\rm d}$ are the maps with the resolution $\texttt{NSIDE}=8$ after the fitting. An fitting example is given in Figure 2 shows the fitting results of some of the parameters.

According to the current understanding of the polarization foregrounds[41, 47], the model for polarized signal includes CMB, synchrotron and thermal dust emission. This suggests that we can apply two other separate constraints on the NILC method. One for removing synchrotron emission signal $\sum_{\nu}w_{\nu}A_{\nu,{\rm s}}=0$ where $A_{\nu,{\rm s}}$, which is the element of mixing matrix $\mathbf{A}$ and given by Equation (3.1) and another for removing thermal dust emission signal $\sum_{\nu}w_{\nu}A_{\nu,{\rm d}}=0$ where $A_{\nu,{\rm d}}$ is given by Equation (3.2). It should be noted that the conversion factors to CMB unit are needed when calculating $A_{\nu,{\rm s}}$ and $A_{\nu,{\rm d}}$. Due to the scarcity of low-frequency observation data, fitting the model (4.1 – 4.3) independently on each pixel results in a substantial error in the estimation of parameter $\beta_{\rm s}$ for all pixels. In this case, the 95% CL for $\beta_{\rm s}$ is about $0.5$. Therefore, in order to mitigate the error, we impose the constraint that parameter $\beta_{\rm s}$ remains constant across all pixels during the fitting process. However, it is important to note that parameter $\beta_{\rm d}$ is allowed to vary across different pixels, as it represents distinct parameters for each pixel. In this case, the signal-to-noise ratio for the parameter $\beta_{\rm s}$ is improved, reaching a value of $0.1$ at a 95% CL. The posterior distribution for $\beta_{\rm s}$ and $\beta_{\rm d}$ are shown in Figure 2.

To better recover the $B$-mode CMB signal, the constrained weights $\mathbf{w}_{jk}$ in Equation (2.5) are calculated on $B$ maps defined as

$\displaystyle B(p)=\sum_{\ell m}a_{\ell m}^{B}Y_{\ell m}(p).$

(4.6)

The $EB$ leakage[23] due to limited sky coverage is also corrected by a so-called template subtraction method[48] when obtain $B$ maps from observational $QU$ maps. We have provided a detailed description of the steps for the EB leakage correction in the appendix. The correction method has been proven to be the optimal blind correction on the pixel domain in the literature[49]. We have outlined the specific correction steps and results in Appendix A. We show that the ILC methods performs slightly better on such leakage-free $B$-family maps[25]. In addition, because of the localization of NILC, the mixing matrix $\mathbf{A}$ can vary with spatial position as $\mathbf{A}_{jk}$. In order to match the NSIDE of $\beta_{\rm d}(p)$ maps with $\texttt{NSIDE}(j)$ for every needlet band $j$, we use map2alm and alm2map routines in healpy⁶⁶6https://github.com/healpy/healpy[50, 51] package to up-grade the resolution.

In order to more intuitively describe how the constrained NILC process is executed from beginning to end, Figure 3 shows the flow chart of the entire data processing process.

After obtaining the clean $B$ map through the constrained NILC method, the post-mapmaking analysis procedures are quite straightforward. In general, the subsequent processing involves the following steps.

Angular power spectra estimation: The pseudo-power spectra are estimated by using the python package of NaMaster⁷⁷7https://github.com/LSSTDESC/NaMaster[52]. In order to avoid instability caused by signal jumps at the edges of the sky regions, a mask that was slightly smaller than the mask shown in Figure 1LABEL:sub@subfig:skypatch by about $2$ degrees is used. The power spectra are binned every $30$ multipole moments from $\ell=35$ to $215$.

Noise debias: We use the $50$ noise only simulations to go through the same constrained NILC processing while keeping the weights unchanged to reconstruct $50$ reconstructed noise only maps, then estimate the mean value the band-powers $N_{\ell}$ over these noise simulations, finally subtract it from the band-power calculated from the resulting map of ILC. After subtraction, the band-power is considered to be noise-unbiased. Also, the $50$ noise only simulation band-powers are used to estimate the noise covariance matrix.

Constraint on $r$: To show how good our method is, we also use the calculated unbiased clean spectra to constrain the tensor to scalar ratio $r$. The full posterior for the individual band powers is non-Gaussian. However, for high enough multipoles (usually $\ell\geq 30$) the central limit theorem justifies a Gaussian approximation[53]. The lowest multipole $\ell$ in our binning scheme is $35$ thus the Gaussian approximation is valid:

$\displaystyle-2\ln\mathcal{L}=(\hat{C}_{\ell_{b}}^{BB}-C_{\ell_{b}}^{BB})^{T}{\rm Cov}(\hat{C}_{\ell_{b}}^{BB},\hat{C}_{\ell_{b}}^{BB})^{-1}(\hat{C}_{\ell_{b}}^{BB}-C_{\ell_{b}}^{BB})+const.$

(4.7)

where the observed power spectrum is denoted as $\hat{C}_{\ell_{b}}^{BB}$. $C_{\ell_{b}}^{BB}$ denotes the theoretical power spectrum which was calculated by CAMB code. The likelihood only concerns the $BB$ spectra at $\ell\geq 30$, where the primordial gravitational wave and lensing effect matters, so we fixed all the parameters to the best fit value from Planck-2018[54], while only free $r$ with prior $r\in[0,1]$. The covariance term is estimated from the noise only simulations. We sample this likelihood for these two parameters using the CosmoMC[55]⁸⁸8https://cosmologist.info/cosmomc/ package. The posteriors are summarized by GetDist[56]⁹⁹9https://pypi.org/project/GetDist/ package.

The foreground residual maps can be obtained by applying the same weights obtained from constrained NILC or standard NILC to the input foreground maps. The resulting and foreground residual $B$-mode maps are shown in Figure 4. When comparing the results of the constrained NILC and standard NILC methods at the map level, it’s easy to find that qualitatively, the resulting map of the constrained NILC exhibits slightly larger fluctuations compared to the standard NILC. However, this trend is reversed when considering the foreground residuals. More quantifiable results can be observed at the power spectrum level. It’s clear that the foreground residual of the constrained NILC method is much smaller. From the power spectra analysis, the bias caused by foreground residual in the constrained NILC method is estimated to be around $r\sim 10^{-2}$. On the other hand, the standard NILC method exhibits a significantly larger foreground residual, estimated to be approximately at the order of $r=10^{-1}$.

One notable observation from the residual foreground maps is that the residual foreground of standard NILC exhibits a higher degree of small-scale structure compared to the residual of constrained NILC. This indicates that the foreground and the instrumental noise are indistinguishable especially on small scale where the noise dominate. The root of the problem lies in the fact that standard NILC lacks any prior information regarding noise and foreground components , as mentioned in Section 2.1. Adding prior information helps to get results with lower foreground residuals. But if the foreground and noise are considered as a whole “junk” term, then the overall residual will become larger. This is why the fluctuations in the resulting map of constrained NILC are greater than that of standard NILC. In a word, adding the extra constraints on NILC can help reduce the bias, but slightly increases uncertainty.

One of the most important scientific goal of observing CMB $B$-mode polarization is to find out the polarized signal produced by the primordial gravitational waves. Figure 5 illustrates the distribution of the best-fit values (left) and the differences between the $95\%$ upper limits and the best-fit values (right) for both the standard NILC and constrained NILC methods, based on $500$ independent simulations. It is evident that the constrained NILC method has a significantly smaller bias, less than $5\times 10^{-3}$, compared to the standard NILC method, which exhibits a bias of approximately $r\sim 0.05$. Additionally, as anticipated, the figure also demonstrates that the difference, which represents the level of uncertainty, increases accordingly. Indeed, for the standard NILC method, the $2\sigma$ uncertainty is approximately $0.025$, while for the constrained NILC method, it increases to around $0.045$. This discrepancy arises because the weights used in the constrained NILC method are optimized specifically for minimizing the foreground residual, rather than minimizing the overall impact of foregrounds and noise. As a result, there is a trade-off between reducing bias and increasing uncertainty in practice. In conclusion, the obtained fitting results align with the analysis of residual foregrounds conducted in the previous section.