Issue 
A&A
Volume 677, September 2023



Article Number  A147  
Number of page(s)  23  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202244824  
Published online  18 September 2023 
Analysis of Needlet Internal Linear Combination performance on Bmode data from suborbital experiments
^{1}
Dipartimento di Fisica, Università di Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Roma, Italy
^{2}
Sezione INFN Roma 2, Via della Ricerca Scientifica 1, 00133 Roma, Italy
email: alessandro.carones@roma2.infn.it
^{3}
Dipartimento di Matematica, Universita’ di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Roma, Italy
Received:
28
August
2022
Accepted:
25
June
2023
Context. The observation of primordial B modes in cosmic microwave background (CMB) polarisation data represents the main scientific goal of most of the future CMB experiments. This signal is predicted to be much lower than polarised Galactic emission (foregrounds) in any region of the sky, pointing to the need for effective component separation methods.
Aims. Among all the techniques, the blind Needlet Internal Linear Combination (NILC) is of great relevance given our current limited knowledge of the Bmode foregrounds. In this work, we explore the possibility of employing NILC for the analysis of B modes reconstructed from partialsky data, specifically addressing the complications that such an application yields such as E–B leakage, needlet filtering, and beam convolution.
Methods. We consider two complementary simulated datasets of future experiments: the balloonborne Short Wavelength Instrument for the Polarisation Explorer (SWIPE) of the Large Scale Polarisation Explorer, which targets the observation of both reionisation and recombination peaks of the primordial CMB Bmode angular power spectrum, and the groundbased Small Aperture Telescope of Simons Observatory, which, instead, is designed to observe only the recombination bump at ℓ ∼ 80. We assessed the performance of the following two alternative techniques to correct for the CMB E–B leakage: the recycling technique and the ZhaoBaskaran method.
Results. We find that both techniques reduce the E–B leakage residuals at a negligible level given the sensitivity of the considered experiments, except for the recycling method in the SWIPE footprint at ℓ < 20. Thus, we implemented two extensions of the pipeline, the iterative B decomposition and the diffusive inpainting, which enabled us to recover the input CMB Bmode power for ℓ ≥ 5. For the considered experiments, we demonstrate that needlet filtering and beam convolution do not affect the CMB Bmode reconstruction. Finally, with an appropriate masking strategy, we find that NILC foregrounds subtraction allows one to achieve sensitivities on the tensortoscalar ratio in agreement with the targets of the considered CMB experiments.
Key words: cosmic background radiation / cosmology: observations / methods: data analysis
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 socalled 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 largescale 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 signaltonoise 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 tensortoscalar ratio r, which is defined as the amplitude of primordial gravitational waves over the one of initial density perturbations. Such a primordial Bmode 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 Bmode polarisation, either from the ground, such as for POLARBEAR (Arnold et al. 2010), QUBIC (Qubic & Battistelli 2011), BICEP (Wu et al. 2016), KeckArray (Keating et al. 2003), LSPESTRIP (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 LSPESWIPE (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 templatefitting 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 multifrequency 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, templatefitting algorithms, such as SEVEM (MartínezGonzález et al. 2003), try to construct internal foreground templates either in pixel or harmonic space using multifrequency 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 socalled 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 mismodelling 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 fullsky satellite data; however, in the near future, new observations of the CMB polarisation will be obtained mainly from groundbased and balloonborne 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 partialsky CMB Bmode 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 balloonborne 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 groundbased 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 balloonbased 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 tensortoscalar 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 6m primary mirror and three 0.5m 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 Bmode 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 multifrequency datasets. In Sect. 3 we describe the ILC methods considered in this work and how to extend the application of the NILC algorithm to partialsky Bmode 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 LSPESWIPE 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: LSPESWIPE and the SOSAT. We do not consider SOLAT, because it will not be devoted to primordial Bmode science. In the case of LSPESWIPE, 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 Bmode 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′.
Instrumental properties of the different CMB experiments considered in this work.
The Q and U maps at the different frequencies are obtained from the coaddition of three separate components: CMB, Galactic emission (synchrotron and thermal dust), and Gaussian white and isotropic noise. We do not include polarised extraGalactic sources, since they are expected to be dominant on angular scales smaller than those of main interest for inflationary Bmode science (ℓ > 100; Puglisi et al. 2018).
The CMB maps are generated with the HEALPix Python package^{1} as random Gaussian realisations from the angular power spectra of the bestfit Planck2018 parameters (Planck Collaboration VI 2020) with tensortoscalar 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^{2} (Thorne et al. 2017). Synchrotron polarised emission is modelled with a simple power law (Rybicki & Lightman 1979):
with X = {Q, U}, the position in the sky and ν the frequency. represents the synchrotron template at a reference frequency ν_{0}. Thermal dust emission is modelled with a modified blackbody (MBB):
where B_{ν}(T) is the blackbody 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 Planck2015 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 9year 23GHz 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 23GHz 7year data (MivilleDeschê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 Bmode 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 LSPESWIPE and SOSAT are shown in Fig. 1. The SOSAT 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 SOSAT 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 LSPESWIPE (left) and SOSAT (right) instruments in Galactic coordinates. The observed region of the sky is shown in grey. 
Fig. 2. Sky area with the highest signaltonoise ratio in the SOSAT patch. This was obtained considering pixels with the largest values in the map of hit counts for a final sky fraction of f_{sky} = 10%. The ILC methods were applied within this patch. 
3. Extension of ILC methods to partialsky 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ándezCobos et al. 2016),

Needlet ILC (NILC; Delabrouille et al. 2009).
All the above ILC techniques consist in linearly combining the N_{ν} multifrequency maps of one or multiple CMB experiments with frequencydependent weights ω:
where X is the Bfield for ILC, P = Q + iU for PILC, and a set of Bmode needlet coefficients for NILC (as discussed in the text below).
The weights are derived by minimising the variance of 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 deproject 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}:
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 . 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:
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 multifrequency 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ándezCobos et al. 2016; Delabrouille et al. 2009). The only assumption of these methods is that the CMB has a known emission law (blackbody 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 Bmode map at frequency i, , are obtained by filtering its harmonic coefficients with a weighting function b_{j}(ℓ) that selects modes at different angular scales for each needlet scale j:
where 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 at each scale j:
The pixeldependent weights are computed as follows:
where the covariance matrix for scale j at pixel 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 , whose width varies with the considered needlet scale j.
The final NILC map is then reconstructed by filtering again the harmonic coefficients in Eq. (7) and summing them all for each ℓ and m:
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 subspace. 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 groundbased and balloonborne experiments, ILC methods cannot be applied straightforwardly to Bmode 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 partialsky Q and U maps (see Sect. 3.1). It affects the ILC and NILC pipelines;

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

degrading cutsky 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 partialsky Stokes parameters is not exact: some modes (ambiguous modes) satisfy both the E and Bconditions simultaneously. When we split the polarisation field into an E and a Bpart, 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 overestimation 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 Bmode 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 leakagecorrected, needletfiltered, and smoothed CMB simulations with those of the input exact Bmode signal, reconstructed with a fullsky 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 pseudoangular power spectrum, which indicates the rotationally invariant variance of the harmonic coefficients of a map:
It is customary to report an analogous quantity: , which would be constant on large scales for a scaleinvariant primordial density perturbations spectrum. In this work, the estimation of the power spectra of the Bmode maps is performed with the NaMaster Python package^{3} 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 spin2 quantity that can be described in terms of the Stokes parameters Q and U as follows:
where denotes the position in the sky. This field can be expanded in spinweighted spherical harmonics:
From the harmonic coefficients in Eq. (12), it is possible to build the polarisation fields E (a scalar map of evenparity) and B (a pseudoscalar field of oddparity):
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 partialsky 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 cutsky is indeed not unique due to the nonlocality 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 Emode 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 partialsky observations, leading to an overestimate of the power in the reconstructed CMB Bmode 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 CMBonly simulations that include lensing from E modes and primordial tensor perturbations with r = 0.01 for the cases of the SWIPE and SOSAT patches (shown in Fig. 1). We have chosen this value of the tensortoscalar 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, and , of leakagecorrected and exact Bmode maps. The exact B modes are reconstructed with a fullsky B decomposition of CMB Q and U simulations. Furthermore, we consider an effective tensortoscalar ratio:
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 socalled E and Bfamily: P_{E} = (Q_{E}, U_{E}) and P_{B} = (Q_{B}, U_{B}), which receive contributions only from E and Bmode harmonic coefficients, respectively. For partialsky 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, 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 leakagecorrected Bfamily . The fit is performed with a simple leastsquares linear regression method. We then obtain the final CMB Bmode map through a B decomposition of . 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 CMBonly simulations with lensing+r = 0.01 for the SWIPE and SOSAT 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 fullsky 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. Bmode 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 N_{sims} = 200 CMBonly 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 fullsky B decomposition of Q and U and then masked for power spectrum estimation; power spectra of leakagecorrected maps (output LC) are shown with circles and diamonds for SWIPE and SOSAT, respectively; the dashed lines are the leakage residuals after the correction; and the dasheddotted lines (if present) are the CMB B modes without any leakage correction (output noLC). On the right: Effective tensortoscalar ratio associated with the absolute difference between leakagecorrected and exact angular power spectra. The top panel shows the results when the recycling method was applied to CMB Bmode maps for the SWIPE (red) and SOSAT (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 LSPESWIPE 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 . Readers can refer to the main text for details. 
The leakage correction is performed on the full SWIPE and SOSAT 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 Bmode power spectrum. It is possible to observe that the residuals after the application of the recycling method are negligible for the SOSAT 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 tensortoscalar ratio associated with the leakage residuals is at the level of r_{leak} ∼ 10^{−4} 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^{−4}, 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 Bmode 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 SOSAT 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 Bfamily, , by performing a B decomposition of the masked . The starting set of Stokes parameters for the iterations, , is the leakagecorrected Bfamily returned by the recycling method, . With such iterations, the pure B modes in the map are preserved, while the residual ambiguous E modes are progressively deprojected, 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 CMBonly simulations when either the recycling method alone (iter 0) or a postprocessing with 1, 2 or 3 iterative B decompositions is applied. The plots highlight how, with this additional procedure, we are able to deproject the residual E–B leakage contamination in the Bmode 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^{−4}) 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 partialsky 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 fullsky analysis, are erroneously estimated as E modes, causing a loss of power in the final Bmode 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 Bmode signal, possibly leading to an underestimate of the input tensortoscalar 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 tensortoscalar ratio fitted on the obtained residuals (see the right panels of Fig. 4) highlights that the first three multipoles (ℓ ≤ 4) suffer a significant BE leakage in the SWIPE case (r_{leak} ∼ 7 × 10^{−3}) 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 underestimate 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 Emode contamination. To assess whether this is the case, we perform the iterative B decomposition on two separate masked initial Bfamilies: and . The former is the set of Stokes parameters reconstructed with a fullsky B decomposition of input Q and U and represents the ideal goal of a leakage correction method. , 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 Bmode maps from the B decomposition performed on the two families after several iterations with those of the initial and . It is possible to observe that for ℓ ≥ 5, iterations do not deproject modes of the B maps obtained from the initial 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 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. LSPESWIPE 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 errorbars highlight the uncertainty on the mean at 1σ, estimated from the dispersion of the angular power spectra of the simulations divided by . Readers can refer to the main text for details. 
Given the previous results, we apply the NILC pipeline on multifrequency Bmode maps reconstructed in the SWIPE patch and processed with the recycling method and three iterations of B decomposition. Hereafter, we refer to it as ritNILC.
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 biLaplacian equation:
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 biLaplacian with the Laplacian and neglecting the Neumann boundary conditions. In this case, a template of residual ambiguous modes in the recycling CMB Bmap 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 CMBonly simulations with lensing+r = 0.01 on the LSPESWIPE 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 underestimate 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 Bmaps reconstructed in the SWIPE patch and corrected with the recycling method and the inpainting technique. Hereafter we refer to it as rinNILC.
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 Bmode CMB anisotropies for ℓ ≥ 5. We generate CMB simulations with different input amplitudes of tensor perturbations and obtain leakagecorrected Bmode maps with the three techniques described above: the standard recycling method, iterative recycling, and inpainting recycling. Then, a tensortoscalar ratio is fitted on the average binned angular power spectrum () of the output CMB Bmode maps for all different cases and excluding the first three multipoles (ℓ ≤ 4). The fit of the tensortoscalar ratio is performed with a Gaussian likelihood (Hamimeche & Lewis 2008; Gerbino et al. 2020):
Fig. 7. LSPESWIPE only. Fitted tensortoscalar ratio on the mean angular power spectrum over 200 CMBonly simulations, including lensing and different values of the amplitude of the primordial tensor power spectrum. We consider the cases of Bmode 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 is the binned angular power spectrum of primordial tensor CMB B modes for a tensortoscalar ratio r = 1, the binned BB lensing angular power spectrum and the covariance matrix associated with . 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 tensortoscalar 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 pseudoscalar polarisation fields are employed:
where and ð_{s} are the spinlowering and raising operators:
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_{ℓ}:
In the case of partialsky observations, the spin operators have to be applied on P_{+} and P_{−} fields multiplied by a window function , (nonzero only in the observed patch of the sky), and the recovered quantities and will be affected by the E–B leakage:
In Zhao & Baskaran (2010), the correction term to reconstruct the pure pseudoscalar field ℬ in the observed region of the sky has been derived as:
where:
with , and .
In this framework, as window functions, we adopt the footprints of SWIPE and SO (shown in Fig. 1) apodised with the ‘C^{1}’ scheme, where any observed pixel is multiplied by a weighting function f:
with , 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 illbehaved nature of W^{−1} and W^{−2} 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 and 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 leakagecorrected CMBonly simulations with lensing+r = 0.01 for the SWIPE and SO patches with that of the exact signal reconstructed with a fullsky 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 tensortoscalar ratio: 10^{−7} ≲ r_{leak} ≲ 10^{−3} (see the right panel in Fig. 8). With the ZB method we are even able to confidently reconstruct the primordial tensor Bmode 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. Bmode 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 CMBonly simulations including lensing and a primordial tensor signal with r = 0.01. Solid lines represent the power of exact Bmaps reconstructed with fullsky EBdecomposition of Q and U; circles and diamonds that of leakagecorrected maps (output LC) with the ZB method; the dashed lines the leakage residuals after the correction; the dasheddotted lines the CMB B modes without any leakage correction (output noLC). The LSPESWIPE (red) and SOSAT (blue) footprints have been considered. On the right: Effective tensortoscalar ratio associated with the absolute difference between leakagecorrected and exact angular power spectra is shown for SWIPE (red) and SOSAT (blue). A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPESWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The errorbars highlight the uncertainty on the mean at 1σ. Readers can refer to the main text for details. 
3.2. Needlet filtering of partialsky maps
The second step in the preprocessing of Bmode data which could introduce distortions in the CMB reconstruction before applying a NILC minimisation is the needlet filtering of partialsky observations. Needlet coefficients, indeed, are computed by convolving a map with a filter in pixelspace. 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 Bmode maps from a fullsky E–B decomposition of Q and U, simulated with Planck 2018 bestfit parameters and lensing+r = 0.01. Then, we apply a forward and an inverse needlet transformation to these Bmode 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, , and those of the input ones, . In analogy to the analysis of E − B leakage correction, we consider an effective tensortoscalar ratio r_{nl}:
which quantifies the absolute error we can make in the partialsky 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 Gaussianrelated weighting function:
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:
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 Bmode 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^{−4} 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 Bmode 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 CMBonly simulations) of the input exact B maps reconstructed with a fullsky EBdecomposition 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 tensortoscalar 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 LSPESWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The errorbars 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: LSPESWIPE complemented with Planck frequencies (SWIPE+PLANCK) and SOSAT.
For SWIPE+Planck it suffices to bring the fullsky Planck maps to the lower SWIPE angular resolution of FWHM = 85′.
For SOSAT, instead, where frequency channels have different angular resolutions, such an operation is not trivial. The convolution for a beam of partialsky 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 Bmode 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 SOSAT patch, of these maps when they are brought to the lowest resolution of SO (FWHM = 91′) fullsky () or after having been masked (). The comparison between these spectra is shown in the bottom panels of Fig. 10 together with an assessment of the residuals through an effective tensortoscalar ratio:
It is possible to observe that the most aggressive convolution, which has to be applied to SO maps, does not considerably affect the Bmode CMB reconstruction with values of r_{smooth} of the order of or below 10^{−4} at all the angular scales of interest for Bmode science.
4. NILC performance: Analysis setup
After having applied a leakage correction method and needlet filtering to input partialsky Bmode multifrequency maps, the NILC application leads to the best blind estimate of the CMB Bmode 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 tensortoscalar 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:

