© The Authors 2023

1. Introduction

Observations of the temperature and polarisation anisotropies of the cosmic microwave background (CMB) have led to the establishment of a cosmological concordance scenario, the Λ cold dark matter (ΛCDM) model, whose parameters are now tightly constrained (see, e.g. Boomerang: MacTavish et al. 2006; WMAP: Hinshaw et al. 2013; Planck: Planck Collaboration VI 2020). CMB polarisation data are usually decomposed into the so-called E and B modes of even and odd parity in the sky, respectively (Kamionkowski et al. 1997; Zaldarriaga & Seljak 1997).

E modes, which are mostly generated on the last scattering surface by scalar perturbations, have been clearly detected by a number of experiments and contribute to confirm the ΛCDM cosmological scenario (see, e.g. DASI: Kovac et al. 2002; WMAP: Spergel et al. 2003; Boomerang: Montroy et al. 2006; Planck: Planck Collaboration VI 2020). B modes, instead, are thought to be associated with at least two independent mechanisms. On small angular scales (typically a few arcminutes), they are mainly sourced by primordial E modes distorted into B modes via weak gravitational lensing by the intervening large-scale structures along the line of sight from the last scattering surface to us. This signal is commonly known as lensing B modes (Zaldarriaga & Seljak 1998) and has been observed by several experiments, although the detection still has a modest signal-to-noise ratio or is made in very small patches of the sky (Hanson et al. 2013; Ade et al. 2014; Planck Collaboration Int. XLI 2016; Ade et al. 2016; Keisler et al. 2015; POLARBEAR Collaboration et al. 2017).

On large scales, B modes are expected to originate from tensor perturbations (primordial gravitational waves) generated in the very early Universe by a phase of cosmic inflation (Kamionkowski et al. 1997). Their magnitude is set by the tensor-to-scalar ratio r, which is defined as the amplitude of primordial gravitational waves over the one of initial density perturbations. Such a primordial B-mode signal still escapes detection so far with an upper limit of r ≤ 0.032 at 95% CL (Tristram et al. 2022). The detection of the primordial gravitational waves’ background would represent a powerful test of cosmic inflation and would be a means to discriminate among numerous theoretical models. The main features in the BB tensor angular power spectrum are the reionisation bump at very large angular scales (ℓ ≲ 10), which is associated with the integral of the linear polarisation generated by quadrupolar anisotropies from the last scattering surface as seen by each free electron after reionisation, and the recombination bump on smaller scales (ℓ ∼ 80), which represents the imprint of the primordial gravitational waves on the physics of the last scattering surface.

Many experiments have been designed or proposed to observe B-mode polarisation, either from the ground, such as for POLARBEAR (Arnold et al. 2010), QUBIC (Qubic & Battistelli 2011), BICEP (Wu et al. 2016), Keck-Array (Keating et al. 2003), LSPE-STRIP (Bersanelli et al. 2012), ACT (Aiola et al. 2020), SPT (Sayre et al. 2020), and Simons Observatory (Ade et al. 2019); from balloons for SPIDER (Filippini et al. 2010) and LSPE-SWIPE (de Bernardis et al. 2012); or from space for LiteBIRD (Matsumura et al. 2014) and PICO (Hanany et al. 2019). The development of effective component separation methods has become one of the crucial aspects of data analysis for all of these experiments. The reason is that the quest for primordial B modes is made much more difficult by the presence of instrumental noise and polarised foregrounds, especially Galactic thermal dust and synchrotron emission.

These methods are usually divided into blind, parametric, and template-fitting techniques. Blind methods, such as ILC (Bennett et al. 2003; Tegmark et al. 2003), NILC (Basak & Delabrouille 2012), FastICA (Maino et al. 2002), and SMICA (Delabrouille et al. 2003), usually linearly combine multi-frequency maps in such a way as to minimise some meaningful statistical quantity of the final map. Parametric methods, such as Commander (Eriksen et al. 2008), FGBuster (Stompor et al. 2009), and FastMEM (Stolyarov et al. 2002), instead, explicitly model the frequency properties of the foregrounds by means of a set of parameters that are fitted to the data. Such methods provide an easy way to characterise and propagate foregrounds’ residual errors, but their effectiveness depends on how reliable the adopted model is. Finally, template-fitting algorithms, such as SEVEM (Martínez-González et al. 2003), try to construct internal foreground templates either in pixel or harmonic space using multi-frequency observations.

There is no clear evidence (especially in polarisation) as to which approach is more effective in reconstructing the CMB signal. Thus, applying several different methods on the same dataset allows for one to perform comparisons and evaluate the robustness of the results.

In this work, we focus our attention on the so-called Internal Linear Combination (ILC) algorithms. They are of great relevance in the context of CMB data analysis because they perform noise and foreground subtraction with minimal prior information on their properties and, hence, they are not prone to systematic errors due to mis-modelling of the spectral energy distribution of the Galactic emission components. The ILC methods will then play a key role in the analysis of future CMB experiments, given our still limited knowledge of the properties of the polarised foregrounds as of yet. One of the drawbacks of such approaches is that the estimation of foreground residuals usually relies on Monte Carlo simulations or other bootstrapping techniques. Throughout the years, several extensions have been proposed; the Needlet Internal Linear Combination (NILC), which performs variance minimisation in needlet space, has proven to be one of the most promising blind techniques. The NILC method has been extensively applied to full-sky satellite data; however, in the near future, new observations of the CMB polarisation will be obtained mainly from ground-based and balloon-borne experiments, which will only observe a portion of the sky. Therefore, in this work, we explore the possibility of employing NILC for the analysis of partial-sky CMB B-mode data, specifically addressing the further challenges that such an extension yields: E–B leakage, needlet filtering, and beam convolution. The importance of accurately characterising the impact of these operations has also recently been pointed out by Zhang et al. (2022). In our analysis, we consider two complementary experiments: (i) the Short Wavelength Instrument for the Polarisation Explorer (SWIPE), which is a balloon-borne telescope that avoids most of the atmospheric contamination, whose scientific goal is the observation of both the reionisation and recombination peaks; and (ii) the Simons Observatory (SO), which is a ground-based experiment targeting the detection of the recombination bump.

The SWIPE telescope (de Bernardis et al. 2012)) is one of the two instruments of the Large Scale Polarisation Explorer (LSPE) experiment (Addamo et al. 2021). It is a balloon-based array of bolometric polarimeters that will survey the sky in three frequency bands centred at 145,  210, and 240 GHz. It is expected to constrain the tensor-to-scalar ratio down to r = 0.015 (at 95% CL). The 145 GHz band will be the main channel for CMB science, while observations at 210 and 240 GHz will monitor thermal dust contamination. The SWIPE telescope is scheduled to be launched from the Svalbard islands and will observe a large sky fraction (around 35%) with an angular resolution of full width at half maximum (FWHM) equal to 85′ for around 15 days, exploiting the optimal observational conditions offered by the Arctic night (at latitude around 78°N).

Simons Observatory consists of a large aperture telescope (LAT) with a 6-m primary mirror and three 0.5-m refracting small aperture telescopes (SATs; Ade et al. 2019). They will be located in the Atacama Desert at an altitude of 5.200 m in Chile’s Parque Astronomico and will collect data in the frequency range between 25 and 280 GHz. It is designed to measure a signal at the r = 0.01 level at a few σ significance or to exclude it at similar significance, using the B-mode amplitude around the recombination bump. The reionisation peak, instead, cannot be observed by SO due to atmospheric contamination that will only allow for cosmological modes with ℓ > 30 to be sampled.

The paper is organised as follows. In Sect. 2, we present the considered simulated multi-frequency datasets. In Sect. 3 we describe the ILC methods considered in this work and how to extend the application of the NILC algorithm to partial-sky B-mode data. In Sect. 4 you can find the procedures adopted to assess the performance of ILC methods. In Sect. 5 we provide the obtained results from the application of NILC to LSPE-SWIPE and SO simulated data. Finally, we report our conclusions in Sect. 6.

2. Simulated datasets

The simulated datasets employed in this analysis include the two aforementioned CMB experiments: LSPE-SWIPE and the SO-SAT. We do not consider SO-LAT, because it will not be devoted to primordial B-mode science. In the case of LSPE-SWIPE, realistic simulated Planck maps have been included in the component separation pipeline to have a broader frequency coverage to better trace the Spectral Energy Distributions (SEDs) of B-mode foregrounds.

The main properties of each dataset are listed in Table 1 (see Addamo et al. 2021; Planck Collaboration I 2020; Ade et al. 2019). All maps are generated by adopting the HEALPix pixelisation scheme (Górski et al. 2005) with N_side = 128, which corresponds to a pixel resolution of 27.5′.

The Q and U maps at the different frequencies are obtained from the co-addition of three separate components: CMB, Galactic emission (synchrotron and thermal dust), and Gaussian white and isotropic noise. We do not include polarised extra-Galactic sources, since they are expected to be dominant on angular scales smaller than those of main interest for inflationary B-mode science (ℓ > 100; Puglisi et al. 2018).

The CMB maps are generated with the HEALPix Python package¹ as random Gaussian realisations from the angular power spectra of the best-fit Planck-2018 parameters (Planck Collaboration VI 2020) with tensor-to-scalar ratio r = 0, unless otherwise specified. Instrumental noise is simulated as white and isotropic Gaussian realisations with standard deviations in each pixel associated with the polarisation sensitivities listed in Table 1.

Galactic emission is generated using the PySM Python package² (Thorne et al. 2017). Synchrotron polarised emission is modelled with a simple power law (Rybicki & Lightman 1979):

$$\begin{array}{c}\hfill X_{\mathrm{s}}(\widehat{n},\nu )=X_{\mathrm{s}}(\widehat{n},{\nu }_0)·(\frac{\nu }{{\nu }_0}{)}^{{\beta }_{\mathrm{s}}(\widehat{n})}\end{array}$$(1)

with X = {Q, U}, $\widehat{n}$ the position in the sky and ν the frequency. $X_{\mathrm{s}}(\widehat{n},{\nu }_0)$ represents the synchrotron template at a reference frequency ν₀. Thermal dust emission is modelled with a modified black-body (MBB):

$$\begin{array}{c}\hfill X_d(\widehat{n},\nu )=X_d(\widehat{n},{\nu }_0)·(\frac{\nu }{{\nu }_0}{)}^{{\beta }_d(\widehat{n})}·\frac{B_{\nu }(T_d(\widehat{n}))}{B_{{\nu }_0}(T_d(\widehat{n}))},\end{array}$$(2)

where B_ν(T) is the black-body spectrum.

We consider two different Galactic models, commonly adopted within the forecast analyses of CMB polarisation experiments:

• d0s0, where β_s = −3,  β_d = 1.54, and T_d = 20 K are assumed to be constant across the sky;

• d1s1, where the spectral indices of synchrotron and thermal dust emission depend on the position in the sky.

In both cases, the dust template is the estimated dust emission at 353 GHz in polarisation from the Planck-2015 analysis (Planck Collaboration X 2016), smoothed with a Gaussian kernel of FWHM = 2.6° and with small scales added by the procedure described in Thorne et al. (2017). The synchrotron template is the WMAP 9-year 23-GHz Q/U map (Bennett et al. 2013), smoothed with a Gaussian kernel of FWHM = 5° and small scales added.

In the d1s1 model, dust temperature and spectral index maps are obtained from the Planck data using the Commander pipeline (Planck Collaboration X 2016). The synchrotron spectral index map is derived using a combination of Haslam 408 MHz observations and WMAP 23-GHz 7-year data (Miville-Deschênes et al. 2008). The d0s0 model is a simplified representation of Galactic polarised emission, because we know that the spectral properties of the foregrounds vary in the sky. However, it still represents an important starting point for assessing the performance of blind methods on polarisation B-mode datasets. On the other hand, the d1s1 is surely closer to what the polarised Galactic emission is expected to be.

The maps of all components are smoothed with the beam of the channel of the considered dataset with the lowest angular resolution and then added together. The sky patches observed by LSPE-SWIPE and SO-SAT are shown in Fig. 1. The SO-SAT footprint has a sky coverage of f_sky = 34%. However, due to the highly inhomogeneous scanning of the telescope, some regions of the sky will be much better observed than others, leading to an effective sky fraction of ∼10% (see Ade et al. 2019). Therefore, to perform a realistic foregrounds subtraction, the ILC methods are applied to the simulated SO-SAT dataset within the reduced patch shown in Fig. 2 and obtained considering only the 10% of the pixels with the highest values in the hit counts map provided by the SO Collaboration. Consistently, the sensitivity values reported in Table 1 refer to a homogeneous hit counts map with a sky fraction of f_sky = 10% (Ade et al. 2019).

Fig. 1. Sky coverage of the LSPE-SWIPE (left) and SO-SAT (right) instruments in Galactic coordinates. The observed region of the sky is shown in grey.

Fig. 2. Sky area with the highest signal-to-noise ratio in the SO-SAT patch. This was obtained considering pixels with the largest values in the map of hit counts for a final sky fraction of fsky = 10%. The ILC methods were applied within this patch.

3. Extension of ILC methods to partial-sky observations

The ILC methodology is one of the most widely used approaches for the reduction of contaminants in CMB observations. It was first adopted in the analysis of the intensity data from the WMAP satellite (Bennett et al. 2003). Then, throughout the years, several extensions of the basic algorithm have been proposed. The methods applicable to polarisation data considered in this work are:

• ILC,

• Polarisation ILC (PILC; Fernández-Cobos et al. 2016),

• Needlet ILC (NILC; Delabrouille et al. 2009).

All the above ILC techniques consist in linearly combining the N_ν multi-frequency maps of one or multiple CMB experiments with frequency-dependent weights ω:

$$\begin{array}{c}\hfill \widehat{X}={\displaystyle \sum _{j=1}^{N_{\nu }}}{\omega }_j·X_j,\end{array}$$(3)

where X is the B-field for ILC, P = Q + iU for PILC, and a set of B-mode needlet coefficients for NILC (as discussed in the text below).

The weights are derived by minimising the variance of $\widehat{X}$ to reduce the contamination of Galactic emission and instrumental noise. Therefore, among the component separation techniques, ILC methods require the smallest number of a priori assumptions. Recently, a further variation of these techniques has been proposed, the constrained Moments ILC (cMILC), where the weights estimation is constrained to de-project some moments of the foregrounds emission (Remazeilles et al. 2021). In this work, however, we do not consider such an extension.

Input maps of Eq. (3) have to be smoothed with a common beam that usually corresponds to the one of the channel with the lowest angular resolution. Each frequency map X_j can be described as the sum of the CMB signal s, foregrounds emission f_j, and instrumental noise n_j:

$$\begin{array}{c}\hfill X_j=a_{j}s+f_j+n_j,\end{array}$$(4)

where a_j are the calibration coefficients. Assuming that the observations are calibrated with respect to the component of interest (the CMB) and are expressed in thermodynamic units, we have a_j = 1,  ∀j. Therefore, to preserve the CMB signal in the output map, the weights must satisfy the further constraint $\sum _{j=1}^{N_{\nu }}{\omega }_j=1$. This last condition, together with Eqs. (3) and (4), enables us to demonstrate that the output solution contains the full CMB signal and some foregrounds and noise residuals:

$$\begin{array}{c}\hfill \widehat{X}={\displaystyle \sum _{j=1}^{N_{\nu }}{\omega }_j·(s+f_j+n_j)=s+\sum _{j=1}^{N_{\nu }}{\omega }_j·(f_j+n_j).}\end{array}$$(5)

In the case of simulated datasets, the estimation of foregrounds and noise residuals is straightforward, while for actual data they have to be estimated through Monte Carlo simulations based on our prior knowledge of the contamination in the input multi-frequency dataset.

The weights which minimise the output variance and whose sum is equal to unity can be analytically estimated for ILC, PILC, and NILC (see Bennett et al. 2003; Fernández-Cobos et al. 2016; Delabrouille et al. 2009). The only assumption of these methods is that the CMB has a known emission law (black-body spectrum) and no correlations with foregrounds or noise.

NILC, which is the main algorithm considered in this work, represents a refinement of the ILC method (Delabrouille et al. 2009), since the effectiveness of variance minimisation is enhanced when performed in the needlet domain. Needlets are a particular wavelet system that guarantees simultaneous localisation in harmonic and pixel space. They have been introduced into the statistical literature by Baldi et al. (2009) and have been applied for the first time to CMB data in Pietrobon et al. (2006).

In practise, the needlet coefficients of an input B-mode map at frequency i, ${\beta }_j^i$, are obtained by filtering its harmonic coefficients $a_{\ell m}^i$ with a weighting function b_j(ℓ) that selects modes at different angular scales for each needlet scale j:

$$\begin{array}{c}\hfill {\beta }_j^i(\widehat{\gamma })={\displaystyle \sum _{\ell ,m}(a_{\ell m}^i·b_j(\ell ))·Y_{\ell m}(\widehat{\gamma }),}\end{array}$$(6)

where $\widehat{\gamma }$ is a direction in the sky. This procedure in harmonic space is equivalent to performing a convolution of the map in real domain.

The shape of the needlet bands is defined by the choice of the harmonic function b whose width is set by a parameter B: lower values of B correspond to a tighter localisation in harmonic space (fewer multipoles entering into any needlet coefficient), whereas larger values result in wider harmonic bands. Commonly adopted constructions of the harmonic function b(ℓ) are the standard (Narcowich et al. 2006; Marinucci et al. 2008), the cosine (Basak & Delabrouille 2012), and the Mexican needlets (Geller & Mayeli 2008).

The input needlet maps are then linearly combined in such a way as to obtain a minimum variance map ${\beta }_j^{\mathrm{NILC}}$ at each scale j:

$$\begin{array}{c}\hfill {\beta }_j^{\mathrm{NILC}}(\widehat{\gamma })={\displaystyle \sum _{i=1}^{N_{\nu }}{\omega }_i^j(\widehat{\gamma })·{\beta }_j^i(\widehat{\gamma })=\sum _{\ell ,m}{a}_{\ell m,j}^{\mathrm{NILC}}·Y_{\ell m}(\widehat{\gamma }).}\end{array}$$(7)

The pixel-dependent weights are computed as follows:

$$\begin{array}{c}\hfill {\omega }_i^j(\widehat{\gamma })=\frac{{\sum }_k{C}_{\mathit{ik}}^j{(\widehat{\gamma })}^{-1}}{{\sum }_{\mathit{ik}}{C}_{\mathit{ik}}^j{(\widehat{\gamma })}^{-1}},\end{array}$$(8)

where the covariance matrix $C_{\mathit{ik}}^j(\widehat{\gamma })=⟨{\beta }_j^i·{\beta }_j^k⟩$ for scale j at pixel $\widehat{\gamma }$ is estimated as the average of the product of needlet coefficients in some space domain 𝒟. This domain is usually a Gaussian window function centred at $\widehat{\gamma }$, whose width varies with the considered needlet scale j.

The final NILC map is then reconstructed by filtering again the harmonic coefficients $a_{\ell m,j}^{\mathrm{NILC}}$ in Eq. (7) and summing them all for each ℓ and m:

$$\begin{array}{c}\hfill X_{\mathrm{NILC}}(\widehat{\gamma })={\displaystyle \sum _{\ell ,m}{a}_{\ell m}^{\mathrm{NILC}}·Y_{\ell m}(\widehat{\gamma })=\sum _{\ell ,m}(\sum _j{a}_{\ell m,j}^{\mathrm{NILC}}·b_j(\ell ))·Y_{\ell m}(\widehat{\gamma }).}\end{array}$$(9)

Minimising the variance (and hence the contamination) separately on different needlet scales leads to a more effective cleaning. On large scales, diffuse Galactic foregrounds dominate over the other components and are better removed in NILC with respect to ILC. The same is expected to happen for the instrumental noise on smaller scales (larger multipoles).

We note that any error on the CMB calibration in Eq. (4) could have a relevant impact on the CMB reconstruction, as discussed in Dick et al. (2010). Furthermore, as noted in Delabrouille & Cardoso (2009) and further discussed in Delabrouille et al. (2009), one of the main limitations of ILC methods is the generation of empirical correlations between the CMB and the contaminants (Delabrouille et al. 2009) induced by the departure of the empirical correlation matrix C_ik from its ensemble average due to the finite size of the sample over which it is estimated.

This can lead to a negative bias in the reconstructed CMB power spectrum, especially at low multipoles, due to the cancellation of CMB modes projected in the foregrounds and noise sub-space. This effect has to be taken into account when estimating the CMB power spectrum. The NILC bias problem in this analysis is further discussed in the Appendix A.

In the case of ground-based and balloon-borne experiments, ILC methods cannot be applied straightforwardly to B-mode data. Indeed, three main complications arise:

• the E–B leakage: some CMB E modes are estimated as B modes in the E–B decomposition performed on partial-sky Q and U maps (see Sect. 3.1). It affects the ILC and NILC pipelines;

• for NILC, needlet filtering of a partial-sky signal can lead to a loss of modes and power, especially close to the borders of the mask (see Sect. 3.2);

• degrading cut-sky polarisation maps to a common angular resolution, analogously to needlet filtering, can lead to errors in the reconstruction of the signal in the pixels close to the border of the observed patch (see Sect. 3.3). It impacts the application of ILC, PILC, and NILC.

The E–B leakage effect is associated with the fact that the E–B decomposition of partial-sky Stokes parameters is not exact: some modes (ambiguous modes) satisfy both the E- and B-conditions simultaneously. When we split the polarisation field into an E- and a B-part, these ambiguous modes can go into either component (Bunn et al. 2003). In the case of CMB, where E modes are much brighter than B modes, this leads to an over-estimation of the power in B modes, especially on large scales.

The needlet filtering, on the other hand, represents a convolution of the data. If B modes are observed only in a portion of the sky, some modes (on large scales) can be lost in the process if the configuration of the needlet bands is not carefully chosen, because the convolution mixes the signal of the pixels close to the border with the null one of the unobserved region.

Both these effects are relevant especially on those angular scales which are most sensitive to the amplitude of primordial tensor perturbations. Thus they must be properly analysed and treated to have a correct reconstruction of the primordial tensor CMB B-mode signal.

Finally, the ILC algorithms require input maps with a common angular resolution (which usually is the one of the frequency channel with the largest beam). The convolution for a beam, in analogy with the needlet filtering, is not trivial if one deals with observations only on a portion of the sky due to the leakage of null signal of the unobserved region within the patch.

The relevance and possible correction of such issues are explored in Sects. 3.1–3.3 by comparing the angular power spectra (evaluated in the footprint of the considered experiments) of leakage-corrected, needlet-filtered, and smoothed CMB simulations with those of the input exact B-mode signal, reconstructed with a full-sky B decomposition of input Q and U maps. In such analysis, we consider the entire footprint observed by SO, given in the right panel of Fig. 1, since this is the region over which the B modes will be reconstructed from Q and U maps. Specifically, we use the estimator of the pseudo-angular power spectrum, which indicates the rotationally invariant variance of the harmonic coefficients of a map:

$$\begin{array}{c}\hfill C_{\ell }=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{m=\ell }|a_{\ell m}{|}^2.}\end{array}$$(10)

It is customary to report an analogous quantity: $D_{\ell }=\frac{\ell (\ell +1)}{2\pi }{C}_{\ell }$, which would be constant on large scales for a scale-invariant primordial density perturbations spectrum. In this work, the estimation of the power spectra of the B-mode maps is performed with the NaMaster Python package³ that extends the estimator in Eq. (10) taking into account the effects of masking and beam convolution (see Alonso et al. 2019).

3.1. Correction of the E–B leakage

The polarisation field of the CMB is a spin-2 quantity that can be described in terms of the Stokes parameters Q and U as follows:

$$\begin{array}{c}\hfill P_{\pm }(\widehat{\gamma })=Q(\widehat{\gamma })\pm iU(\widehat{\gamma }),\end{array}$$(11)

where $\widehat{\gamma }$ denotes the position in the sky. This field can be expanded in spin-weighted spherical harmonics:

$$\begin{array}{c}\hfill P_{\pm }(\widehat{\gamma })={\displaystyle \sum _{\ell m}{a}_{\pm 2,\ell m}{}_{\pm 2}{Y}_{\ell m}(\widehat{\gamma }).}\end{array}$$(12)

From the harmonic coefficients in Eq. (12), it is possible to build the polarisation fields E (a scalar map of even-parity) and B (a pseudo-scalar field of odd-parity):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E_{\ell m}=& -\frac{1}{2}(a_{2,\ell m}+a_{-2,\ell m})\hfill \\ \hfill B_{\ell m}=& -\frac{1}{2i}(a_{2,\ell m}-a_{-2,\ell m}).\hfill \\ \hfill \end{array}\end{array}$$(13)

The construction of maps from E_ℓm and B_ℓm represents the E and B decomposition of a set of Q and U maps.

For partial-sky observations, we do not have information on Q and U on the entire sky and, therefore, Eqs. (11)–(13) lead to an incorrect reconstruction of E and B maps. The E–B decomposition on the cut-sky is indeed not unique due to the non-locality of the transformation and, therefore, some E modes will be interpreted as B modes and vice versa (ambiguous modes; see Lewis et al. 2001). The CMB E-mode power is much greater at all multipoles than that of B modes and therefore the issue of leakage of the E modes into the B modes is the most relevant for partial-sky observations, leading to an over-estimate of the power in the reconstructed CMB B-mode map.

Several different methods have been proposed to address the leakage correction problem (Lewis et al. 2001; Bunn et al. 2003; Kim & Naselsky 2010; Liu et al. 2019; Zhao & Baskaran 2010; Kim 2011; Ghosh et al. 2021). In this work, we consider:

• the recycling method introduced in Liu et al. (2019);

• the procedure presented in Zhao & Baskaran (2010), hereafter referred to as ZB method.

We test these techniques on 200 CMB-only simulations that include lensing from E modes and primordial tensor perturbations with r = 0.01 for the cases of the SWIPE and SO-SAT patches (shown in Fig. 1). We have chosen this value of the tensor-to-scalar ratio because it is close to the upper bound targeted at the 95% confidence level (CL) by LSPE in the case of no detection (Addamo et al. 2021) and represents the amplitude expected to be observed at high significance by SO (Ade et al. 2019).

The performance of these methods has been evaluated by comparing the angular power spectra, $D_{\ell }^{\mathit{BB}}$ and $D_{\ell ,in}^{\mathit{BB}}$, of leakage-corrected and exact B-mode maps. The exact B modes are reconstructed with a full-sky B decomposition of CMB Q and U simulations. Furthermore, we consider an effective tensor-to-scalar ratio:

$$\begin{array}{c}\hfill r_{\mathrm{leak}}=\frac{|D_{\ell ,\mathrm{in}}^{\mathit{BB}}-D_{\ell }^{\mathit{BB}}|}{D_{\ell ,r=1}^{\mathit{BB}}},\end{array}$$

which quantifies the absolute error we can make in the estimation of the amplitude of primordial tensor perturbations due to the presence of leakage residuals at the different angular scales. This parameter accounts for the error in the reconstruction of both the tensor modes and the lensing signal. We detail the recycling and ZB methods in the following sections.

3.1.1. The recycling method

The recycling method has been introduced in Liu et al. (2019). The procedure is summarised in Fig. 3. It consists in decomposing the masked polarisation field P = (Q, U) into the so-called E- and B-family: P_E = (Q_E, U_E) and P_B = (Q_B, U_B), which receive contributions only from E- and B-mode harmonic coefficients, respectively. For partial-sky observations, P_B is largely affected by the presence of ambiguous E modes due to the E–B leakage. To correct for this effect, in the recycling method, $P_B^{\prime }={(Q_B,U_B)}^{\prime }$ is constructed through the B decomposition of masked P_E and is then used as a template of the E–B leakage contamination. This template is linearly fitted to P_B and removed from it, thus providing a leakage-corrected B-family $P_B^{″}={(Q_B,U_B)}^{″}$. The fit is performed with a simple least-squares linear regression method. We then obtain the final CMB B-mode map through a B decomposition of $P_B^{″}$. The leakage correction of the recycling method is not exact, due to our lack of knowledge of the Q and U parameters in the unobserved region of the sky, and therefore some residuals are still present in the reconstructed map.