computing the polarised intensity P^{2} = Q^{2} + U^{2} 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};

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 SOSAT 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 rescaling 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 smallscale 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 methodspecific. 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 illbehaved estimation of the power spectrum. We adopt the C^{1} 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 SOSAT) that we consider for the computation of the angular power spectra.
In the case of the NILC application on SWIPE+PlanckBmode 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 SOSAT, 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 Bmode 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 () can be written as follows:
where and are the angular power spectra of residual foregrounds and noise. The noise bias, , can be subtracted at the power spectrum level by estimating its contribution through Monte Carlo simulations or by computing crossspectra between uncorrelated splits of the data (e.g. maps from subsets of detectors or from partialmission 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 , we fitted r on binned assuming a Gaussian likelihood (Hamimeche & Lewis 2008; Gerbino et al. 2020):
where is the primordial tensor CMB Bmode angular power spectrum for a tensortoscalar ratio r = 1 and 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 inverseWishart likelihood function.
The covariance matrix is estimated from the angular power spectra of 200 NILC simulated solutions as:
where A_{L} quantifies our ability to delens the output map with A_{L} = 1 for no delensing and A_{L} = 0 for full delensing. 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 crossterms. In this work, we assume A_{L} = 0.5 (which is the level that can be achieved with multitracers techniques; see e.g. Yu et al. 2017).
5. Results
In the processing of future groundbased and balloonborne 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 Bmode signal in an ideal case where we assume to be able to perfectly reconstruct the needlet and pixel coefficients of Bmode data even for partialsky 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 leakagecorrected with the recycling or ZB method. In this case, also the needlet filtering and the beam convolution are performed on partialsky maps.
In all the following analyses, the input CMB Bmode 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 SOSAT.
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 SOSAT 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 partialsky B modes data: Ideal case
At the moment, PILC (see Sect. 3) is the only blind method available to easily subtract foregrounds emission from groundbased and balloonborne 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 Bmode signal in needlet and pixel space for partialsky observations. This is explicitly obtained by reconstructing from simulated fullsky Q and U maps fullsky Bmode 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 fullsky reconstruction of Bmode 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 Bmode maps from partialsky 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 Bmode data (reconstructed on the fullsky 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 (due to the presence of the Emode signal) than in B and therefore the PILC weights are less capable of tracing and subtracting the frequencydependent contaminants in the Bmode 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 cutsky observations given the sensitivity to the tensortoscalar 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 Bmode data for groundbased and balloonborne experiments.
5.2. NILC results for realistic SWIPE and SO case studies
A realistic application of NILC on partialsky Bmode maps can be summarised with the following steps:

Correct the E–B leakage effect in the multifrequency 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:

rNILC, for maps corrected with the recycling method,

ritNILC, when the recycling method and iterative B decomposition are applied,

rinNILC, for maps corrected with the recycling method and diffusive inpainting,

ZBNILC, when the ZB technique is performed.
We have verified that all the leakagecorrection 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 SOSAT, 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 ritNILC, rinNILC, and ZBNILC to the SWIPE+Planck dataset, while rNILC and ZBNILC to SOSAT simulated maps. In the case of the application of ritNILC and rinNILC 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 SOSAT 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 tensortoscalar ratio r_{fgds} 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 multifrequency Bmode maps are constructed fullsky and then masked. The considered cases are SWIPE+Planck (left) with both d0s0 (red) and d1s1 (green) Galactic models and SOSAT 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,

SOSAT 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 tensortoscalar 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^{−2} or greater. When, instead, NILC is applied on ideal simulated Bmode data of SOSAT, 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 tensortoscalar 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 (rNILC)
In this section we report the results obtained when the NILC algorithm is applied to multifrequency Bmode maps that include CMB, noise, and foregrounds, and have been corrected with the recycling method (rNILC). We note that, in the case of the LSPE+Planck dataset, the maps have also been postprocessed either with the iterative B decomposition with three iterations (ritNILC) or with the diffusive inpainting (rinNILC).
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 tensortoscalar ratio fitted on the foregrounds residuals in the case of half delensing (A_{L} = 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 (ritNILC and rinNILC). On the bottom: Results when the ZB (blue) and hybrid ZBNILC (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 f_{sky} 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 righthand side of Fig. 13 shows that the effective tensortoscalar ratio associated with foregrounds residuals has a posterior distribution with an upper limit at 95% CL of ∼1.5 × 10^{−2}, 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 Bmode 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 tensortoscalar ratio.
The angular power spectra of the residuals are shown in the upper left panel of Fig. 14. Both the ritNILC and rinNILC cleaning methods still allow one to possibly detect both peaks in the Bmode primordial spectrum. The posterior distribution of an effective tensortoscalar ratio fitted on the foregrounds residuals has an upper limit of around r ≲ 2.5 × 10^{−2} 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 tensortoscalar ratio fitted to the foregrounds residuals in the case of half delensing (A_{L} = 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 (ritNILC and rinNILC). On the bottom: Results when the ZB (blue) and hybrid ZBNILC (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 f_{sky} of 12%. 
The very same pipeline can be applied to SOSAT 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 SOSAT we have considered only a Galactic emission simulated with the d1s1 model. The residuals and the posterior distribution of the tensortoscalar 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 rNILC 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 partialsky extension of the NILC pipeline to Bmode 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 SOSAT 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 tensortoscalar ratio fitted to the foregrounds residuals in the case of half delensing (A_{L} = 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 (rNILC). On the bottom: Results when the ZB (blue) and hybrid ZBNILC (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 f_{sky} of 8%. 
5.2.2. ZBNILC
The ZBNILC pipeline consists of implementing the variance minimisation in needlet domain after having applied the ZB leakage correction method on multifrequency Bmode 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 Bmode signal is recovered only at the end when the angular power spectra are computed.
We consider the same datasets and cases of the rNILC 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 ZBNILC foregrounds residuals are much higher than those of the recycling NILC. This is caused by the fact that, in ZBNILC, the input maps are convolved for the harmonic filter (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, ZBNILC 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 ZBcorrected maps, we have implemented the Hybrid ZBNILC (HZBNILC). HZBNILC 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 Bmode maps processed with the standard recycling method.
The foregrounds and noise residuals of the HZBNILC 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 tensortoscalar ratio fitted on the foregrounds residuals of 5 × 10^{−3} 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 rNILC is motivated by the inclusion in this analysis of the lowest multipoles ℓ ≤ 4.
The upper limit on the tensortoscalar ratio obtained with HZBNILC 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 delens our Bmode 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 ZBNILC 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 HZBNILC). 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 tensortoscalar 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 rNILC pipelines is given by the possibility of including the first three multipoles in the likelihood.
The same trends are obtained when ZBNILC and Hybrid ZBNILC have been applied to the simulated SOSAT dataset with Galactic emission simulated with the d1s1 model (see the bottom panel of Fig. 15). The angular power spectrum of the HZBNILC 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 HZBNILC foregrounds residuals lead to an upper bound in the posterior distribution of the effective tensortoscalar ratio of r < 0.0028 (68% CL) which is compatible with the SO target (r ≤ 3 × 10^{−3} at 1σ).
We observe that for SOSAT the difference in the amplitudes of foregrounds residuals on large scales between ZB and HZBNILC is much less evident than the one for SWIPE+Planck. This is due to the fact that, without modes with ℓ ≤ 30, the calculation of HZBNILC weights in the first needlet band is more affected by noise contamination on intermediate scales. This condition is close to that experienced by the ZBNILC 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 Bmode data to the commonly used parametric component separation methods, which are more subject to systematic biases in the case of mismodelling of the foreground properties. However, most future surveys will be balloonborne or groundbased and will observe only a portion of the sky, whereas, at present, NILC has been successfully applied to fullsky data from either Planck or WMAP.
In this work, we explore the possibility to extend the NILC formalism to future Bmode partialsky 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: LSPESWIPE and SOSAT. 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 LSPESWIPE, 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 Bmode foregrounds.
When dealing with partialsky 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 Bmode angular power spectrum on Monte Carlo simulations, also quantifying a possible bias on the estimate of the tensortoscalar ratio. We find that:

The recycling method is able to reduce the E–B leakage contamination in the Bmode maps at a negligible level for SO at all the angular scales of interest (ℓ≥30), while for LSPESWIPE 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 postprocessing 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.
Capability to reconstruct the input angular power spectrum of CMB B modes with lensing and primordial tensor perturbations at different angular scales for the two experiments considered in this work (SWIPE and SOSAT).
Once the above preprocessing steps were tested, the performance of the NILC cleaning method was assessed on simulated multifrequency Bmode datasets which include CMB with only lensing, instrumental noise, and Galactic foreground emission. For LSPESWIPE 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 SOSAT, 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 partialsky Bmode 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 (ritNILC) for SWIPE+Planck,

recycling + diffusive inpainting (rinNILC) for SWIPE+Planck,

recycling (rNILC) for SOSAT, and

ZB (ZBNILC) for both SOSAT and SWIPE+Planck.
For all cases, we quantified the ability of the pipeline to mitigate foreground contamination in terms of an effective tensortoscalar 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 LSPESWIPE and SO, respectively, as in the official forecasts of the two experiments (see Addamo et al. 2021 and Ade et al. 2019).
Upper bounds and constraints of an effective tensortoscalar ratio fitted on the average angular power spectrum of foreground residuals given by the application of NILC for the listed cases.
We find that SOSAT 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 SOSAT 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 LSPESWIPE.
When the ZBNILC pipeline is applied (see Sect. 5.2.2), the amplitude of the foreground residuals is much greater than that obtained with rNILC. This can be explained by the fact that in this case the input maps are convolved with the harmonic filter (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 HZBNILC, where the E–B leakage effect is corrected with the ZB method, while the NILC weights are estimated on the Bmode maps corrected with the recycling technique.
We note that the upper limit on the tensortoscalar ratio obtained with the HZBNILC 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 delens our Bmode 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 Bmode data collected by experiments that will observe the microwave sky from the ground and balloons. Taking into account realworld issues due to the incomplete sky coverage and stateoftheart 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 MCNILC (Carones et al. 2022), to multifrequency Bmode data from groundbased and balloonborne experiments.
Acknowledgments
We thank Luca Pagano, Francesco Piacentini, and Giuseppe Puglisi for helpful comments. MM and NV acknowledge support by ASI/COSMOS grant n. 201624H.0 and ASI/LiteBIRD grant n. 20209HH.0. This work has been supported by the ASIINFN Grant Agreement N. 202143HH.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
 Addamo, G., Ade, P. A. R., Baccigalupi, C., et al. 2021, J. Cosmol. Asytropart. Phys., 2021, 008 [Google Scholar]
 Ade, P. A. R., Akiba, Y., Anthony, A. E., et al. 2014, Phys. Rev. Lett., 113, 021301 [NASA ADS] [CrossRef] [Google Scholar]
 Ade, P. A. R., Ahmed, Z., Aikin, E. W., et al. (Keck Array and BICEP Collaboration) 2016, ApJ, 833, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmol. Asytropart. Phys., 2019, 056 [CrossRef] [Google Scholar]
 Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmol. Asytropart. Phys., 2020, 047 [CrossRef] [Google Scholar]
 Alonso, D., Sanchez, J., Slosar, A., & Dark, E. S. C. 2019, MNRAS, 484, 4127 [NASA ADS] [CrossRef] [Google Scholar]
 Arnold, K., Ade, P. A. R., Anthony, A. E., et al. 2010, Proc. SPIE, 7741, 77411E [NASA ADS] [CrossRef] [Google Scholar]
 Baldi, P., Kerkyacharian, G., Marinucci, D., & Picard, D. 2009, Ann. Stat., 37, 1150 [NASA ADS] [Google Scholar]
 Basak, S., & Delabrouille, J. 2012, MNRAS, 419, 1163 [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 97 [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 Bersanelli, M., Mennella, A., Morgante, G., et al. 2012, SPIE Conf. Ser., 8446, 7 [Google Scholar]
 Bunn, E. F., Zaldarriaga, M., Tegmark, M., & de OliveiraCosta, A. 2003, Phys. Rev. D, 67, 023501 [Google Scholar]
 Carones, A., Migliaccio, M., Puglisi, G., et al. 2022, arXiv eprints [arXiv:2212.04456] [Google Scholar]
 de Bernardis, P., Aiola, S., Amico, G., et al. 2012, SPIE Conf. Ser., 8452, 3 [Google Scholar]
 Delabrouille, J., & Cardoso, J. F. 2009, Lecture Notes Phys., 665, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J. F., & Patanchon, G. 2003, MNRAS, 346, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J.F., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dick, J., Remazeilles, M., & Delabrouille, J. 2010, MNRAS, 401, 1602 [CrossRef] [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [Google Scholar]
 FernándezCobos, R., MarcosCaballero, A., Vielva, P., MartínezGonzález, E., & Barreiro, R. B. 2016, MNRAS, 459, 441 [CrossRef] [Google Scholar]
 Filippini, J. P., Ade, P. A. R., Amiri, M., et al. 2010, Proc. SPIE, 7741, 77411N [NASA ADS] [CrossRef] [Google Scholar]
 Geller, D., & Mayeli, A. 2008, Math. Z., 262, 895 [Google Scholar]
 Gerbino, M., Lattanzi, M., Migliaccio, M., et al. 2020, Front. Phys., 8, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Ghosh, S., Delabrouille, J., Zhao, W., & Santos, L. 2021, J. Cosmol. Astropart. Phys., 2021, 036 [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [Google Scholar]
 Grain, J., Tristram, M., & Stompor, R. 2009, Phys. Rev. D, 79, 123515 [Google Scholar]
 Hamimeche, S., & Lewis, A. 2008, Phys. Rev. D, 77, 103013 [NASA ADS] [CrossRef] [Google Scholar]
 Hanany, S., Alvarez, M., Artis, E., et al. 2019, Bull. Am. Astron. Soc., 51, 194 [Google Scholar]
 Hanson, D., Hoover, S., Crites, A., et al. 2013, Phys. Rev. Lett., 111, 141301 [CrossRef] [Google Scholar]
 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [Google Scholar]
 Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. Lett., 78, 2058 [Google Scholar]
 Keating, B. G., Ade, P. A. R., Bock, J. J., et al. 2003, SPIE Conf. Ser., 4843, 284 [NASA ADS] [Google Scholar]
 Keisler, R., Hoover, S., Harrington, N., et al. 2015, ApJ, 807, 151 [Google Scholar]
 Kim, J. 2011, A&A, 531, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kim, J., & Naselsky, P. 2010, A&A, 519, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kovac, J. M., Leitch, E. M., Pryke, C., et al. 2002, Nature, 420, 772 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., Challinor, A., & Turok, N. 2001, Phys. Rev. D, 65, 023505 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, H., Creswell, J., von Hausegger, S., & Naselsky, P. 2019, Phys. Rev. D, 100, 023538 [Google Scholar]
 MacTavish, C. J., Ade, P. A. R., Bock, J. J., et al. 2006, ApJ, 647, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Maino, D., Farusi, A., Baccigalupi, C., et al. 2002, MNRAS, 334, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Marinucci, D., Pietrobon, D., Balbi, A., et al. 2008, MNRAS, 383, 539 [NASA ADS] [CrossRef] [Google Scholar]
 MartínezGonzález, E., Diego, J. M., Vielva, P., & Silk, J. 2003, MNRAS, 345, 1101 [Google Scholar]
 Matsumura, T., Akiba, Y., Borrill, J., et al. 2014, J. Low Temp. Phys., 176, 733 [Google Scholar]
 MivilleDeschênes, M. A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montroy, T. E., Ade, P. A. R., Bock, J. J., et al. 2006, ApJ, 647, 813 [NASA ADS] [CrossRef] [Google Scholar]
 Narcowich, F. J., Petrushev, P., & Ward, J. D. 2006, SIAM J. on Math. Anal., 38, 574 [CrossRef] [Google Scholar]
 Pietrobon, D., Balbi, A., & Marinucci, D. 2006, Phys. Rev. D, 74, 043524 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLI. 2016, A&A, 596, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 POLARBEAR Collaboration (Ade, P. A. R., et al.) 2017, ApJ, 848, 121 [Google Scholar]
 Puglisi, G., Galluzzi, V., Bonavera, L., et al. 2018, ApJ, 858, 85 [Google Scholar]
 Qubic, C., Battistelli, E., et al. 2011, Astropart. Phys., 34, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Rotti, A., & Chluba, J. 2021, MNRAS, 503, 2478 [NASA ADS] [CrossRef] [Google Scholar]
 Rybicki, D. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley) [Google Scholar]
 Sayre, J. T., Reichardt, C. L., Henning, J. W., et al. 2020, Phys. Rev. D, 101, 122003 [Google Scholar]
 Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175 [Google Scholar]
 Stolyarov, V., Hobson, M. P., Ashdown, M. A. J., & Lasenby, A. N. 2002, MNRAS, 336, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Stompor, R., Leach, S., Stivoli, F., & Baccigalupi, C. 2009, MNRAS, 392, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M., de OliveiraCosta, A., & Hamilton, A. J. 2003, Phys. Rev. D, 68, 123523 [Google Scholar]
 Thorne, B., Dunkley, J., Alonso, D., & Næss, S. 2017, MNRAS, 469, 2821 [Google Scholar]
 Tristram, M., Banday, A. J., Górski, K. M., et al. 2022, Phys. Rev. D, 105, 083524 [Google Scholar]
 Wu, W. L. K., Ade, P. A. R., Ahmed, Z., et al. 2016, J. Low Temp. Phys., 184, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, B., Hill, J. C., & Sherwin, B. D. 2017, Phys. Rev. D, 96, 123511 [NASA ADS] [CrossRef] [Google Scholar]
 Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830 [Google Scholar]
 Zaldarriaga, M., & Seljak, U. 1998, Phys. Rev. D, 58, 023003 [Google Scholar]
 Zhang, Z., Liu, Y., Li, S.Y., et al. 2022, J. Cosmol. Asytropart. Phys., 2022, 044 [CrossRef] [Google Scholar]
 Zhao, W., & Baskaran, D. 2010, Phys. Rev. D, 82, 023001 [Google Scholar]
Appendix A: The NILC bias
The NILC weights are estimated employing the empirical covariance matrix of the input multifrequency 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 deprojection 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 Bmode data, one would like to extract cosmological information from the angular power spectrum () of the BB cleaned map. This quantity is composed of several terms:
where , and are the angular power spectra of CMB, residuals of foregrounds and noise, while the other terms represent the corresponding correlations among these components.
With denoising techniques (crossspectra of halfmission solutions or Monte Carlo simulations) one can eliminate , while the foregrounds contamination () can be marginalised at the likelihood level. Therefore, our final estimate of the angular power spectrum of the NILC CMB solution reads:
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:
which should be compatible with zero at all angular scales.
The relative bias is shown in Figs. A.1 and A.2 for the different cases considered in Sect. 5: the application of rNILC to the LSPE+Planck dataset assuming the d0s0 or d1s1 sky model (see the top panels of Figs. 13 and 14) and to SOSAT dataset with the d1s1 Galactic model (see the top panels of Fig. 15); the application of ZB and Hybrid ZBNILC 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 , 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: .
Fig. A.1. Bias in the CMB Bmode power spectrum reconstruction for SWIPE+Planck dataset. The ratio between the bias (see Eq. A.3) and the input CMB power spectrum . 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 , where N_{sims} is the number of simulations. Top: Results of the application of ritNILC (blue) and rinNILC (red). Bottom: Bias from ZB (blue) and hybrid ZBNILC (red). The Galactic emission was simulated with the d0s0 (left) and the d1s1 (right) sky models. 
Fig. A.2. Bias in the CMB Bmode power spectrum reconstruction for SOSAT dataset. The ratio between the bias (see Eq. A.3) and the input CMB power spectrum . 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 . Left: Results of the application of rNILC. Right: Bias from the ZB (blue) and hybrid ZBNILC (red) methods. In all cases input foregrounds emission was simulated assuming the d1s1 sky model. 
All Tables
Instrumental properties of the different CMB experiments considered in this work.
Capability to reconstruct the input angular power spectrum of CMB B modes with lensing and primordial tensor perturbations at different angular scales for the two experiments considered in this work (SWIPE and SOSAT).
Upper bounds and constraints of an effective tensortoscalar ratio fitted on the average angular power spectrum of foreground residuals given by the application of NILC for the listed cases.
All Figures
Fig. 1. Sky coverage of the LSPESWIPE (left) and SOSAT (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 signaltonoise ratio in the SOSAT patch. This was obtained considering pixels with the largest values in the map of hit counts for a final sky fraction of f_{sky} = 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. Bmode 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 N_{sims} = 200 CMBonly 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 fullsky B decomposition of Q and U and then masked for power spectrum estimation; power spectra of leakagecorrected maps (output LC) are shown with circles and diamonds for SWIPE and SOSAT, respectively; the dashed lines are the leakage residuals after the correction; and the dasheddotted lines (if present) are the CMB B modes without any leakage correction (output noLC). On the right: Effective tensortoscalar ratio associated with the absolute difference between leakagecorrected and exact angular power spectra. The top panel shows the results when the recycling method was applied to CMB Bmode maps for the SWIPE (red) and SOSAT (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 LSPESWIPE 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 . 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. LSPESWIPE 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 errorbars highlight the uncertainty on the mean at 1σ, estimated from the dispersion of the angular power spectra of the simulations divided by . Readers can refer to the main text for details. 

In the text 
Fig. 7. LSPESWIPE only. Fitted tensortoscalar ratio on the mean angular power spectrum over 200 CMBonly simulations, including lensing and different values of the amplitude of the primordial tensor power spectrum. We consider the cases of Bmode 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. Bmode 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 CMBonly simulations including lensing and a primordial tensor signal with r = 0.01. Solid lines represent the power of exact Bmaps reconstructed with fullsky EBdecomposition of Q and U; circles and diamonds that of leakagecorrected maps (output LC) with the ZB method; the dashed lines the leakage residuals after the correction; the dasheddotted lines the CMB B modes without any leakage correction (output noLC). The LSPESWIPE (red) and SOSAT (blue) footprints have been considered. On the right: Effective tensortoscalar ratio associated with the absolute difference between leakagecorrected and exact angular power spectra is shown for SWIPE (red) and SOSAT (blue). A binning scheme of Δℓ = 3 for ℓ ≤ 4, Δℓ = 6 for 5 ≤ ℓ ≤ 28, and Δℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPESWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The errorbars 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 Bmode 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 CMBonly simulations) of the input exact B maps reconstructed with a fullsky EBdecomposition 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 tensortoscalar 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 LSPESWIPE patch, while a constant Δℓ = 15 for those estimated in the SO footprint. The errorbars 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 Bmode data (reconstructed on the fullsky 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 tensortoscalar ratio r_{fgds} 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 multifrequency Bmode maps are constructed fullsky and then masked. The considered cases are SWIPE+Planck (left) with both d0s0 (red) and d1s1 (green) Galactic models and SOSAT 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 tensortoscalar ratio fitted on the foregrounds residuals in the case of half delensing (A_{L} = 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 (ritNILC and rinNILC). On the bottom: Results when the ZB (blue) and hybrid ZBNILC (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 f_{sky} 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 tensortoscalar ratio fitted to the foregrounds residuals in the case of half delensing (A_{L} = 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 (ritNILC and rinNILC). On the bottom: Results when the ZB (blue) and hybrid ZBNILC (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 f_{sky} of 12%. 

In the text 
Fig. 15. NILC residuals for SOSAT 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 tensortoscalar ratio fitted to the foregrounds residuals in the case of half delensing (A_{L} = 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 (rNILC). On the bottom: Results when the ZB (blue) and hybrid ZBNILC (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 f_{sky} of 8%. 

In the text 
Fig. A.1. Bias in the CMB Bmode power spectrum reconstruction for SWIPE+Planck dataset. The ratio between the bias (see Eq. A.3) and the input CMB power spectrum . 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 , where N_{sims} is the number of simulations. Top: Results of the application of ritNILC (blue) and rinNILC (red). Bottom: Bias from ZB (blue) and hybrid ZBNILC (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 Bmode power spectrum reconstruction for SOSAT dataset. The ratio between the bias (see Eq. A.3) and the input CMB power spectrum . 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 . Left: Results of the application of rNILC. Right: Bias from the ZB (blue) and hybrid ZBNILC (red) methods. In all cases input foregrounds emission was simulated assuming the d1s1 sky model. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.