Fig. 3. Recycling method for E–B leakage correction.

The performance of the method on 200 CMB-only simulations with lensing+r = 0.01 for the SWIPE and SO-SAT patches is shown in the upper panels of Fig. 4, where the average angular power spectrum of the corrected B maps is compared to that of the exact signal reconstructed with a full-sky B decomposition of the simulated Q and U maps. As can be seen in the top left panel, if not corrected, the E–B leakage highly contaminates the estimation of the angular power spectrum over most of the multipole range of interest.

Fig. 4. B-mode CMB power spectrum reconstruction after the E–B leakage correction with the application of the recycling method and the extensions presented in Sect. 3.1.1. On the left: Mean angular power spectrum over Nsims = 200 CMB-only simulations that include lensing and a primordial tensor signal with r = 0.01. Solid lines represent the power spectrum of exact B maps (input), which were reconstructed with a full-sky B decomposition of Q and U and then masked for power spectrum estimation; power spectra of leakage-corrected maps (output LC) are shown with circles and diamonds for SWIPE and SO-SAT, respectively; the dashed lines are the leakage residuals after the correction; and the dashed-dotted lines (if present) are the CMB B modes without any leakage correction (output noLC). On the right: Effective tensor-to-scalar ratio associated with the absolute difference between leakage-corrected and exact angular power spectra. The top panel shows the results when the recycling method was applied to CMB B-mode maps for the SWIPE (red) and SO-SAT (blue) footprints; the middle panel reports those obtained when an iterative B decomposition with no (red), one (green), two (cyan), or three (magenta) iterations were performed. In the bottom panel, we compare the performance of iterative decomposition with three iterations (magenta) and diffusive inpainting (orange) when correcting the residual leakage. The middle and bottom panels present results only for SWIPE. The adopted binning scheme is Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 for the spectra computed in the LSPE-SWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The error bars highlight the uncertainty on the mean at 1σ, which was estimated from the dispersion of the angular power spectra of the simulations divided by $\sqrt{N_{\mathrm{sims}}}$. Readers can refer to the main text for details.

The leakage correction is performed on the full SWIPE and SO-SAT footprints; the angular power spectra are estimated, instead, masking the 4% of pixels closest to the borders, where the residual leakage is expected to still affect the reconstruction of the CMB B-mode power spectrum. It is possible to observe that the residuals after the application of the recycling method are negligible for the SO-SAT case in the entire multipole range of interest (30 ≲ ℓ ≲ 130), if compared to the upper bound on r targeted by the experiment (r ≲ 0.003 at 1σ, Ade et al. 2019). Indeed, the effective tensor-to-scalar ratio associated with the leakage residuals is at the level of r_leak ∼ 10⁻⁴ at all angular scales.

For SWIPE, the recombination bump can be exactly reconstructed given the sensitivity of the instrument (Δr ≲ 0.015 at 95% CL) with r_leak ∼ 10⁻⁴, while large angular scales still suffer relevant contamination due to residual leakage. We also observe an increasing trend of r_leak for ℓ ≥ 100, which is caused by two main effects:

• the angular resolution of the input CMB maps: on those multipoles, indeed, the instrumental beam begins to have an impact and, therefore, the correction coefficient of the recycling method is less sensitive to the E–B leakage contamination on those smaller angular scales;

• the CMB lensing B-mode signal becomes dominant while the expected primordial tensor spectrum decreases, and even a small relative error in the reconstruction results in a larger value of r_leak.

Given the capability of LSPE to measure primordial B modes on the largest angular scales, it would be desirable to lower the E–B leakage residuals in that multipole range with a further processing of the maps corrected with the recycling method. To this end, we implement two separate extensions:

• the iterative B decomposition,

• the diffusive inpainting.

Both proposed techniques are not applied to the SO-SAT dataset, since in this case optimal results in correcting the E–B leakage effect are already obtained with the standard recycling method.

Iterative B decomposition. is summarised in Fig. 5 and consists of iteratively reconstructing a B-family, $P_B^{(n)}={(Q_B,U_B)}^{(n)}$, by performing a B decomposition of the masked $P_B^{(n-1)}={(Q_B,U_B)}^{(n-1)}$. The starting set of Stokes parameters for the iterations, $P_B^{(1)}={(Q_B,U_B)}^{(1)}$, is the leakage-corrected B-family returned by the recycling method, $P_B^{″}={(Q_B,U_B)}^{″}$. With such iterations, the pure B modes in the map are preserved, while the residual ambiguous E modes are progressively de-projected, because interpreted as E in one of the iterations.

Fig. 5. Iterative B decomposition extension of the recycling E–B leakage correction method.

The results of the iterative B decomposition are shown in the second row of Fig. 4, where we plot the average angular power spectra of the 200 CMB-only simulations when either the recycling method alone (iter 0) or a post-processing with 1,  2 or 3 iterative B decompositions is applied. The plots highlight how, with this additional procedure, we are able to de-project the residual E–B leakage contamination in the B-mode CMB maps at most of the angular scales of interest, obtaining, when three iterations are performed, residuals at a much lower level (r_leak ∼ 10⁻⁴) than the LSPE targeted sensitivity already for ℓ > 4. On the other hand, iterative B decomposition does not lead to an improvement in the reconstruction for ℓ ≤ 4. This does not mean that we have not further removed residual E–B leakage contamination on those scales, but highlights another unavoidable phenomenon of the QU to EB decomposition of partial-sky observations. Both the first (in the standard recycling method) and the following B decompositions (in the iterations) of the input Q and U maps also suffer from B–E leakage. In each harmonic transformation, some ambiguous modes, which would be identified as B modes in a full-sky analysis, are erroneously estimated as E modes, causing a loss of power in the final B-mode map. This effect is particularly relevant on large scales and when E–B leakage residuals are very low. This explains our inability to correctly reconstruct the input power at ℓ ≤ 4 and the larger (but still irrelevant) errors in the map at ℓ ≥ 20 when several iterations are performed with respect to the case with just one. In these cases, the B–E leakage can become the limiting factor in our ability to reconstruct the input CMB B-mode signal, possibly leading to an under-estimate of the input tensor-to-scalar ratio. This is particularly meaningful for those experiments (such as LSPE) targeting the observation of the reionisation peak. It is important, then, to determine which multipoles are critically affected by this phenomenon and to exclude them from the cosmological analysis.

The effective tensor-to-scalar ratio fitted on the obtained residuals (see the right panels of Fig. 4) highlights that the first three multipoles (ℓ ≤ 4) suffer a significant B-E leakage in the SWIPE case (r_leak ∼ 7 × 10⁻³) and must be excluded. The hexadecupole has been estimated as the largest multipole to reject because including those modes in the second bin of the angular power spectra of the reconstructed CMB maps in Fig. 4 leads to a significant under-estimate of its value if compared to the input one.

The applicability of iterative B decomposition may present a drawback if the use of a large number of iterations leads to a considerable power loss due to B–E leakage along with the subtraction of residual ambiguous E-mode contamination. To assess whether this is the case, we perform the iterative B decomposition on two separate masked initial B-families: $P_B^{\mathrm{only}B}$ and $P_B^{\mathrm{leak}}$. The former is the set of Stokes parameters reconstructed with a full-sky B decomposition of input Q and U and represents the ideal goal of a leakage correction method. $P_B^{\mathrm{leak}}$, instead, is composed of Q and U maps sourced only by the residual E–B leakage after the application of the recycling method. In Fig. 6, we compare the angular power spectra of the output B-mode maps from the B decomposition performed on the two families after several iterations with those of the initial $P_B^{\mathrm{only}B}$ and $P_B^{\mathrm{leak}}$. It is possible to observe that for ℓ ≥ 5, iterations do not de-project modes of the B maps obtained from the initial $P_B^{\mathrm{only}B}$ and at the same time are able to subtract most of the leakage residuals (60% of the power already after the first iteration). At ℓ < 5, instead, the iterative B decomposition suffers loss of power in $P_B^{\mathrm{only}B}$ due to B–E leakage, but these scales would have been excluded anyway from the cosmological analysis of the NILC solution because an analogous significant loss is observed even without any iteration. Therefore, Fig. 6 highlights the benefit of employing the iterative B decomposition.

Fig. 6. LSPE-SWIPE only. Ratio between the mean angular power spectrum of either pure B modes only (solid lines) or E–B leakage residuals (dashed lines) after having applied the iterative B decomposition correction and that of corresponding input maps. The cases with 1 (red), 2 (blue), 4 (magenta), and 9 (green) iterations are shown. 200 different realisations of CMB have been considered, including lensing and a primordial tensor signal with r = 0.01. The trend of dashed lines highlights the reduction of E–B leakage residuals obtained by performing the iterations. The trend of solid lines, instead, the loss of power suffered due to the B–E leakage when applying the iterative B decomposition. A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed. The error-bars highlight the uncertainty on the mean at 1σ, estimated from the dispersion of the angular power spectra of the simulations divided by $\sqrt{N_{\mathrm{sims}}}$. Readers can refer to the main text for details.

Given the previous results, we apply the NILC pipeline on multi-frequency B-mode maps reconstructed in the SWIPE patch and processed with the recycling method and three iterations of B decomposition. Hereafter, we refer to it as rit-NILC.

Diffusive inpainting. The second proposed method to improve leakage removal on large scales for SWIPE is diffusive inpainting. This technique has already been introduced in Liu et al. (2019) as an alternative to the recycling method. It is based on the fact that ambiguous modes ψ satisfy the spherical bi-Laplacian equation:

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2({\mathrm{\nabla }}^2+2)\psi =0\end{array}$$(14)

subject to homogeneous Neumann and Dirichlet boundary conditions at the edge of the observed patch. An approximate solution of the equation above can be obtained by replacing the bi-Laplacian with the Laplacian and neglecting the Neumann boundary conditions. In this case, a template of residual ambiguous modes in the recycling CMB B-map is given by diffusive inpainting. The procedure consists in imposing as boundary conditions the values of the pixels of the recycling CMB solution at the edges of the patch and replacing iteratively inner sky pixels in the footprint with the average of their neighbours. The obtained template of the residuals is then subtracted from the input recycling CMB map.

The results of the application of the diffusive inpainting on the solutions of the recycling method for 200 CMB-only simulations with lensing+r = 0.01 on the LSPE-SWIPE patch are shown in the bottom panels of Fig. 4. They are compared with those obtained with iterative B decomposition and three iterations. Even in this case, we are able to significantly reduce the leakage residuals with respect to the standard recycling method for ℓ ≥ 5, while the angular scales at ℓ < 5 are too affected by the B–E leakage and have to be excluded from the subsequent cosmological analysis of the NILC solution. As for iterative B decomposition, just including ℓ = 4 in the second bin of the angular power spectra of the reconstructed CMB maps leads to a significant under-estimate of its value with respect to that of input B modes. The comparison of r_leak in Fig. 4 when diffusive inpainting or iterative B decomposition is applied highlights that the latter performs better on the largest scales (5 ≤ ℓ ≤ 20), while they are equivalent for intermediate and small ones (ℓ > 20).

Given these results, the NILC method is also applied to B-maps reconstructed in the SWIPE patch and corrected with the recycling method and the inpainting technique. Hereafter we refer to it as rin-NILC.

Finally, Fig. 7 shows that, independently of the amplitude of the input tensor perturbations, iterative B decomposition with three iterations and diffusive inpainting allow us to perfectly reconstruct the B-mode CMB anisotropies for ℓ ≥ 5. We generate CMB simulations with different input amplitudes of tensor perturbations and obtain leakage-corrected B-mode maps with the three techniques described above: the standard recycling method, iterative recycling, and inpainting recycling. Then, a tensor-to-scalar ratio is fitted on the average binned angular power spectrum ($C_{{\ell }_b}^{\mathit{BB}}$) of the output CMB B-mode maps for all different cases and excluding the first three multipoles (ℓ ≤ 4). The fit of the tensor-to-scalar ratio is performed with a Gaussian likelihood (Hamimeche & Lewis 2008; Gerbino et al. 2020):

$$\begin{array}{c}\hfill -2log\mathcal{L}(r)={\displaystyle \sum _{{\ell }_b,{\ell }_b^{\prime }}(C_{{\ell }_b}^{\mathit{BB}}-r{C}_{{\ell }_b}^{r=1}-C_{{\ell }_b}^{\mathrm{lens}})M_{{\ell }_b{\ell }_b^{\prime }}^{-1}(C_{{\ell }_b}^{\mathit{BB}}-r{C}_{{\ell }_b}^{r=1}-C_{{\ell }_b}^{\mathrm{lens}}),}\end{array}$$(15)

Fig. 7. LSPE-SWIPE only. Fitted tensor-to-scalar ratio on the mean angular power spectrum over 200 CMB-only simulations, including lensing and different values of the amplitude of the primordial tensor power spectrum. We consider the cases of B-mode maps reconstructed with: the recycling method (RM, red), the iterative recycling method with three iterations (magenta), the inpainting recycling method (orange), and without any E–B leakage contamination (ideal, cyan). The input values of r are represented with vertical dashed grey lines. The used likelihood is shown in Eq. (15) and the covariance is given only by the input CMB signal (ideal case). A binning scheme of Δℓ = 15 is employed and the first three multipoles (ℓ ≤ 4) are excluded. The error bars indicate the bounds at 1σ obtained from the posterior distribution.

where $C_{{\ell }_b}^{r=1}$ is the binned angular power spectrum of primordial tensor CMB B modes for a tensor-to-scalar ratio r = 1, $C_{{\ell }_b}^{\mathrm{lens}}$ the binned BB lensing angular power spectrum and $M_{{\mathcal{l}}_b{\mathcal{l}}_b^{\prime }}$ the covariance matrix associated with $C_{{\ell }_b}^{\mathit{BB}}$. A binning scheme of Δ_ℓ = 15 is adopted to make the angular power spectrum Gaussianly distributed.

We have considered an input r that varies between the sensitivity of LSPE and the current upper bounds reported in Tristram et al. (2022). In every case, with both iterative recycling and inpainting recycling methods, we are able to recover the same results obtained in an ideal case without any E–B leakage contamination in the maps. Instead, with the standard recycling method, a bias is observable, especially for lower values of the input tensor-to-scalar ratio, because of the residual E–B leakage contamination on large scales.

3.1.2. The ZB method

The second technique we consider for the correction of E–B leakage is the ZB method introduced in Zhao & Baskaran (2010). In this case, two different but related definitions for scalar and pseudo-scalar polarisation fields are employed:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{E}(\widehat{\gamma })=& -\frac{1}{2}[{\overline{ð}}_1{\overline{ð}}_2{P}_{+}(\widehat{\gamma })+{ð}_1{ð}_2{P}_{-}(\widehat{\gamma })]\hfill \\ \hfill \mathcal{B}(\widehat{\gamma })=& -\frac{1}{2i}[{\overline{ð}}_1{\overline{ð}}_2{P}_{+}(\widehat{\gamma })-{ð}_1{ð}_2{P}_{-}(\widehat{\gamma })],\hfill \end{array}\end{array}$$(16)

where ${\overline{ð}}_{\mathrm{s}}$ and ð_s are the spin-lowering and raising operators:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\overline{ð}}_{\mathrm{s}}=-{(sin\theta )}^{-s}(\frac{\partial }{\partial \theta }-\frac{i}{sin\theta }\frac{\partial }{\partial \varphi })[{sin}^s\theta ]\\ \hfill {ð}_{\mathrm{s}}=-{(sin\theta )}^{-s}(\frac{\partial }{\partial \theta }+\frac{i}{sin\theta }\frac{\partial }{\partial \varphi })[{sin}^s\theta ].\end{array}\end{array}$$(17)

The harmonic coefficients ℰ_ℓm and ℬ_ℓm of ℰ and ℬ are related to E_ℓm and B_ℓm of Eq. (13) by a ℓ-dependent numerical factor N_ℓ:

$$\begin{array}{c}\hfill {\mathcal{E}}_{\ell m}=N_{\ell }{E}_{\ell m},{\mathcal{B}}_{\ell m}=N_{\ell }{B}_{\ell m},N_{\ell }=\sqrt{(\ell +2)!/(\ell -2)!}.\end{array}$$(18)

In the case of partial-sky observations, the spin operators have to be applied on P₊ and P₋ fields multiplied by a window function $W(\widehat{\gamma })$, (non-zero only in the observed patch of the sky), and the recovered quantities $\stackrel{\sim }{\mathcal{E}}$ and $\stackrel{\sim }{\mathcal{B}}$ will be affected by the E–B leakage:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \stackrel{\sim }{\mathcal{E}}(\widehat{\gamma })=& -\frac{1}{2}[{\overline{ð}}_1{\overline{ð}}_2(P_{+}(\widehat{\gamma })W(\widehat{\gamma }))+{ð}_1{ð}_2(P_{-}(\widehat{\gamma })W(\widehat{\gamma }))]\hfill \\ \hfill \stackrel{\sim }{\mathcal{B}}(\widehat{\gamma })=& -\frac{1}{2i}[{\overline{ð}}_1{\overline{ð}}_2(P_{+}(\widehat{\gamma })W(\widehat{\gamma }))-{ð}_1{ð}_2(P_{-}(\widehat{\gamma })W(\widehat{\gamma }))].\hfill \end{array}\end{array}$$(19)

In Zhao & Baskaran (2010), the correction term to reconstruct the pure pseudo-scalar field ℬ in the observed region of the sky has been derived as:

$$\begin{array}{c}\hfill \mathcal{B}=\stackrel{\sim }{\mathcal{B}}{W}^{-1}+ct·W^{-2},\end{array}$$(20)

where:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill ct=& U[3cot\theta W{W}_x+W(W_{\mathit{xx}}-W_{\mathit{yy}})-2(W_x^2-W_y^2)]\hfill \\ \hfill & -Q[2cot\theta W{W}_y+2W{W}_{\mathit{xy}}-4W_x{W}_y]\hfill \\ \hfill & -2W_y[{(QW)}_x+{(UW)}_y]+2W_x[{(UW)}_x-{(QW)}_y],\hfill \end{array}\end{array}$$(21)

with $F_x=\frac{\partial F}{\partial \theta },F_{\mathit{y}}=\frac{\partial F}{sin\theta \partial \varphi },F_{\mathit{xx}}=\frac{{\partial }^{2}F}{\partial {\theta }^2},F_{\mathit{y}}=\frac{\partial F^2}{{sin}^2\theta \partial {\varphi }^2}$, and $F_{x\mathit{y}}=\frac{\partial F^2}{sin\theta \partial \varphi \partial \theta }$.

In this framework, as window functions, we adopt the footprints of SWIPE and SO (shown in Fig. 1) apodised with the ‘C¹’ scheme, where any observed pixel is multiplied by a weighting function f:

$$\begin{array}{c}\hfill f=\{\begin{array}{cc}x-sin(2\pi x)/(2\pi )\hfill & x<1\hfill \\ 1\hfill & \mathrm{otherwise}\hfill \end{array}\end{array}$$(22)

with $x=\sqrt{(1-cos\theta )/(1-cos({\theta }_{\ast }))}$, being θ_* the apodisation scale and θ the angular separation between the pixel itself and the closest unobserved one. Such an apodisation scheme ensures that the total uncertainty of the recovered power spectra is (nearly) minimal (see Grain et al. 2009). In our analysis, we adopt an apodisation length of 5° for both SWIPE and SO. In the ZB method, masks have to be apodised because spinorial derivatives of the polarisation field are badly estimated close to the borders of the patch if one adopts a binary mask.

The reconstruction of the pure ℬ field with this procedure is exact. However, the ill-behaved nature of W⁻¹ and W⁻² may complicate it, especially in the proximity of the unobserved region. Therefore, in this analysis, we have excluded all pixels where W ≤ 0.01 with a loss of sky fraction that corresponds to almost 1% for both footprints. In practise $\stackrel{\sim }{\mathcal{E}}$ and $\stackrel{\sim }{\mathcal{B}}$ are obtained by computing E_ℓm and B_ℓm of Eq. (13) with an harmonic decomposition of masked Q and U maps and then filtering them with N_ℓ.

Figure 8 compares the mean angular power spectrum of 200 leakage-corrected CMB-only simulations with lensing+r = 0.01 for the SWIPE and SO patches with that of the exact signal reconstructed with a full-sky E–B decomposition of the simulated Q and U maps. In these cases, all angular power spectra have been computed employing masks obtained from the exclusion of pixels with values lower than 0.01 in the original apodised patches and apodising them with a 7° apodisation scale, leading to final sky fractions of f_sky = 32% and 29%, respectively, for SWIPE and SO. In both footprints, the power of the residuals is lower across the entire multipole range of interest than the sensitivities targeted by the experiments, with associated values of the effective tensor-to-scalar ratio: 10⁻⁷ ≲ r_leak ≲ 10⁻³ (see the right panel in Fig. 8). With the ZB method we are even able to confidently reconstruct the primordial tensor B-mode power spectrum at ℓ < 5. This is motivated by the fact that, in this case, in contrast to the application of the recycling method, the correction term in Eq. (21) is exact and we can employ apodised masks which reduce both the E–B and the B–E leakage effects. Thus, these large angular scales are included in the cosmological analysis of the NILC CMB solutions when the ZB method is adopted.

Fig. 8. B-mode CMB power spectrum reconstruction after the E–B leakage correction with the application of the ZB method. On the left: Mean angular power spectrum over 200 CMB-only simulations including lensing and a primordial tensor signal with r = 0.01. Solid lines represent the power of exact B-maps reconstructed with full-sky EB-decomposition of Q and U; circles and diamonds that of leakage-corrected maps (output LC) with the ZB method; the dashed lines the leakage residuals after the correction; the dashed-dotted lines the CMB B modes without any leakage correction (output noLC). The LSPE-SWIPE (red) and SO-SAT (blue) footprints have been considered. On the right: Effective tensor-to-scalar ratio associated with the absolute difference between leakage-corrected and exact angular power spectra is shown for SWIPE (red) and SO-SAT (blue). A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPE-SWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The error-bars highlight the uncertainty on the mean at 1σ. Readers can refer to the main text for details.

3.2. Needlet filtering of partial-sky maps

The second step in the pre-processing of B-mode data which could introduce distortions in the CMB reconstruction before applying a NILC minimisation is the needlet filtering of partial-sky observations. Needlet coefficients, indeed, are computed by convolving a map with a filter in pixel-space. If observations only cover a partial fraction of the sky, the signal in the pixels close to the borders of the observed patch is mixed with the null values of the unobserved region by the convolution procedure. Thus, the estimation of needlet coefficients can be highly affected by this operation. The wider the needlet band is in pixel space, the larger this effect will be.

To assess the impact of needlet filtering, we generate 200 CMB B-mode maps from a full-sky E–B decomposition of Q and U, simulated with Planck 2018 best-fit parameters and lensing+r = 0.01. Then, we apply a forward and an inverse needlet transformation to these B-mode maps, properly masked to recover the SWIPE and SO footprints. We assess the goodness of the reconstruction by comparing the angular power spectra of maps obtained after a forward and inverse needlet transform, $D_{\ell ,\mathrm{nl}}^{\mathit{BB}}$, and those of the input ones, $D_{\ell ,\mathrm{in}}^{\mathit{BB}}$. In analogy to the analysis of E − B leakage correction, we consider an effective tensor-to-scalar ratio r_nl:

$$\begin{array}{c}\hfill r_{\mathrm{nl}}=|D_{\ell ,\mathrm{in}}^{\mathit{BB}}-D_{\ell ,\mathrm{nl}}^{\mathit{BB}}|/D_{\ell ,r=1}^{\mathit{BB}},\end{array}$$

which quantifies the absolute error we can make in the partial-sky reconstruction of both the tensor modes and the lensing signal due to the needlet filtering.

In this analysis and in the rest of the paper, we adopt Mexican needlets (Geller & Mayeli 2008), which are characterised by the following Gaussian-related weighting function:

$$\begin{array}{c}\hfill b_m(\frac{\ell }{B^j})=(\frac{\ell }{B^j}{)}^p·exp(-\frac{1}{2}·(\frac{\ell }{B^j}{)}^2).\end{array}$$(23)

We have set B = 1.5 and p = 1. The localisation properties of these wavelets in the real domain are better than those of standard and cosine needlets, because their tails are Gaussian, and therefore the convolution procedure is expected to be less affected by border effects.

The NILC CMB reconstruction can be distorted, especially on the largest scales, by two phenomena:

• if the first needlet harmonic band includes only few low multipoles, the localisation of the filter would be very poor and the needlet coefficients estimation would be largely impacted by the null values of the unobserved pixels;

• if only few modes are sampled by the first needlet band, the empirical covariance matrix of Eq. (8) would be highly uncertain and chance correlations between the CMB and other components could lead to the loss of some CMB power in the NILC reconstruction (Delabrouille et al. 2009).

To avoid these effects, we add together the first 11 bands of the employed Mexican needlets in the following way:

$$\begin{array}{c}\hfill b^{\mathit{new}}(\ell )=\sqrt{{\displaystyle \sum _{j=j_{\min }}^{j_{\max }}}{b}_j^2(\ell )},\end{array}$$(24)

leading to the configuration of needlet filters shown in Fig. 9.

Fig. 9. Needlet configuration employed in this work. Left panel: Mexican needlets bands with B = 1.5 in the harmonic domain. Right: Corresponding needlet filters in real space of the first five bands, plotted with the same colours of the associated harmonic bands.

Top panels of Fig. 10 show that filtering masked B-mode maps with such a needlet configuration does not introduce any significant error in the reconstruction of the CMB signal on a partial sky. Indeed, the residuals between the input maps and those obtained after a forward and an inverse needlet transform are at the level of r_nl ≲ 10⁻⁴ for both SWIPE and SO at all angular scales (see the upper right panel of Fig. 10).

Fig. 10. Goodness of the reconstruction of the CMB B-mode signal (lensing+r = 0.01) after needlet filtering (top) in the SWIPE (red) and SO (blue) patches or after smoothing CMB maps with input angular resolution of FWHM = 9′ to bring them to FWHM = 91′ in the SO footprint (bottom). On the left: Mean angular power spectra (over 200 CMB-only simulations) of the input exact B maps reconstructed with a full-sky EB-decomposition of Q and U (solid lines); those of the reconstructed maps after needlet filtering or beam convolution (dots and diamonds); and the corresponding residuals (dashed lines). The plot on the right represents the effective tensor-to-scalar ratios associated with the absolute differences between output power spectra after needlet filtering and beam convolution and exact input ones in the SWIPE (red) and SO (blue) patches. A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPE-SWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The error-bars highlight the uncertainty on the mean of the simulations at 1σ.

3.3. Smoothing effect

All ILC methods require input frequency maps with a common angular resolution, which is usually that of the channel with the largest beam. As anticipated in Sect. 2, we consider two datasets in this analysis: LSPE-SWIPE complemented with Planck frequencies (SWIPE+PLANCK) and SO-SAT.

For SWIPE+Planck it suffices to bring the full-sky Planck maps to the lower SWIPE angular resolution of FWHM = 85′.

For SO-SAT, instead, where frequency channels have different angular resolutions, such an operation is not trivial. The convolution for a beam of partial-sky maps can lead to an error in the CMB reconstruction, as in the case of needlet filtering, because the values of the pixels close to the border of the patch may be influenced by the null values outside the observed region.

To assess the relevance of this effect, we generate 200 CMB B-mode simulated maps with lensing+r = 0.01 at the highest angular resolution of SO (FWHM = 9′) and compare the average angular power spectrum, computed in the SO-SAT patch, of these maps when they are brought to the lowest resolution of SO (FWHM = 91′) full-sky ($D_{\ell ,\mathrm{in}}^{\mathit{BB}}$) or after having been masked ($D_{\ell ,\mathrm{smooth}}^{\mathit{BB}}$). The comparison between these spectra is shown in the bottom panels of Fig. 10 together with an assessment of the residuals through an effective tensor-to-scalar ratio:

$$\begin{array}{c}\hfill r_{\mathrm{smooth}}=|D_{\ell ,\mathrm{in}}^{\mathit{BB}}-D_{\ell ,\mathrm{smooth}}^{\mathit{BB}}|/D_{\ell ,r=1}^{\mathit{BB}}.\end{array}$$

It is possible to observe that the most aggressive convolution, which has to be applied to SO maps, does not considerably affect the B-mode CMB reconstruction with values of r_smooth of the order of or below 10⁻⁴ at all the angular scales of interest for B-mode science.

4. NILC performance: Analysis setup

After having applied a leakage correction method and needlet filtering to input partial-sky B-mode multi-frequency maps, the NILC application leads to the best blind estimate of the CMB B-mode signal in the considered patch. Such a solution, as shown in Eq. (5), is contaminated by residuals of the polarised foregrounds emission, especially close to the Galactic plane. For this reason, when computing the angular power spectrum for a cosmological analysis, the expected most contaminated regions are masked, lowering the effective sky fraction. The masking strategies adopted in this work are described in Sect. 4.1.

In this work, where we analyse simulated data, we can precisely determine the expected level of foregrounds and noise residual contamination in the NILC CMB solution by combining the input foregrounds and noise simulations at different frequencies with the NILC weights. Therefore, it is possible to easily assess the impact of these residuals on the estimate of the tensor-to-scalar ratio following the approach detailed in Sect. 4.2.

In addition to the residual foregrounds contamination, the other main issue of the NILC CMB reconstruction is the loss of modes induced by the empirical correlation of the reconstructed CMB with other components, especially on large scales, which might lead to a negative bias in the estimated power spectrum. However, we have verified that this effect is negligible in our analysis, given the adopted configuration of needlet bands (see Appendix A).

4.1. Foregrounds masks

In order to perform a cosmological analysis of a NILC CMB solution, masking the most foregrounds contaminated regions is usually needed before estimating the angular power spectrum. We have adopted two distinct and alternative methods to generate such masks:

1. computing the polarised intensity P² = Q² + U² of simulated foregrounds maps at a CMB frequency channel of the considered dataset, smoothing it with a Gaussian beam of FWHM = 12.5°, thresholding it to obtain masks with decreasing f_sky;

2. considering directly the average foreground residuals of the method under study, thresholding this map to obtain binary masks with decreasing f_sky, smoothing them with a Gaussian beam of FWHM = 8°, finally thresholding the smoothed masks to get binary masks with f_sky equal to the original ones.

In both cases, the considered SO-SAT patch is shown in Fig. 2, and corresponds, as already discussed in Sect. 2, to the region where the ILC methods are applied. The rationale of the first masking strategy is to have common masks for all the considered methods (NILC, ILC, and PILC). In such a way, we can assess and compare their foregrounds residuals in the same patch of the sky. These masks are obtained by assuming to have perfect knowledge of the polarised Galactic emission across the sky at a CMB channel. This is not possible when dealing with real data. However, similar masks can be obtained from foreground templates built by re-scaling with appropriate SEDs the Q and U maps of the considered dataset at the lowest (for synchrotron) and highest (for thermal dust) frequencies. The smoothing procedure is performed on the map to be thresholded in order to attenuate the small-scale fluctuations, which in real data are mainly induced by noise and CMB, and thus to avoid masks with a patchy structure.

Instead, the second masking approach is more method-specific. It is based on the assumption to be able in the future to better predict the distribution and intensity of foregrounds residuals in the sky even thanks to future surveys devoted to the studies of Galactic emission at microwave frequencies. This masking strategy is adopted to compute the angular power spectra of NILC CMB solutions and of its residuals when the E–B leakage is corrected with the recycling or ZB method.

When the ZB method is applied for the E–B leakage correction, the mask used for the power spectrum computation is apodised. Indeed, when dealing with second spinorial derivatives of a field on the sphere, the use of a binary mask would lead to an ill-behaved estimation of the power spectrum. We adopt the C¹ scheme (see Eq. (22)) with an apodisation length of 12° and 2° for SWIPE and SO, respectively. The choice of the apodisation lengths in the two cases is influenced by the minimum multipole (ℓ = 2 for SWIPE+Planck, ℓ = 30 for SO-SAT) that we consider for the computation of the angular power spectra.

In the case of the NILC application on SWIPE+PlanckB-mode maps corrected with the recycling method, we exclude the 4% of the pixels closest to the patch border, where some residual leakage could affect the estimate of the power spectra. This step, instead, is not needed for SO-SAT, because the reduced patch of Fig. 2 already excludes the borders of the original footprint. In this second masking strategy, the sky fractions for each case have been determined by comparing the level of the foregrounds residuals with the primordial tensor signal targeted by the experiment.

4.2. Estimates of the residuals contamination

When NILC is applied on simulated datasets, we have access to its foregrounds and noise residuals in the B-mode solution. Therefore, it is possible to compare their angular power spectra with a primordial BB power spectrum of tensor modes with r equal to the target of the experiment under consideration.

If cross correlations among components are negligible (as proved in the Appendix A), the power spectrum of the cleaned map ($C_{\ell }^{\mathrm{out}}$) can be written as follows:

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{out}}=C_{\ell }^{\mathit{cmb}}+C_{\ell }^{\mathit{fgds}}+C_{\ell }^{\mathrm{noi}},\end{array}$$(25)

where $C_{\ell }^{\mathit{fgds}}$ and $C_{\ell }^{\mathrm{noi}}$ are the angular power spectra of residual foregrounds and noise. The noise bias, $C_{\ell }^{\mathrm{noi}}$, can be subtracted at the power spectrum level by estimating its contribution through Monte Carlo simulations or by computing cross-spectra between uncorrelated splits of the data (e.g. maps from subsets of detectors or from partial-mission observations). Hence, the main systematic bias to the final angular power spectrum will be given by the residual Galactic contamination.

To assess the impact of this residual term on the estimation of r from $C_{\ell }^{\mathrm{out}}$, we fitted r on binned $C_{{\ell }_b}^{\mathit{fgds}}$ assuming a Gaussian likelihood (Hamimeche & Lewis 2008; Gerbino et al. 2020):

$$\begin{array}{c}\hfill -2log\mathcal{L}(r)={\displaystyle \sum _{{\ell }_b,{\ell }_b^{\prime }}(C_{{\ell }_b}^{\mathit{fgds}}-r{C}_{{\ell }_b}^{r=1})M_{{\ell }_b{\ell }_b^{\prime }}^{-1}(C_{{\ell }_b^{\prime }}^{\mathit{fgds}}-r{C}_{{\ell }_b^{\prime }}^{r=1}),}\end{array}$$(26)

where $C_{{\ell }_b}^{r=1}$ is the primordial tensor CMB B-mode angular power spectrum for a tensor-to-scalar ratio r = 1 and $M_{{\mathcal{l}}_b{\mathcal{l}}_b^{\prime }}$ is the covariance matrix for a fiducial cosmological model with r = 0. An aggressive binning makes the angular power spectrum Gaussianly distributed and mitigates correlations among different modes (induced by the use of a mask). In this work, we adopt a constant binning of Δ_ℓ = 15. We have tested the robustness of the obtained constraints by both varying the binning scheme and even considering an inverse-Wishart likelihood function.

The covariance matrix is estimated from the angular power spectra of 200 NILC simulated solutions as:

$$\begin{array}{c}\hfill M_{{\ell }_b{\ell }_b^{\prime }}=\mathrm{Cov}(C_{{\ell }_b}^{\mathrm{out}}-(1-A_L)·C_{{\ell }_b}^{\mathrm{lens}},C_{{\ell }_b^{\prime }}^{\mathrm{out}}-(1-A_L)·C_{{\ell }_b}^{\mathrm{lens}}),\end{array}$$(27)

where A_L quantifies our ability to de-lens the output map with A_L = 1 for no de-lensing and A_L = 0 for full de-lensing. Equation (27) accounts for the cosmic variance associated with the residual lensing signal, the sample variance of the residual foreground and noise angular power spectra and all the cross-terms. In this work, we assume A_L = 0.5 (which is the level that can be achieved with multi-tracers techniques; see e.g. Yu et al. 2017).

5. Results

In the processing of future ground-based and balloon-borne CMB polarisation experiments, NILC represents an effective alternative to the commonly used parametric methods, given the very low number of assumptions on which the algorithm is based. In Sect. 5.1, we assess the ability of NILC to recover the input CMB B-mode signal in an ideal case where we assume to be able to perfectly reconstruct the needlet and pixel coefficients of B-mode data even for partial-sky observations. We compare these results with those obtained with the application of ILC and PILC.

Once we have verified that it is worth applying NILC on B modes, we proceed to report in Sect. 5.2 the results of the application of NILC to a more realistic case with input maps leakage-corrected with the recycling or ZB method. In this case, also the needlet filtering and the beam convolution are performed on partial-sky maps.

In all the following analyses, the input CMB B-mode signal has no tensor modes and only lensing. We have verified that the same conclusions hold when injecting a primordial inflationary tensor signal.

We consider two different datasets:

• the three channels of SWIPE, complemented with the seven polarised Planck frequency maps (SWIPE+Planck);

• the six simulated frequency maps of the SO-SAT.

The elements of the covariance matrix of Eq. (8) associated with the Planck channels are computed considering only the observed SWIPE region (shown in Fig. 1). For SO, the considered patch is shown in Fig. 2 and represents a faithful approximation of the SO-SAT sky coverage under the assumption of homogeneous noise. For the application of NILC, the used configuration of needlet bands in all cases is the same of the test performed in Sect. 3.2: Mexican needlets with B = 1.5 and the first eleven bands added together.

5.1. NILC for partial-sky B modes data: Ideal case

At the moment, PILC (see Sect. 3) is the only blind method available to easily subtract foregrounds emission from ground-based and balloon-borne CMB polarisation data. However, it works with polarisation intensity data, which account for both E and B modes; therefore, the application of ILC and NILC algorithms could represent a valid alternative, since they are sensitive only to Galactic contamination in B modes.

We compare the performance of PILC, ILC, and NILC by looking at the angular power spectra of their foregrounds and noise residuals in an ideal case, where we assume to be able to perfectly recover the B-mode signal in needlet and pixel space for partial-sky observations. This is explicitly obtained by reconstructing from simulated full-sky Q and U maps full-sky B-mode maps, which are then smoothed to a common angular resolution and needlet filtered, and only at the end they are masked according to the patch of the considered experiments. In all cases, a model of foregrounds with spatially varying spectral indices is assumed for polarised dust and synchrotron emissions (d1s1).

This ideal case with a full-sky reconstruction of B-mode maps is considered just to compare the performance of the different blind methods and to assess the masking strategy needed to reach the scientific targets of the experiments under study. Then in Sects. 5.2.1 and 5.2.2 we report the results for a realistic approach with a reconstruction of B-mode maps from partial-sky Q and U observations.

In Fig. 11, the angular power spectra of foregrounds and noise residuals from the application of PILC, ILC, and NILC are shown. They have been computed using Galactic masks retaining the 15% (for SWIPE) and 7.5% (for SO) of the sky and generated with the first procedure described in Sect. 4.1.

Fig. 11. Mean angular power spectra (over 200 simulations) of foregrounds (solid lines) and noise (dashed lines) residuals when PILC (green), ILC (blue) or NILC (red) are applied on B-mode data (reconstructed on the full-sky and then masked) of SWIPE+Planck (left) or SO (right) frequency maps. The spectra have been computed considering the 15% (for SWIPE) and 7.5% (for SO) of the sky which is expected to be less contaminated by Galactic emission in the CMB frequency channels of SWIPE and SO. The input foregrounds emission is simulated with the d1s1 sky model. A binning scheme of Δℓ = 10 is employed.

We start by considering the SWIPE dataset complemented with Planck frequencies (SWIPE+Planck) as introduced in Sect. 2. A method which minimises the variance directly of B modes (ILC) instead of polarised intensity (PILC) leads to lower residuals of Galactic emission and instrumental noise at all angular scales. The CMB, indeed, is much brighter in $P=\sqrt{Q^2+U^2}$ (due to the presence of the E-mode signal) than in B and therefore the PILC weights are less capable of tracing and subtracting the frequency-dependent contaminants in the B-mode sky. Furthermore, extending the procedure of variance minimisation to the needlet domain (NILC) yields significantly lower foregrounds residuals, especially on large scales, with respect to ILC and PILC at the price of a mild increase of the noise residuals. This justifies the implementation of the NILC algorithm for cut-sky observations given the sensitivity to the tensor-to-scalar ratio in that multipole range.

The NILC effectiveness is particularly relevant for SWIPE because we have access to information at very low multipoles. For SO, which will observe only modes at ℓ ≥ 30, indeed, NILC foregrounds residuals are comparable with those of ILC and PILC. Indeed, the range of multipoles where NILC performs better in subtracting Galactic contamination is highly dependent on the sensitivity and angular resolution of the considered dataset. However, as expected, the power of total residuals is lower for NILC than for ILC and PILC across all multipoles, because the output variance is separately minimised at different needlet scales.

The above results support the effort to extend the application of NILC on B-mode data for ground-based and balloon-borne experiments.

5.2. NILC results for realistic SWIPE and SO case studies

A realistic application of NILC on partial-sky B-mode maps can be summarised with the following steps:

• Correct the E–B leakage effect in the multi-frequency dataset, which include CMB, noise, and Galactic foregrounds

• where needed, bring all the frequency maps to the same resolution

• filter the maps with needlet bands

• perform the variance minimisation independently on the different needlet scales

• reconstruct a cleaned CMB map in real space.

Depending on the technique employed to correct the E–B leakage, we can identify several distinct pipelines:

• r-NILC, for maps corrected with the recycling method,

• rit-NILC, when the recycling method and iterative B decomposition are applied,

• rin-NILC, for maps corrected with the recycling method and diffusive inpainting,

• ZB-NILC, when the ZB technique is performed.

We have verified that all the leakage-correction methods presented in Sect. 3.1 perform analogously in the correction of the E–B leakage even in the presence of instrumental noise and Galactic foregrounds.

We subsequently report the NILC results obtained when recycling (and its extensions) and when the ZB methods are applied, for both SWIPE+Planck and SO-SAT, in Sects. 5.2.1 and 5.2.2, respectively. Considering the results obtained in Sect. 3.1 regarding the effectiveness of the different methods in correcting E–B leakage in the two footprints, we apply rit-NILC, rin-NILC, and ZB-NILC to the SWIPE+Planck dataset, while r-NILC and ZB-NILC to SO-SAT simulated maps. In the case of the application of rit-NILC and rin-NILC to SWIPE+Planck, the minimum considered multipole is ℓ = 5 both in the plots and for the cosmological analysis, given the significant loss of power on larger angular scales due to B–E leakage, as found in Sect. 3.1.1. In the SO-SAT analysis, the leakage correction methods are performed within the full SO footprint (shown in the right panel of Fig. 1), since CMB data will be available in that region. The NILC pipelines, instead, are applied considering the patch of Fig. 2, which retains 10% of the sky and takes into account the highly inhomogeneous scanning strategy of the telescope.

The reference sky fraction for the cosmological analysis in each case is established by looking at the results in Fig. 12. In this analysis, we perform a fit of r on the power spectra of NILC foregrounds residuals (see Sect. 4.2) in the ideal case considered in Sect. 5.1 adopting the second masking approach described in Sect. 4.1 and varying the sky fraction. Figures 13–15 provide a detailed comparison of the two approaches. We analyse simulated data of:

Fig. 12. Trend of the tensor-to-scalar ratio rfgds fitted on the NILC foregrounds residuals angular power spectra with the likelihood of Eq. (26) when the sky fraction of the employed mask varies. The adopted masking strategy is the second one described in Sect. 4.1: we mask the most contaminated regions by foregrounds residuals. In this case, the input multi-frequency B-mode maps are constructed full-sky and then masked. The considered cases are SWIPE+Planck (left) with both d0s0 (red) and d1s1 (green) Galactic models and SO-SAT dataset (right) assuming only the d1s1 sky model. The grey shaded areas show the targeted sensitivity of the LSPE and SO experiments. The angular power spectra in the likelihood are binned with Δℓ = 15. The error bars indicate the bounds at 2σ (left) and 1σ (right) obtained from the posterior distributions.

• SWIPE+Planck with d0s0 Galactic model,

• SWIPE+Planck with d1s1 Galactic model,

• SO-SAT with d1s1 Galactic model.

The case with constant spectral indices for polarised thermal dust and synchrotron emission is considered to have a direct comparison with the recent forecast of the sensitivity on r for the LSPE experiment published in Addamo et al. (2021).

It is possible to observe that when NILC is applied on d0s0 simulations of the SWIPE+Planck dataset, no masking of the foregrounds residuals is needed to achieve the scientific goal of the mission in terms of sensitivity on the tensor-to-scalar ratio r. When a more realistic Galactic model (d1s1) is instead considered for the same dataset, an aggressive masking strategy is needed. The computation of the power spectrum in approximately the 12% less contaminated fraction of the sky allows us to obtain an upper bound at 95% CL on r lower than the sensitivity targeted by the experiment. Considering larger portions of the sky with such a blind method would lead to a bias on r at the level of 1.5 × 10⁻² or greater. When, instead, NILC is applied on ideal simulated B-mode data of SO-SAT, with a sky fraction of 8%, the foregrounds contamination is sensibly lower than the primordial tensor signal targeted by the experiment, even assuming a Galactic foreground model with varying spectral indices across the sky.

In the following, for all cases, we show the posterior distributions of an effective tensor-to-scalar ratio fitted on the foregrounds residuals average power spectrum, as described in Sect. 4.2. When r_fgds = 0 is within 2σ, we report upper bounds on r_fgds at 95% CL (for SWIPE) or 68% CL (for SO) to be compared with those in the official forecasts of the two experiments (see Addamo et al. 2021 and Ade et al. 2019).

5.2.1. Recycling NILC (r-NILC)

In this section we report the results obtained when the NILC algorithm is applied to multi-frequency B-mode maps that include CMB, noise, and foregrounds, and have been corrected with the recycling method (r-NILC). We note that, in the case of the LSPE+Planck dataset, the maps have also been post-processed either with the iterative B decomposition with three iterations (rit-NILC) or with the diffusive inpainting (rin-NILC).

Recently, updated forecasts on r for the LSPE experiment have been published in Addamo et al. (2021) considering a Galactic model with constant spectral indices for polarised thermal dust and synchrotron. Following the same approach, we first test the recycling NILC on the SWIPE+Planck dataset adopting the d0s0 model.

The mean angular power spectrum of foregrounds and noise residuals in this case are shown in the upper left panel of Fig. 13. The adopted binning scheme is Δℓ = 6 for 5 ≤ ℓ ≤ 28 and Δℓ = 15 for ℓ ≥ 29. The foregrounds residuals are confidently below both the reionisation and recombination bumps for r = 0.015 (the LSPE targeted upper bound at 2σ in case of no detection). The angular power spectra have been estimated without any masking of the most contaminated regions by residuals of the Galactic emission.

Fig. 13. NILC residuals for SWIPE+Planck dataset with the d0s0 Galactic model. On the left: Mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations. On the right: Posterior distributions of an effective tensor-to-scalar ratio fitted on the foregrounds residuals in the case of half de-lensing (AL = 0.5). For the estimation of the posterior a binning scheme of Δℓ = 15 has been used to make the angular power spectrum Gaussianly distributed (see Sect. 4.2). On the top: Blue and red lines represent the results when E − B leakage is corrected with, respectively, the iterative and inpainting recycling methods (rit-NILC and rin-NILC). On the bottom: Results when the ZB- (blue) and hybrid ZB-NILC (red) methods are applied. The adopted binning scheme is Δℓ = 6 for 5 ≤ ℓ ≤ 28 and Δℓ = 15 for ℓ ≥ 29 on the top left, while Δℓ = 6 for 2 ≤ ℓ ≤ 19 and Δℓ = 15 for ℓ ≥ 20 on the bottom left. In both cases the chosen binning strategy is just for visualisation purposes. Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.015, 0.03]. The lower bound represents the targeted upper limit at 95% CL by LSPE in case of no detection, the upper one the targeted detection at 99.7% CL. The spectra have been estimated without any masking of the regions most contaminated by foregrounds. The final fsky is 32% in the case on top, because 4% of the pixels closest to the borders are masked to avoid residual E–B leakage effects, and 28% in the case on the bottom due to the needed apodisation of the mask.

The plot on the upper right-hand side of Fig. 13 shows that the effective tensor-to-scalar ratio associated with foregrounds residuals has a posterior distribution with an upper limit at 95% CL of ∼1.5 × 10⁻², which is fully in accordance with the LSPE target. This constraint has been obtained by excluding the multipoles ℓ ≤ 4 where a significant negative bias is observed in the reconstructed CMB B-mode power spectrum due to the B–E leakage effect (see Sect. 3.1.1).

When a more complicated sky model is considered (d1s1 with anisotropic spectral indices for dust and synchrotron), the NILC weights try to trace locally several different SEDs at the same time, leading to a more noisy and less effective removal of the Galactic emission in each pixel. Therefore, in this case, a more aggressive masking strategy (based on the second approach described in Sect. 4.1) is needed to exclude the most contaminated regions. We have found in the analysis of Fig. 12 that a sky fraction of f_sky = 12% is enough to avoid biases in the estimate of the tensor-to-scalar ratio.

The angular power spectra of the residuals are shown in the upper left panel of Fig. 14. Both the rit-NILC and rin-NILC cleaning methods still allow one to possibly detect both peaks in the B-mode primordial spectrum. The posterior distribution of an effective tensor-to-scalar ratio fitted on the foregrounds residuals has an upper limit of around r ≲ 2.5 × 10⁻² at 2σ, which is mainly sourced by the contribution of the noise residuals to the covariance matrix of Eq. (27) and the exclusion of the first three multipoles from the cosmological analysis.

Fig. 14. NILC residuals for SWIPE+Planck dataset with the d1s1 Galactic model. On the left: Mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations. On the right: Posterior distributions of an effective tensor-to-scalar ratio fitted to the foregrounds residuals in the case of half de-lensing (AL = 0.5). For the estimation of the posteriors, a binning scheme of Δℓ = 15 has been used to make the angular power spectrum Gaussianly distributed (see Sect. 4.2). On the top: Blue and red lines represent the results when E–B leakage is corrected with, respectively, the iterative and inpainting recycling methods (rit-NILC and rin-NILC). On the bottom: Results when the ZB- (blue) and hybrid ZB-NILC (red) methods are applied. The adopted binning scheme is Δℓ = 6 for 5 ≤ ℓ ≤ 28 and Δℓ = 15 for ℓ ≥ 29 on the top left, while Δℓ = 6 for 2 ≤ ℓ ≤ 19 and Δℓ = 15 for ℓ ≥ 20 on the bottom left. In both cases the chosen binning strategy is just for visualisation purposes. Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.015, 0.03]. The lower bound represents the targeted upper limit at 95% CL by LSPE in case of no detection, the upper one the targeted detection at 99.7% CL. The spectra have been estimated masking the regions most contaminated by foregrounds for a final fsky of 12%.

The very same pipeline can be applied to SO-SAT simulations. In this case, the standard recycling method is adopted, given its effectiveness in removing E–B leakage contamination at the angular scales of interest for SO (see Sect. 3.1.1). The target of SO is the observation of only the recombination bump, as the experiment is not sensitive to multipoles ℓ < 30 due to atmospheric contamination. Therefore, the needlet bands and the cleaning algorithm do not take into account modes on angular scales larger than almost 6°.

For SO-SAT we have considered only a Galactic emission simulated with the d1s1 model. The residuals and the posterior distribution of the tensor-to-scalar ratio are shown in the top panels of Fig. 15. The goal of SO is a detection at high significance of r = 0.01. The application of r-NILC leads to foregrounds residuals that are below such a level in all the multipole range of interest, yielding an upper bound in the posterior distribution of r < 0.0027 (68% CL). This result is fully in agreement with the constraint targeted by the SO Collaboration (δr ∼ 0.003, Ade et al. 2019). However, in this analysis, we have considered input simulations with simple white and isotropic noise, since we simply aim to validate a partial-sky extension of the NILC pipeline to B-mode data from future experiments. A proper SO forecast should also consider a 1/f noise component, which is, however, expected to have only a mild impact on the effectiveness of the cleaning method and on the width of the posterior distribution. The above results for the SO case have been obtained considering the cleanest 8% of the sky, as suggested by the analysis reported in Fig. 12.

Fig. 15. NILC residuals for SO-SAT dataset with the d1s1 Galactic model. On the left: Mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations. On the right: Posterior distributions of an effective tensor-to-scalar ratio fitted to the foregrounds residuals in the case of half de-lensing (AL = 0.5). Both for the plots of the angular power spectra and the estimation of the posterior a binning scheme of Δℓ = 15 has been used. On the top: Results when E–B leakage is corrected with the recycling method (r-NILC). On the bottom: Results when the ZB- (blue) and hybrid ZB-NILC (red) methods are applied. Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.003, 0.01]. The lower bound represents the targeted upper limit at 68% CL by SO in case of no detection, the upper one the targeted detection at few σ significance. The spectra have been estimated masking the regions most contaminated by foregrounds for a final fsky of 8%.

5.2.2. ZB-NILC

The ZB-NILC pipeline consists of implementing the variance minimisation in needlet domain after having applied the ZB leakage correction method on multi-frequency B-mode maps that include CMB, noise, and foregrounds. Therefore, the cleaning algorithm is applied to ℬ maps of Eq. (16) instead of the usual B field. The standard B-mode signal is recovered only at the end when the angular power spectra are computed.

We consider the same datasets and cases of the r-NILC analysis to closely compare the results of the two methodologies by looking at Figs. 13–15. The mean angular power spectra (over 200 simulations) of the foregrounds and noise residuals are shown with blue solid and dashed lines, respectively.

In the bottom left panel of Fig. 13, the results for the Planck+SWIPE dataset with the d0s0 Galactic model are reported. It is possible to observe that the ZB-NILC foregrounds residuals are much higher than those of the recycling NILC. This is caused by the fact that, in ZB-NILC, the input maps are convolved for the harmonic filter $N_{\ell }=\sqrt{(\ell +2)!/(\ell -2)!}$ (see Sect. 3.1 for details), which boosts the power on smaller scales. Therefore, even the weights of the first needlet band are sensitive much more to contamination on intermediate scales (where noise has a not negligible contribution) than to that on the largest ones (where diffused foregrounds dominate). As a consequence, ZB-NILC weights are much less capable of removing Galactic polarised emission in ℬ maps.

To obtain an acceptable level of foregrounds residuals even in the case of the application of the NILC pipeline to ZB-corrected maps, we have implemented the Hybrid ZB-NILC (HZB-NILC). HZB-NILC consists of correcting the E–B leakage effect in frequency maps with the ZB method, but then combining them using the NILC weights estimated from the B-mode maps processed with the standard recycling method.

The foregrounds and noise residuals of the HZB-NILC are shown with solid and dashed red lines in the bottom panels of Figs. 13–15. With such an implementation, we can recover foregrounds and noise residuals at a level comparable to recycling NILC, but at the same time we can exploit the greater leakage correction capabilities of the ZB technique on the largest angular scales (ℓ ≤ 4).

In this case, for SWIPE+Planck with the d0s0 model, we obtain an upper bound on an effective tensor-to-scalar ratio fitted on the foregrounds residuals of 5 × 10⁻³ at 2σ that is fully compatible with the LSPE target (see bottom right panel of Fig. 13). The improved constraint on r with respect to r-NILC is motivated by the inclusion in this analysis of the lowest multipoles ℓ ≤ 4.

The upper limit on the tensor-to-scalar ratio obtained with HZB-NILC is sensibly lower than that obtained in Addamo et al. (2021) even though the same Galactic model is assumed. Such a difference can be explained taking into account three key factors: (i) in this work we do not take into account the scanning strategy and, thus, the anisotropic distribution of the instrumental noise, which may play a role in the performance of the component separation; (ii) we assume to be able to de-lens our B-mode power spectrum at a 50% level (Yu et al. 2017), thus reducing the impact of lensing B modes to the total uncertainty in the r reconstruction; (iii) for this simplistic Galactic model the minimisation of the variance across the different needlet scales can in principle outperform the parametric methods, which fit spectral parameters in different regions of the sky.

The ZB and hybrid ZB-NILC methods are then applied to the LSPE+Planck dataset assuming the more realistic d1s1 model of the Galaxy. The corresponding results are shown in the lower panels of Fig. 14. We again observe different capabilities of the two methods in subtracting Galactic contamination (in favour of the HZB-NILC). However, in this case, the differences between the power spectra are less evident because the considered sky components are intrinsically more difficult to subtract. With a more aggressive Galactic mask with respect to the d0s0 analysis, which retains the 12% of the sky, both the reionisation and recombination peaks with r = 0.015 would be observed. The upper bound on an effective tensor-to-scalar ratio fitted on the foregrounds residuals is r < 0.016 (95% CL) for the Hybrid ZB technique, which is in accordance with the LSPE target in case of no detection. As above, the lower constraint with respect to r-NILC pipelines is given by the possibility of including the first three multipoles in the likelihood.

The same trends are obtained when ZB-NILC and Hybrid ZB-NILC have been applied to the simulated SO-SAT dataset with Galactic emission simulated with the d1s1 model (see the bottom panel of Fig. 15). The angular power spectrum of the HZB-NILC Galactic residuals is below the curve of the primordial tensor signal targeted by SO (r = 0.01) for ℓ > 30, if we observe the cleanest 8% of the sky. The HZB-NILC foregrounds residuals lead to an upper bound in the posterior distribution of the effective tensor-to-scalar ratio of r < 0.0028 (68% CL) which is compatible with the SO target (r ≤ 3 × 10⁻³ at 1σ).

We observe that for SO-SAT the difference in the amplitudes of foregrounds residuals on large scales between ZB- and HZB-NILC is much less evident than the one for SWIPE+Planck. This is due to the fact that, without modes with ℓ ≤ 30, the calculation of HZB-NILC weights in the first needlet band is more affected by noise contamination on intermediate scales. This condition is close to that experienced by the ZB-NILC pipeline.

6. Conclusions

The main goal of future CMB experiments will be the detection of polarisation B modes generated by primordial tensor perturbations to definitely confirm the inflationary scenario. Given its low number of assumptions, NILC represents an effective alternative for the analysis of future B-mode data to the commonly used parametric component separation methods, which are more subject to systematic biases in the case of mis-modelling of the foreground properties. However, most future surveys will be balloon-borne or ground-based and will observe only a portion of the sky, whereas, at present, NILC has been successfully applied to full-sky data from either Planck or WMAP.

In this work, we explore the possibility to extend the NILC formalism to future B-mode partial-sky observations, specifically addressing the complications that such an application yields such as the E–B leakage, needlet filtering, and beam convolution (see Sect. 3). Their impact on the reconstruction of CMB B modes has been assessed in Monte Carlo simulations for two complementary experiments: LSPE-SWIPE and SO-SAT. The former aims at a detection with high significance of both reionisation and recombination peaks with r = 0.03 or to set an upper bound of r = 0.015 at a 95% confidence level; the latter, instead, is designed to be able to measure a signal at the r = 0.01 level at a few σ significances, or exclude it at a similar significance, using the amplitude of the B modes around the recombination bump. In the case of LSPE-SWIPE, realistic simulated Planck maps have been included in the component separation pipeline to have a broader frequency coverage and better trace the spectral properties of B-mode foregrounds.

When dealing with partial-sky observations, the estimation of B modes from CMB polarisation measurements is challenging due to the mixing of E and B modes. In order to correct for this effect, we implemented two different techniques which are to be applied in pixel space: the recycling (Liu et al. 2019) and the ZB (Zhao & Baskaran 2010) method. We have tested their ability to recover the CMB B-mode angular power spectrum on Monte Carlo simulations, also quantifying a possible bias on the estimate of the tensor-to-scalar ratio. We find that:

• The recycling method is able to reduce the E–B leakage contamination in the B-mode maps at a negligible level for SO at all the angular scales of interest (ℓ≥30), while for LSPE-SWIPE it is only able to do so at ℓ ≳ 15 (see Fig. 4).

• The ZB method, instead, allowed us to recover the input CMB signal at all angular scales for both the considered experiments (see Fig. 8). However, this technique involves some complications: the need to work with the ℬ field of Eq. (16) and to apodise the mask.

• The E–B leakage residuals of the recycling method on large angular scales (ℓ ≲ 15) in the SWIPE patch can be further reduced for ℓ ≥ 5 by applying either of the following two post-processing prescriptions:

  • diffusive inpainting (Liu et al. 2019)

  • iterative B decomposition, which we introduced for the first time in this paper. Ambiguous modes, which are still present in the map, were filtered out by iteratively applying the B decomposition of the Q and U CMB maps reconstructed with the recycling method.

• The reconstruction of the CMB B modes obtained with the recycling method and its extensions fails at ℓ < 5 due to the unavoidable loss of modes caused by the leakage of B modes into E modes. Therefore, when this method is applied, the first three multipoles have to be excluded from the final cosmological analysis.

Needlet filtering and beam convolution performed on incomplete sky observations can also potentially affect the reconstruction of CMB B modes, as values of the pixels close to the border of the observed patch are influenced by the null ones of the unobserved regions. However, for the experiments under consideration, we have verified that these operations have a negligible impact. We summarise the ability to recover the angular power spectrum of CMB B modes at different angular scales after applying the E–B leakage correction methods, needlet filtering, and beam convolution in Table 2.

Once the above pre-processing steps were tested, the performance of the NILC cleaning method was assessed on simulated multi-frequency B-mode datasets which include CMB with only lensing, instrumental noise, and Galactic foreground emission. For LSPE-SWIPE we have considered two different foreground models: one in which the spectral indices of the synchrotron and dust emissions are assumed to be position independent (d0s0), and another in which they vary across the sky (d1s1). For SO-SAT, instead, we generated simulations only with the d1s1 model, as in Ade et al. (2019).

The purpose of this paper is to introduce and validate a complete extension of the NILC pipeline, which is to be applied on partial-sky B-mode data for future experiments. A detailed forecast of the performance of the considered CMB experiments should take some additional realistic effects into account, such as the scanning strategy of the instrument or the 1/f noise contamination.

Simulated maps have been corrected for the E–B leakage with the following methods:

• recycling + iterative decomposition (rit-NILC) for SWIPE+Planck,

• recycling + diffusive inpainting (rin-NILC) for SWIPE+Planck,

• recycling (r-NILC) for SO-SAT, and

• ZB (ZB-NILC) for both SO-SAT and SWIPE+Planck.

For all cases, we quantified the ability of the pipeline to mitigate foreground contamination in terms of an effective tensor-to-scalar ratio, which was fitted on the angular power spectrum of foreground residuals through a Gaussian likelihood (see Eq. (26)) after the application of an appropriate masking strategy. The derived upper bounds are summarised in Table 3. They are reported at 95% and 68% CL for LSPE-SWIPE and SO, respectively, as in the official forecasts of the two experiments (see Addamo et al. 2021 and Ade et al. 2019).

We find that SO-SAT bounds are within the goal of the experiment regardless of the adopted E–B leakage correction approach. In this work, as mentioned above, we did not incorporate a 1/f noise component in the SO-SAT simulations, which is necessary for an accurate description of the experiment. However, this component is expected to only mildly affect the performance of the method and the obtained constraints, as also proved in Ade et al. (2019). For SWIPE+Planck with the d1s1 foregrounds model, when the recycling method is applied to correct the E–B leakage, the obtained upper bound on r is greater than the target of the experiment. This result is caused by the need to exclude the first three multipoles where the B–E leakage has a significant impact on the CMB power spectrum reconstruction. For all other cases, instead, the constraints are also in agreement with the expected sensitivity for LSPE-SWIPE.

When the ZB-NILC pipeline is applied (see Sect. 5.2.2), the amplitude of the foreground residuals is much greater than that obtained with r-NILC. This can be explained by the fact that in this case the input maps are convolved with the harmonic filter $N_{\ell }=\sqrt{(\ell +2)!/(\ell -2)!}$ (see Sect. 3.1.2 for details), which boosts the power on small scales where noise is dominant. This reduces our capability to trace and remove Galactic contamination at low multipoles. To solve this problem, we have implement the HZB-NILC, where the E–B leakage effect is corrected with the ZB method, while the NILC weights are estimated on the B-mode maps corrected with the recycling technique.

We note that the upper limit on the tensor-to-scalar ratio obtained with the HZB-NILC pipeline and assuming the Galactic model d0s0 is sensibly lower than the one obtained in Addamo et al. (2021). This difference can be motivated by taking into account several factors: (i) in this work, we neither consider the scanning strategy and thus nor the anisotropic distribution of the instrumental noise, which may play a role in the performance of the component separation; (ii) we assume to be able to de-lens our B-mode power spectrum at a 50% level (Yu et al. 2017), thus reducing the impact of lensing B modes to the total uncertainty in the r reconstruction; and (iii) for the simplistic d0s0 Galactic model, the minimisation of the variance in the needlet domain can in principle outperform parametric methods, which fit spectral parameters in different regions of the sky.

In this paper, we have presented the development and validation of a pipeline, based on the NILC component separation method, for the analysis of future CMB B-mode data collected by experiments that will observe the microwave sky from the ground and balloons. Taking into account real-world issues due to the incomplete sky coverage and state-of-the-art simulations of Galactic foregrounds, our results demonstrate the effectiveness of NILC, which emerges as a valid alternative to parametric component separation techniques for these kinds of experiments. Furthermore, this study permits one to easily extend the application of other ILC based techniques, such as cMILC (Remazeilles et al. 2021) and MC-NILC (Carones et al. 2022), to multi-frequency B-mode data from ground-based and balloon-borne experiments.

Acknowledgments

We thank Luca Pagano, Francesco Piacentini, and Giuseppe Puglisi for helpful comments. MM and NV acknowledge support by ASI/COSMOS grant n. 2016-24-H.0 and ASI/LiteBIRD grant n. 2020-9-HH.0. This work has been supported by the ASI-INFN Grant Agreement N. 2021-43-HH.0. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory. Part of this work was also supported by the InDark and LiteBIRD INFN projects. D.M. acknowledges support from the MIUR Excellence Project awarded to the Department of Mathematics, Università di Roma Tor Vergata, CUP E83C18000100006.

References

Appendix A: The NILC bias

The NILC weights are estimated employing the empirical covariance matrix of the input multi-frequency channels (see Eq. 8). If the covariance matrix is not correctly estimated due to the finite size of the domain over which it is computed, empirical correlations between the CMB modes and the residual contaminants are generated, leading to the de-projection of some of the CMB signal together with foregrounds and noise. This phenomenon is known as NILC bias and can cause a negative bias in the reconstructed CMB angular power spectrum, especially on the largest scales.

To avoid such an effect, enough modes have to be sampled in the computation of the covariance matrix to preserve the input CMB signal. In general, this problem is tackled by identifying a proper domain to compute the sample average of the product of needlet coefficients. Furthermore, usually the first needlet band is built from the sum in quadrature of several harmonic filters at low multipoles (see Eq. 24).

In practise, after the application of a component separation method on the B-mode data, one would like to extract cosmological information from the angular power spectrum ($C_{\ell }^{\mathit{out}}$) of the BB cleaned map. This quantity is composed of several terms:

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{out}}=C_{\ell }^{\mathit{cmb}}+C_{\ell }^{\mathit{fgds}}+C_{\ell }^{\mathit{noi}}+2·C_{\ell }^{c-f}+2·C_{\ell }^{c-n}+2·C_{\ell }^{n-f},\end{array}$$(A.1)

where $C_{\ell }^{\mathit{cmb}},C_{\ell }^{\mathit{fgds}}$, and $C_{\ell }^{\mathit{noi}}$ are the angular power spectra of CMB, residuals of foregrounds and noise, while the other terms represent the corresponding correlations among these components.

With de-noising techniques (cross-spectra of half-mission solutions or Monte Carlo simulations) one can eliminate $C_{\ell }^{\mathit{noi}}$, while the foregrounds contamination ($C_{\ell }^{\mathit{fgds}}$) can be marginalised at the likelihood level. Therefore, our final estimate of the angular power spectrum of the NILC CMB solution reads:

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }^{\mathit{out}}=C_{\ell }^{\mathit{cmb}}+2·C_{\ell }^{c-f}+2·C_{\ell }^{c-n}+2·C_{\ell }^{n-f}.\end{array}$$(A.2)

If NILC is properly implemented, the correlation terms should be very low and the CMB well reconstructed.

In this paper, we analyse simulated data. Thus, it is possible to assess the goodness of the CMB angular power spectrum reconstruction by computing:

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{bias}}=C_{\ell }^{\mathit{out}}-C_{\ell }^{\mathit{fgds}}-C_{\ell }^{\mathit{noi}}-C_{\ell }^{\mathit{cmb}},\end{array}$$(A.3)

which should be compatible with zero at all angular scales.

The relative bias $C_{\ell }^{\mathit{bias}}/C_{\ell }^{\mathit{cmb}}$ is shown in Figs. A.1 and A.2 for the different cases considered in Sect. 5: the application of r-NILC to the LSPE+Planck dataset assuming the d0s0 or d1s1 sky model (see the top panels of Figs. 13 and 14) and to SO-SAT dataset with the d1s1 Galactic model (see the top panels of Fig. 15); the application of ZB- and Hybrid ZB-NILC in the same cases (see the bottom panels of Figs. 13, 14, and 15). We can observe that the bias is always fully compatible with zero given the uncertainty on the reconstructed mean $C_{\ell }^{\mathit{bias}}$, which is estimated by dividing the dispersion of the angular power spectra of the bias among the different CMB reconstructions by the square root of the number of simulations: $\sigma (C_{\ell }^{\mathit{bias}})/\sqrt{N_{\mathit{sims}}}$.

Fig. A.1. Bias in the CMB B-mode power spectrum reconstruction for SWIPE+Planck dataset. The ratio between the bias $C_{\ell }^{\mathit{bias}}$ (see Eq. A.3) and the input CMB power spectrum $C_{\ell }^{\mathit{cmb}}$. All angular power spectra were computed as the average among 200 different simulated NILC CMB reconstructions and adopting the masks employed in Sect. 5. The shaded regions display the relative uncertainty computed as $\sigma (C_{\ell }^{\mathit{bias}})/(C_{\ell }^{\mathit{cmb}}·\sqrt{N_{\mathit{sims}}})$, where Nsims is the number of simulations. Top: Results of the application of rit-NILC (blue) and rin-NILC (red). Bottom: Bias from ZB (blue) and hybrid ZB-NILC (red). The Galactic emission was simulated with the d0s0 (left) and the d1s1 (right) sky models.

Fig. A.2. Bias in the CMB B-mode power spectrum reconstruction for SO-SAT dataset. The ratio between the bias $C_{\ell }^{\mathit{bias}}$ (see Eq. A.3) and the input CMB power spectrum $C_{\ell }^{\mathit{cmb}}$. All angular power spectra are computed as the average among 200 different simulated NILC CMB reconstructions and adopting the masks employed in Sect. 5. The shaded regions display the relative uncertainty computed as $\sigma (C_{\ell })/(C_{\ell }^{\mathit{cmb}}·\sqrt{N_{\mathit{sims}}})$. Left: Results of the application of r-NILC. Right: Bias from the ZB (blue) and hybrid ZB-NILC (red) methods. In all cases input foregrounds emission was simulated assuming the d1s1 sky model.

All Tables

All Figures

	Fig. 1. Sky coverage of the LSPE-SWIPE (left) and SO-SAT (right) instruments in Galactic coordinates. The observed region of the sky is shown in grey.
In the text

	Fig. 2. Sky area with the highest signal-to-noise ratio in the SO-SAT patch. This was obtained considering pixels with the largest values in the map of hit counts for a final sky fraction of fsky = 10%. The ILC methods were applied within this patch.
In the text

	Fig. 3. Recycling method for E–B leakage correction.
In the text

Fig. 4. B-mode CMB power spectrum reconstruction after the E–B leakage correction with the application of the recycling method and the extensions presented in Sect. 3.1.1. On the left: Mean angular power spectrum over Nsims = 200 CMB-only simulations that include lensing and a primordial tensor signal with r = 0.01. Solid lines represent the power spectrum of exact B maps (input), which were reconstructed with a full-sky B decomposition of Q and U and then masked for power spectrum estimation; power spectra of leakage-corrected maps (output LC) are shown with circles and diamonds for SWIPE and SO-SAT, respectively; the dashed lines are the leakage residuals after the correction; and the dashed-dotted lines (if present) are the CMB B modes without any leakage correction (output noLC). On the right: Effective tensor-to-scalar ratio associated with the absolute difference between leakage-corrected and exact angular power spectra. The top panel shows the results when the recycling method was applied to CMB B-mode maps for the SWIPE (red) and SO-SAT (blue) footprints; the middle panel reports those obtained when an iterative B decomposition with no (red), one (green), two (cyan), or three (magenta) iterations were performed. In the bottom panel, we compare the performance of iterative decomposition with three iterations (magenta) and diffusive inpainting (orange) when correcting the residual leakage. The middle and bottom panels present results only for SWIPE. The adopted binning scheme is Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 for the spectra computed in the LSPE-SWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The error bars highlight the uncertainty on the mean at 1σ, which was estimated from the dispersion of the angular power spectra of the simulations divided by $\sqrt{N_{\mathrm{sims}}}$. Readers can refer to the main text for details.

In the text

	Fig. 5. Iterative B decomposition extension of the recycling E–B leakage correction method.
In the text

Fig. 6. LSPE-SWIPE only. Ratio between the mean angular power spectrum of either pure B modes only (solid lines) or E–B leakage residuals (dashed lines) after having applied the iterative B decomposition correction and that of corresponding input maps. The cases with 1 (red), 2 (blue), 4 (magenta), and 9 (green) iterations are shown. 200 different realisations of CMB have been considered, including lensing and a primordial tensor signal with r = 0.01. The trend of dashed lines highlights the reduction of E–B leakage residuals obtained by performing the iterations. The trend of solid lines, instead, the loss of power suffered due to the B–E leakage when applying the iterative B decomposition. A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed. The error-bars highlight the uncertainty on the mean at 1σ, estimated from the dispersion of the angular power spectra of the simulations divided by $\sqrt{N_{\mathrm{sims}}}$. Readers can refer to the main text for details.

In the text

Fig. 7. LSPE-SWIPE only. Fitted tensor-to-scalar ratio on the mean angular power spectrum over 200 CMB-only simulations, including lensing and different values of the amplitude of the primordial tensor power spectrum. We consider the cases of B-mode maps reconstructed with: the recycling method (RM, red), the iterative recycling method with three iterations (magenta), the inpainting recycling method (orange), and without any E–B leakage contamination (ideal, cyan). The input values of r are represented with vertical dashed grey lines. The used likelihood is shown in Eq. (15) and the covariance is given only by the input CMB signal (ideal case). A binning scheme of Δℓ = 15 is employed and the first three multipoles (ℓ ≤ 4) are excluded. The error bars indicate the bounds at 1σ obtained from the posterior distribution.

In the text

Fig. 8. B-mode CMB power spectrum reconstruction after the E–B leakage correction with the application of the ZB method. On the left: Mean angular power spectrum over 200 CMB-only simulations including lensing and a primordial tensor signal with r = 0.01. Solid lines represent the power of exact B-maps reconstructed with full-sky EB-decomposition of Q and U; circles and diamonds that of leakage-corrected maps (output LC) with the ZB method; the dashed lines the leakage residuals after the correction; the dashed-dotted lines the CMB B modes without any leakage correction (output noLC). The LSPE-SWIPE (red) and SO-SAT (blue) footprints have been considered. On the right: Effective tensor-to-scalar ratio associated with the absolute difference between leakage-corrected and exact angular power spectra is shown for SWIPE (red) and SO-SAT (blue). A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPE-SWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The error-bars highlight the uncertainty on the mean at 1σ. Readers can refer to the main text for details.

In the text

	Fig. 9. Needlet configuration employed in this work. Left panel: Mexican needlets bands with B = 1.5 in the harmonic domain. Right: Corresponding needlet filters in real space of the first five bands, plotted with the same colours of the associated harmonic bands.
In the text

Fig. 10. Goodness of the reconstruction of the CMB B-mode signal (lensing+r = 0.01) after needlet filtering (top) in the SWIPE (red) and SO (blue) patches or after smoothing CMB maps with input angular resolution of FWHM = 9′ to bring them to FWHM = 91′ in the SO footprint (bottom). On the left: Mean angular power spectra (over 200 CMB-only simulations) of the input exact B maps reconstructed with a full-sky EB-decomposition of Q and U (solid lines); those of the reconstructed maps after needlet filtering or beam convolution (dots and diamonds); and the corresponding residuals (dashed lines). The plot on the right represents the effective tensor-to-scalar ratios associated with the absolute differences between output power spectra after needlet filtering and beam convolution and exact input ones in the SWIPE (red) and SO (blue) patches. A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPE-SWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The error-bars highlight the uncertainty on the mean of the simulations at 1σ.

In the text

Fig. 11. Mean angular power spectra (over 200 simulations) of foregrounds (solid lines) and noise (dashed lines) residuals when PILC (green), ILC (blue) or NILC (red) are applied on B-mode data (reconstructed on the full-sky and then masked) of SWIPE+Planck (left) or SO (right) frequency maps. The spectra have been computed considering the 15% (for SWIPE) and 7.5% (for SO) of the sky which is expected to be less contaminated by Galactic emission in the CMB frequency channels of SWIPE and SO. The input foregrounds emission is simulated with the d1s1 sky model. A binning scheme of Δℓ = 10 is employed.

In the text

Fig. 12. Trend of the tensor-to-scalar ratio rfgds fitted on the NILC foregrounds residuals angular power spectra with the likelihood of Eq. (26) when the sky fraction of the employed mask varies. The adopted masking strategy is the second one described in Sect. 4.1: we mask the most contaminated regions by foregrounds residuals. In this case, the input multi-frequency B-mode maps are constructed full-sky and then masked. The considered cases are SWIPE+Planck (left) with both d0s0 (red) and d1s1 (green) Galactic models and SO-SAT dataset (right) assuming only the d1s1 sky model. The grey shaded areas show the targeted sensitivity of the LSPE and SO experiments. The angular power spectra in the likelihood are binned with Δℓ = 15. The error bars indicate the bounds at 2σ (left) and 1σ (right) obtained from the posterior distributions.

In the text

Fig. 13. NILC residuals for SWIPE+Planck dataset with the d0s0 Galactic model. On the left: Mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations. On the right: Posterior distributions of an effective tensor-to-scalar ratio fitted on the foregrounds residuals in the case of half de-lensing (AL = 0.5). For the estimation of the posterior a binning scheme of Δℓ = 15 has been used to make the angular power spectrum Gaussianly distributed (see Sect. 4.2). On the top: Blue and red lines represent the results when E − B leakage is corrected with, respectively, the iterative and inpainting recycling methods (rit-NILC and rin-NILC). On the bottom: Results when the ZB- (blue) and hybrid ZB-NILC (red) methods are applied. The adopted binning scheme is Δℓ = 6 for 5 ≤ ℓ ≤ 28 and Δℓ = 15 for ℓ ≥ 29 on the top left, while Δℓ = 6 for 2 ≤ ℓ ≤ 19 and Δℓ = 15 for ℓ ≥ 20 on the bottom left. In both cases the chosen binning strategy is just for visualisation purposes. Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.015, 0.03]. The lower bound represents the targeted upper limit at 95% CL by LSPE in case of no detection, the upper one the targeted detection at 99.7% CL. The spectra have been estimated without any masking of the regions most contaminated by foregrounds. The final fsky is 32% in the case on top, because 4% of the pixels closest to the borders are masked to avoid residual E–B leakage effects, and 28% in the case on the bottom due to the needed apodisation of the mask.

In the text

Fig. 14. NILC residuals for SWIPE+Planck dataset with the d1s1 Galactic model. On the left: Mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations. On the right: Posterior distributions of an effective tensor-to-scalar ratio fitted to the foregrounds residuals in the case of half de-lensing (AL = 0.5). For the estimation of the posteriors, a binning scheme of Δℓ = 15 has been used to make the angular power spectrum Gaussianly distributed (see Sect. 4.2). On the top: Blue and red lines represent the results when E–B leakage is corrected with, respectively, the iterative and inpainting recycling methods (rit-NILC and rin-NILC). On the bottom: Results when the ZB- (blue) and hybrid ZB-NILC (red) methods are applied. The adopted binning scheme is Δℓ = 6 for 5 ≤ ℓ ≤ 28 and Δℓ = 15 for ℓ ≥ 29 on the top left, while Δℓ = 6 for 2 ≤ ℓ ≤ 19 and Δℓ = 15 for ℓ ≥ 20 on the bottom left. In both cases the chosen binning strategy is just for visualisation purposes. Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.015, 0.03]. The lower bound represents the targeted upper limit at 95% CL by LSPE in case of no detection, the upper one the targeted detection at 99.7% CL. The spectra have been estimated masking the regions most contaminated by foregrounds for a final fsky of 12%.

In the text

Fig. 15. NILC residuals for SO-SAT dataset with the d1s1 Galactic model. On the left: Mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations. On the right: Posterior distributions of an effective tensor-to-scalar ratio fitted to the foregrounds residuals in the case of half de-lensing (AL = 0.5). Both for the plots of the angular power spectra and the estimation of the posterior a binning scheme of Δℓ = 15 has been used. On the top: Results when E–B leakage is corrected with the recycling method (r-NILC). On the bottom: Results when the ZB- (blue) and hybrid ZB-NILC (red) methods are applied. Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.003, 0.01]. The lower bound represents the targeted upper limit at 68% CL by SO in case of no detection, the upper one the targeted detection at few σ significance. The spectra have been estimated masking the regions most contaminated by foregrounds for a final fsky of 8%.

In the text

Fig. A.1. Bias in the CMB B-mode power spectrum reconstruction for SWIPE+Planck dataset. The ratio between the bias $C_{\ell }^{\mathit{bias}}$ (see Eq. A.3) and the input CMB power spectrum $C_{\ell }^{\mathit{cmb}}$. All angular power spectra were computed as the average among 200 different simulated NILC CMB reconstructions and adopting the masks employed in Sect. 5. The shaded regions display the relative uncertainty computed as $\sigma (C_{\ell }^{\mathit{bias}})/(C_{\ell }^{\mathit{cmb}}·\sqrt{N_{\mathit{sims}}})$, where Nsims is the number of simulations. Top: Results of the application of rit-NILC (blue) and rin-NILC (red). Bottom: Bias from ZB (blue) and hybrid ZB-NILC (red). The Galactic emission was simulated with the d0s0 (left) and the d1s1 (right) sky models.

In the text

Fig. A.2. Bias in the CMB B-mode power spectrum reconstruction for SO-SAT dataset. The ratio between the bias $C_{\ell }^{\mathit{bias}}$ (see Eq. A.3) and the input CMB power spectrum $C_{\ell }^{\mathit{cmb}}$. All angular power spectra are computed as the average among 200 different simulated NILC CMB reconstructions and adopting the masks employed in Sect. 5. The shaded regions display the relative uncertainty computed as $\sigma (C_{\ell })/(C_{\ell }^{\mathit{cmb}}·\sqrt{N_{\mathit{sims}}})$. Left: Results of the application of r-NILC. Right: Bias from the ZB (blue) and hybrid ZB-NILC (red) methods. In all cases input foregrounds emission was simulated assuming the d1s1 sky model.

In the text