Analysis of NILC performance on B-modes data from sub-orbital experiments

,


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 ΛCDM model, whose parameters are now tightly constrained (see, e.g., Boomerang: MacTavish et al. 2006;WMAP: Hinshaw et al. 2013;Planck: Planck Collaboration et al. 2020b).CMB polarisation data are usually decomposed into the so-called E and B modes, of even and odd parity, respectively, in the sky (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 et al. 2020b).B modes, instead, are thought to be associated to at least two independent mechanisms.On small angular scales (typically a few arcminutes), they are mainly sourced by primordial E modes distorted into B modes via weak gravitational lensing by the intervening large-scale structures along the line of sight from the last scattering surface to us.This signal is commonly known as lensing B modes Send offprint requests to: alessandro.carones@roma2.infn.it(Zaldarriaga & Seljak 1998) and has been observed by several experiments, although the detections are still with a modest signal-to-noise ratio or in very small patches of sky (Hanson et al. 2013;Ade et al. 2014;Planck Collaboration et al. 2016b;BICEP2 Collaboration et al. 2016;Keisler et al. 2015;POLARBEAR Collaboration et al. 2017).On large scales, B modes are expected to originate from tensor perturbations (primordial gravitational waves) generated in the very early Universe by a phase of cosmic inflation (Kamionkowski et al. 1997).Their magnitude is set by the tensor-to-scalar ratio r, which is defined as the amplitude of primordial gravitational waves over the one of initial density perturbations.Such a primordial B-mode signal still escapes detection so far with an upper limit of r ≤ 0.032 at 95 % CL (Tristram et al. 2022).The detection of the primordial gravitational waves background would represent a powerful test of cosmic inflation and 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 to the integral of the linear polarisation generated by quadrupolar anisotropies from the last scattering surface as seen by each free electron after reionisation, and the recombination bump on smaller scales (ℓ ∼ 80), which represents the imprint of the primordial gravitational waves on the physics of the last scattering surface.Many experiments have been designed or proposed to observe B-mode polarisation, either from the ground: POLARBEAR (Arnold et al. 2010), QUBIC (Qubic Collaboration et al. 2011), BICEP (Wu et al. 2016), Keck-Array (Keating et al. 2003), LSPE-STRIP (Bersanelli et al. 2012), ACT (Aiola et al. 2020), SPT (Sayre et al. 2020), Simons Observatory (Ade et al. 2019); from balloons: SPI-DER (Filippini et al. 2010), LSPE-SWIPE (de Bernardis et al. 2012); or from space: LiteBIRD (Matsumura et al. 2014), PICO (Hanany et al. 2019).The development of effective component separation methods has become one of the crucial aspects of data analysis for all these experiments.The reason being that the quest for primordial B modes is made much more difficult by the presence of instrumental noise and polarised foregrounds, especially Galactic thermal dust and synchrotron emission.These methods are usually divided into blind, parametric and template-fitting techniques.Blind methods, such as ILC (Bennett et al. 2003, Tegmark et al. 2003), NILC (Basak & Delabrouille 2012), FastICA (Maino et al. 2002), SMICA (Delabrouille et al. 2003), usually linearly combine multi-frequency maps in such a way as to minimise some meaningful statistical quantity of the final map.Parametric methods, such as Commander (Eriksen et al. 2008), FG-Buster (Stompor et al. 2009), FastMEM (Stolyarov et al. 2002), instead, explicitly model the frequency properties of the foregrounds by means of a set of parameters that are fitted to the data.Such methods provide an easy way to characterise and propagate foregrounds residual errors, but their effectiveness depends on how reliable the adopted model is.Finally, template-fitting algorithms, such as SEVEM (Martínez-González et al. 2003), try to construct internal foregrounds templates either in pixel or harmonic space using multi-frequency observations.There is no clear evidence (especially in polarisation) on which approach is more effective in reconstructing the CMB signal.Thus, applying several different methods on the same data-set allows performing comparisons and evaluating the robustness of the results.In this work, we focus our attention on the so-called Internal Linear Combination (ILC) algorithms.They are of great relevance in the context of CMB data-analysis because they perform noise and foregrounds 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.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.One of the drawbacks of such approaches is that the estimation of foregrounds residuals usually relies on Monte Carlo simulations or other bootstrapping techniques.Throughout the years, several extensions have been proposed; the Needlet ILC (NILC), which performs variance minimisation in needlet space, has proven to be one of the most promising blind techniques.NILC has been extensively applied to full-sky satellite data, but in the near future new observations of the CMB polarisation will be obtained mainly from ground-based and balloon-borne experiments, which will only observe a portion of the sky.Therefore, in this work, we explore the possibility of employing NILC for the analysis of partial-sky CMB B-mode data, specifically addressing the further challenges that such an extension yields: E-B leakage, needlet filter-ing, 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: the LSPE-SWIPE, a balloon-borne telescope that avoids most of the atmospheric contamination, whose scientific goal is the observation of both the reionisation and recombination peaks -The Simons Observatory, which is a ground-based experiment targeting the detection of the recombination bump.
The Short Wavelength Instrument for the Polarisation Explorer (SWIPE, de Bernardis et al. 2012) is one of the two instruments of the LSPE (Large Scale Polarisation Explorer) 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, expecting to constrain the tensor-to-scalar ratio down to r = 0.015 (at 95% CL).The SWIPE 145 GHz band will be the main channel for CMB science, while observations at 210 and 240 GHz will monitor thermal dust contamination.SWIPE 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) 85 ′ for around 15 days, exploiting the optimal observational conditions offered by the Arctic night (at latitude around 78 • N).
The Simons Observatory (SO) consists of a Large Aperture Telescope (LAT) with a 6-meter primary mirror and three 0.5-meter refracting Small Aperture Telescopes (SATs) (Ade et al. 2019).SO will be located in the Atacama Desert at an altitude of 5.200 metres in Chile's Parque Astronomico and will collect data in the frequency range between 25 and 280 GHz.It is designed to measure a signal at the r = 0.01 level at a few σ significance or to exclude it at similar significance, using the B-mode amplitude around the recombination bump.The reionisation peak, instead, cannot be observed by SO due to atmospheric contamination that will allow to sample only cosmological modes with ℓ > 30.
The paper is organised as follows.In Sect.2, we present the considered simulated multi-frequency data-sets; in Sect. 3 we describe the ILC methods considered in this work and how to extend the application of the NILC algorithm to partial-sky B-mode data; in Sect. 4 you can find the procedures adopted to assess the performance of ILC methods; in Sect. 5 the obtained results from the application of NILC to LSPE-SWIPE and SO simulated data; finally, we report our conclusions in Sect.6.
the HEALPix pixelisation scheme (Górski et al. 2005) with N side = 128, which corresponds to a pixel resolution of 27.5 ′ .The Q and U maps at the different frequencies are obtained from the co-addition of three separate components: CMB, Galactic emission (synchrotron and thermal dust) and Gaussian white and isotropic noise.We do not include polarised extra-Galactic sources, since they are expected to be dominant on angular scales smaller than those of main interest for inflationary B-mode science (ℓ > 100) (Puglisi et al. 2018).
The CMB maps are generated with the HEALPix Python package1 as random Gaussian realisations from the angular power spectra of the best-fit Planck-2018 parameters (Planck Collaboration et al. 2020b) with tensor-to-scalar ratio r = 0, unless otherwise specified.Instrumental noise is simulated as white and isotropic Gaussian realisations with standard deviations in each pixel associated to the polarisation sensitivities listed in Table 1.Galactic emission is generated using the PySM Python package2 (Thorne et al. 2017).Synchrotron polarised emission is modelled with a simple power law (Rybicki & Lightman 1979): with X = {Q, U}, n the position in the sky and ν the frequency.X s (n, ν 0 ) represents the synchrotron template at a reference frequency ν 0 .Thermal dust emission is modelled with a modified black-body (MBB): where B ν (T ) is the black-body spectrum.
We consider two different Galactic models, commonly adopted within the forecast analyses of CMB polarisation experiments: -d0s0, where β s = −3, β d = 1.54 and T d = 20 K are assumed to be constant across the sky; -d1s1, where the spectral indices of synchrotron and thermal dust emission depend on the position in the sky.
In both cases, the dust template is the estimated dust emission at 353 GHz in polarisation from the Planck-2015 analysis (Planck Collaboration et al. 2016a), smoothed with a Gaussian kernel of FWHM=2.6 • and with small scales added by the procedure described in Thorne et al. (2017).The synchrotron template is the WMAP 9-year 23-GHz Q/U map (Bennett et al. 2013), smoothed with a Gaussian kernel of FWHM=5 • and small scales added.
In the d1s1 model, dust temperature and spectral index maps are obtained from the Planck data using the Commander pipeline (Planck Collaboration et al. 2016a).The synchrotron spectral index map is derived using a combination of Haslam 408 -MHz observations and WMAP 23-GHz 7-year data (Miville-Deschênes et al. 2008).The d0s0 model is a simplified representation of Galactic polarised emission, because we know that the spectral properties of the foregrounds vary in the sky.However, it still represents an important starting point for assessing the performance of blind methods on polarisation B-mode data-sets.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 data-set with the lowest angular resolution and then added together.The sky patches observed by LSPE-SWIPE and SO SAT are shown in Fig. 1.The SO SAT footprint has a sky coverage of f sky = 34%.However, due to the highly inhomogeneous scanning of the telescope, some regions of the sky will be much better observed than others, leading to an effective sky fraction of ∼ 10% (see Ade et al. 2019).Therefore, to perform a realistic foregrounds subtraction, the ILC methods are applied to the simulated SO SAT dataset within the reduced patch shown in Fig. 2 and obtained considering only the 10% of the pixels with the highest values in the hit counts map provided by the SO Collaboration.Consistently, the sensitivity values reported in Table 1 refer to a homogeneous hit counts map with a sky fraction of f sky = 10% (Ade et al. 2019).

Extension of ILC methods to partial-sky observations
The Internal Linear Combination (ILC) method 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 Wilkinson Microwave Anisotropy Probe (WMAP) satellite (Bennett et al. 2003).Then, throughout the years, several extensions of the basic algorithm have been proposed.The methods applicable to polarisation data considered in this work are: -ILC, -Polarisation ILC (PILC) (Fernández-Cobos et al. 2016), -Needlet ILC (NILC) (Delabrouille et al. 2009).
All the above ILC methods consist in linearly combining the N ν multi-frequency maps of one or multiple CMB experiments with frequency-dependent weights ω: where X is the B-field for ILC, P = Q + iU for PILC and a set of B-mode needlet coefficients for NILC (as discussed in the text below).The weights are derived by minimising the variance of X to reduce the contamination of Galactic emission and instrumental noise.Therefore, among the component separation techniques, ILC methods require the smallest number of a priori assumptions.Recently, a further variation of these techniques has been proposed, the constrained Moments ILC (cMILC), where the weights estimation is constrained to de-project some moments of the foregrounds emission (Remazeilles et al. 2021).In this work, however, we will 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 N ν j=1 ω j = 1.This last condition, together with Eqs. 3 and 4, enables us to demonstrate that the output solution contains the full CMB signal and some foregrounds and noise residuals: (5) In the case of simulated data-sets, the estimation of foregrounds and noise residuals is straightforward, while for actual data they have to be estimated through Monte Carlo simulations based on our prior knowledge of the contamination in the input multi-frequency data-set.
The weights which minimise the output variance and whose sum is equal to unity can be analytically estimated for ILC, PILC and NILC (see Bennett et al. 2003,Fernández-Cobos et al. 2016,Delabrouille et al. 2009).The only assumption of these methods is that the CMB has a known emission law (black-body spectrum) and no correlations with foregrounds or noise.NILC, which is the main algorithm considered in this work, represents a refinement of the ILC method (Delabrouille et al. 2009), since the effectiveness of variance minimisation is enhanced when performed in the needlet domain.
Needlets are a particular wavelet system that guarantees simultaneous localisation in harmonic and pixel space.They have been introduced into the statistical literature by Baldi et al. (2009) and have been applied for the first time to CMB data in Pietrobon et al. (2006).
In practise, the needlet coefficients of an input B-mode map at frequency i, β i j , are obtained by filtering its harmonic coefficients a i ℓm 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 Article number, page 4 of 23 a way as to obtain a minimum variance map β NILC j at each scale j: The pixel-dependent weights are computed as follows: where the covariance matrix C j ik (γ) = ⟨β i j • β k j ⟩ for scale j at pixel γ is estimated as the average of the product of needlet coefficients in some space domain D. 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 a NILC ℓm, j 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 sub-space.This effect has to be taken into account when estimating the CMB power spectrum.The NILC bias problem in this analysis is further discussed in the Appendix A.
In the case of ground-based and balloon-borne experiments, ILC methods cannot be applied straightforwardly to Bmode data.Indeed, three main complications arise: - The E-B leakage effect is associated to the fact that the E-B decomposition of partial-sky Stokes parameters is not exact: some modes (ambiguous modes) satisfy both the Eand B-conditions simultaneously.When we split the polarisation field into an E-and a B-part, these ambiguous modes can go into either component (Bunn et al. 2003).In the case of CMB, where E modes are much brighter than B modes, this leads to an over-estimation of the power in B modes, especially on large scales.The needlet filtering, on the other hand, represents a convolution of the data.If B modes are observed only in a portion of the sky, some modes (on large scales) can be lost in the process if the configuration of the needlet bands is not carefully chosen, because the convolution mixes the signal of the pixels close to the border with the null one of the unobserved region.Both these effects are relevant especially on those angular scales which are most sensitive to the amplitude of primordial tensor perturbations, and thus they must be properly analysed and treated to have a correct reconstruction of the primordial tensor CMB B-mode signal.
Finally, the ILC algorithms require input maps with a common angular resolution (which usually is the one of the frequency channel with the largest beam).The convolution for a beam, in analogy with the needlet filtering, is not trivial if one deals with observations only on a portion of the sky due to the leakage of null signal of the unobserved region within the patch.The relevance and possible correction of such issues are explored in Sects.3.1, 3.2 and 3.3 by comparing the angular power spectra (evaluated in the footprint of the considered experiments) of leakage-corrected, needlet-filtered, and smoothed CMB simulations with those of the input exact B-mode signal, reconstructed with a full-sky Bdecomposition of input Q and U maps.In such analysis, we consider the entire footprint observed by SO, given in the right panel of Fig. 1, since this is the region over which the B modes will be reconstructed from Q and U maps.Specifically, we use the estimator of the pseudo-angular power spectrum, which indicates the rotationally invariant variance of the harmonic coefficients of a map: It is customary to report an analogous quantity: 2π C ℓ , which would be constant on large scales for a scale-invariant primordial density perturbations spectrum.In this work, the estimation of the power spectra of the B-mode maps is performed with the NaMaster Python package3 that extends the estimator in Eq. 10 taking into account the effects of masking and beam convolution (see Alonso et al. 2019).

Correction of the E-B leakage
The polarisation field of the CMB is a spin-2 quantity that can be described in terms of the Stokes parameters Q and U as follows: Work of Liu, Creswell et al. (https://arxiv.org/pdf/1811.04691.pdf) Alessandro Carones (UTOV) NILC for cut-skies May 12, 2022 7 / 28 Fig. 3. Recycling method for E-B leakage correction where γ denotes the position in the sky.This field can be expanded in spin-weighted spherical harmonics: From the harmonic coefficients in Eq. 12, it is possible to build the polarisation fields E (a scalar map of even-parity) and B (a pseudo-scalar field of odd-parity): The construction of maps from E ℓm and B ℓm represents the E-and B-decomposition of a set of Q and U maps.
For partial-sky observations, we do not have information on Q and U on the entire sky and, therefore, Eqs.11, 12 and 13 lead to an incorrect reconstruction of E and B maps.
The E-B decomposition on the cut-sky is indeed not unique due to the non-locality of the transformation and, therefore, some E modes will be interpreted as B modes and vice versa (ambiguous modes) (see Lewis et al. 2001).The CMB 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 partial-sky observations, leading to an over-estimate of the power in the reconstructed CMB B-mode map.Several different methods have been proposed to address the leakage correction problem (Lewis et al. 2001;Bunn et al. 2003;Kim & Naselsky 2010;Liu et al. 2019;Zhao & Baskaran 2010;Kim 2011;Ghosh et al. 2021).In this work, we consider: the recycling method introduced in Liu et al. ( 2019) the ZB method presented in Zhao & Baskaran (2010).
We test these techniques on 200 CMB-only simulations that include lensing from E modes and primordial tensor perturbations with r = 0.01 for the cases of the SWIPE and SO SAT patches (shown in Fig. 1).We have chosen this value of the tensor-to-scalar ratio because it is close to the upper bound targeted at the 95% confidence level (CL) by LSPE in the case of no detection (Addamo et al. 2021) and represents the amplitude expected to be observed at high significance by SO (Ade et al. 2019).
The performance of these methods has been evaluated by comparing the angular power spectra, D BB ℓ and D BB ℓ,in , of leakage-corrected and exact B-mode maps.The exact B modes are reconstructed with a full-sky B-decomposition of CMB Q and U simulations.Furthermore, we consider an effective tensor-to-scalar ratio: , 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.

The recycling method
The recycling method has been introduced in Liu et al. (2019).The procedure is summarised in Fig. 3.It consists in decomposing the masked polarisation field P = (Q, U) into the so-called E-and B-family: ), which receive contributions only from Eand B-mode 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 leakage-corrected B-family The fit is performed with a simple leastsquares linear regression method.We then obtain the final CMB B-mode map through a B-decomposition of P ′′ B .The leakage correction of the recycling method is not exact, due to our lack of knowledge of the Q and U parameters in the unobserved region of the sky, and therefore some residuals are still present in the reconstructed map.The performance of the method on 200 CMB-only simulations with lensing+r=0.01for the SWIPE and SO SAT patches is shown in the upper panels of Fig. 4, where the average angular power spectrum of the corrected B maps is compared to that of the exact signal reconstructed with a full-sky B-decomposition of the simulated Q and U maps.As can be seen in the top left panel, if not corrected, the E-B leakage highly contaminates the estimation of the angular power spectrum over most of the multipole range of interest.The leakage correction is performed on the full SWIPE and SO SAT footprints; the angular power spectra are estimated, instead, masking the 4% of pixels closest to the borders, signal with r = 0.01.Solid lines represent the power spectrum of exact B-maps (input), reconstructed with a full-sky B-decomposition of Q and U and then masked for power spectrum estimation; circles and diamonds that of leakage-corrected maps (output LC) for SWIPE and SO SAT, respectively; the dashed lines the leakage residuals after the correction; the dashed-dotted lines (if present) the CMB B modes without any leakage correction (output noLC).On the right: the effective tensor-to-scalar ratio associated to the absolute difference between leakage-corrected and exact angular power spectra.Top panel shows the results when the recycling method is applied to CMB B-mode maps for the SWIPE (red) and SO SAT (blue) footprints; middle panel those obtained when iterative B-decomposition with no (red), 1 (green), 2 (cyan) or 3 (magenta) iterations is performed; in the bottom panel, we compare the performance of iterative decomposition with three iterations (magenta) and diffusive inpainting (orange) in correcting the residual leakage.The middle and bottom panels present results only for SWIPE.The adopted binning scheme is ∆ℓ = 3 for ℓ ≤ 4, ∆ℓ = 6 for 5 ≤ ℓ ≤ 28 and ∆ℓ = 15 for ℓ ≥ 29 for the spectra computed in the LSPE-SWIPE patch, while a constant ∆ℓ = 15 for those estimated in the SO footprint.The error-bars highlight the uncertainty on the mean at 1σ, estimated from the dispersion of the angular power spectra of the simulations divided by √ N sims .See text for details.
where the residual leakage is expected to still affect the reconstruction of the CMB B-mode power spectrum.It is possible to observe that the residuals after the application of the recycling method are negligible for the SO SAT case in the entire multipole range of interest (30 ≲ ℓ ≲ 130), if compared to the upper bound on r targeted by the experiment (r ≲ 0.003 at 1σ, Ade et al. 2019).Indeed, the effective tensor-to-scalar ratio associated to 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 B-mode signal becomes dominant while the expected primordial tensor spectrum decreases, and even a small relative error in the reconstruction results in a larger value of r leak .Given the capability of LSPE to measure primordial B modes on the largest angular scales, it would be desirable to lower the E-B leakage residuals in that multipole range with a further processing of the maps corrected with the recycling method.To this end, we implement two separate extensions: the iterative B-decomposition the diffusive inpainting.
Both proposed techniques are not applied to the SO SAT data-set, since in this case optimal results in correcting the E-B leakage effect are already obtained with the standard recycling method.
Iterative B-decomposition is summarised in Fig. 5 and consists of iteratively reconstructing a B-family, (n) , by performing a B-decomposition of the masked P (n−1) B = (Q B , U B ) (n−1) .The starting set of Stokes parameters for the iterations, P (1) (1) , is the leakage-corrected B-family 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 de-projected, because interpreted as E in one of the iterations.The results of the iterative B-decomposition are shown in the second row of Fig. 4, where we plot the average angular power spectra of the 200 CMB-only simulations when either the recycling method alone (iter 0) or a post-processing with 1, 2 or 3 iterative B-decompositions is applied.The plots highlight how, with this additional procedure, we are able to de-project the residual E-B leakage contamination in the B-mode CMB maps at most of the angular scales of interest, obtaining, when three iterations are performed, residuals at a much lower level (r leak ∼ 10 −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 full-sky analysis, are erroneously estimated as E modes, causing a loss of power in the final B-mode map.This effect is particularly relevant on large scales and when E-B leakage residuals are very low.This explains our inability to correctly reconstruct the input power at ℓ ≤ 4 and the larger (but still irrelevant) errors in the map at ℓ ≥ 20 when several iterations are performed with respect to the case with just one.In these cases, the B-E leakage can become the limiting factor in our ability to reconstruct the input CMB B-mode signal, possibly leading to an under-estimate of the input tensorto-scalar ratio.This is particularly meaningful for those experiments (such as LSPE) targeting the observation of the reionisation peak.It is important, then, to determine which multipoles are critically affected by this phenomenon and to exclude them from the cosmological analysis.The effective tensor-to-scalar ratio fitted on the obtained residuals (see the right panels of Fig. 4) highlights that the first three multipoles (ℓ ≤ 4) suffer a significant B-E leakage in the SWIPE case (r leak ∼ 7 • 10 −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 under-estimate of its value if compared to the input one.The applicability of iterative B-decomposition may present a drawback if the use of a large number of iterations leads to a considerable power loss due to B-E leakage along with the subtraction of residual ambiguous E-mode contamination.To assess whether this is the case, we perform the iterative B-decomposition on two separate masked initial Bfamilies: P onlyB B and P leak B .The former is the set of Stokes parameters reconstructed with a full-sky B-decomposition of input Q and U and represents the ideal goal of a leakage correction method.P leak B , instead, is composed of Q and U maps sourced only by the residual E-B leakage after the application of the recycling method.In Fig. 6, we compare the angular power spectra of the output B-mode maps from the B-decomposition performed on the two families after several iterations with those of the initial P onlyB B and P leak B .It is possible to observe that for ℓ ≥ 5, iterations do not deproject modes of the B maps obtained from the initial P onlyB B and at the same time are able to subtract most of the leakage residuals (60% of the power already after the first iteration).At ℓ < 5, instead, the iterative B-decomposition suffers loss of power in P onlyB B due to B-E leakage, but these scales would have been excluded anyway from the cosmological analysis of the NILC solution because an analogous significant loss is observed even without any iteration.Therefore, Fig. 6 highlights the benefit of employing the iterative Bdecomposition.Given the previous results, we apply the NILC pipeline on multi-frequency B-mode maps reconstructed in the SWIPE patch and processed with the recycling method and three iterations of B-decomposition.Hereafter, we refer to it as It is based on the fact that ambiguous modes ψ satisfy the spherical bi-Laplacian 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 bi-Laplacian with the Laplacian and neglecting the Neumann boundary conditions.In this case, a template of residual ambiguous modes in the recycling CMB B-map is given by diffusive inpainting.The procedure consists in imposing as boundary conditions the values of the pixels of the recycling CMB solution at the edges of the patch and replacing iteratively inner sky pixels in the footprint with the average of their neighbours.The obtained template of the residuals is then subtracted from the input recycling CMB map.
The results of the application of the diffusive inpainting on the solutions of the recycling method for 200 CMB-only simulations with lensing+r=0.01 on the LSPE-SWIPE patch are shown in the bottom panels of Fig. 4.They are compared with those obtained with iterative Bdecomposition 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 Bdecomposition, just including ℓ = 4 in the second bin of the angular power spectra of the reconstructed CMB maps leads to a significant under-estimate of its value with respect to that of input B modes.The comparison of r leak in Fig. 4 when diffusive inpainting or iterative B-decomposition is applied highlights that the latter performs better on the largest scales (5 ≤ ℓ ≤ 20), while they are equivalent for intermediate and small ones (ℓ > 20).Given these results, the NILC method is also applied to Bmaps reconstructed in the SWIPE patch and corrected with the recycling method and the inpainting technique.Hereafter we refer to it as rin-NILC.Finally, Fig. 7 shows that, independently of the amplitude of the input tensor perturbations, iterative B-decomposition with three iterations and diffusive inpainting allow us to perfectly reconstruct the B-mode CMB anisotropies for ℓ ≥ 5.
We generate CMB simulations with different input amplitudes of tensor perturbations and obtain leakage-corrected B-mode maps with the three techniques described above: the standard recycling method, iterative recycling, and inpainting recycling.Then, a tensor-to-scalar ratio is fitted on the average binned angular power spectrum (C BB ℓ b ) of the output CMB B-mode maps for all different cases and excluding the first three multipoles (ℓ ≤ 4).The fit of the tensor-to-scalar ratio is performed with a Gaussian likelihood (Hamimeche & Lewis 2008;Gerbino et al. 2020): where C r=1 ℓ b is the binned angular power spectrum of primordial tensor CMB B modes for a tensor-to-scalar ratio r = 1, C lens ℓ b the binned BB lensing angular power spectrum and M ℓ b ℓ ′ b the covariance matrix associated to C BB ℓ b .A binning scheme of ∆ ℓ = 15 is adopted to make the angular power spectrum Gaussianly distributed.We have considered an input r that varies between the sensitivity of LSPE and the current upper bounds reported in Tristram et al. (2022).In every case, with both iterative recycling and inpainting recycling methods, we are able to recover the same results obtained in an ideal case without any E-B leakage contamination in the maps.Instead, with the standard recycling method, a bias is observable, especially for lower values of the input tensor-to-scalar ratio, because of the residual E-B leakage contamination on large scales.

The ZB method
The second technique we consider for the correction of E-B leakage is the ZB method introduced in Zhao & Baskaran (2010).In this case, two different but related definitions for scalar and pseudo-scalar polarisation fields are employed: where ðs and ð s are the spin-lowering and raising operators: The harmonic coefficients E ℓm and B ℓm of E and B are related to E ℓm and B ℓm of Eq. 13 by a ℓ-dependent numerical factor N ℓ : In the case of partial-sky observations, the spin operators have to be applied on P + and P − fields multiplied by a window function W(γ) (non-zero only in the observed patch of the sky), and the recovered quantities Ẽ and B will be affected by the E-B leakage: In Zhao & Baskaran (2010), the correction term to reconstruct the pure pseudo-scalar field B in the observed region of the sky has been derived as: where: with sin θ∂ϕ∂θ .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 : , 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 B field with this procedure is exact.However, the ill-behaved 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 B are obtained by computing E ℓm and B ℓm of Eq. 13 with an harmonic decomposition of masked Q and U maps and then filtering them with N ℓ .Fig. 8 compares the mean angular power spectrum of 200 leakage-corrected CMB-only simulations with lensing+r = 0.01 for the SWIPE and SO patches with that of the exact signal reconstructed with a full-sky E-B decomposition of the simulated Q and U maps.In these cases, all angular power spectra have been computed employing masks obtained from the exclusion of pixels with values lower than 0.01 in the original apodised patches and apodising them with a 7 • apodisation scale, leading to final sky fractions of f sky = 32% and 29%, respectively, for SWIPE and SO.In both footprints, the power of the residuals is lower across the entire multipole range of interest than the sensitivities targeted by the experiments, with associated values of the effective tensor-to-scalar ratio: 10 −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 B-mode power spectrum at ℓ < 5.This is motivated by the fact that, in this case, in contrast to the application of the recycling method, the correction term in Eq.21 is exact and we can employ apodised masks which reduce both the E-B and the B-E leakage effects.Thus, these large angular scales are included in the cosmological analysis of the NILC CMB solutions when the ZB method is adopted.

Needlet filtering of partial-sky maps
The second step in the pre-processing of B-mode data which could introduce distortions in the CMB reconstruction before applying a NILC minimisation is the needlet filtering of partial-sky observations.Needlet coefficients, indeed, are computed by convolving a map with a filter in pixel-space.If observations only cover a partial fraction of the sky, the signal in the pixels close to the borders of the observed patch is mixed with the null values of the unobserved region by the convolution procedure.Thus, the estimation of needlet coefficients can be highly affected by this operation.The wider the needlet band is in pixel space, the larger this effect will be.To assess the impact of needlet filtering, we generate 200 CMB B-mode maps from a full-sky E-B decomposition of Q and U, simulated with Planck 2018 best-fit parameters and lensing+r=0.01.Then, we apply a forward and an inverse needlet transformation to these B-mode maps, properly masked to recover the SWIPE and SO footprints.We assess the goodness of the reconstruction by comparing the angular power spectra of maps obtained after a forward and inverse needlet transform, D BB ℓ,nl , and those of the input ones, D BB ℓ,in .In analogy to the analysis of E-B leakage correction, we consider an effective tensor-to-scalar ratio r nl : , which quantifies the absolute error we can make in the partial-sky reconstruction of both the tensor modes and the lensing signal due to the needlet filtering.In this analysis and in the rest of the paper, we adopt mexican needlets (Geller & Mayeli 2008), which are characterised by the following Gaussian-related weighting function: 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.
Top panels of Fig. 10 show that filtering masked B-mode maps with such a needlet configuration does not introduce any significant error in the reconstruction of the CMB signal on a partial sky.Indeed, the residuals between the input maps and those obtained after a forward and an inverse needlet transform are at the level of r nl ≲ 10 −4 for both SWIPE and SO at all angular scales (see the upper right panel of Fig. 10).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: the mean angular power spectra (over 200 CMB-only simulations) of the input exact B-maps reconstructed with a full-sky EB-decomposition of Q and U (solid lines); those of the reconstructed maps after needlet filtering or beam convolution (dots and diamonds); the corresponding residuals (dashed lines).The plot on the right represents the effective tensor-toscalar ratios associated to the absolute differences between output power spectra after needlet filtering and beam convolution and exact input ones in the SWIPE (red) and SO (blue) patches.A binning scheme of ∆ℓ = 3 for ℓ ≤ 4, ∆ℓ = 6 for 5 ≤ ℓ ≤ 28 and ∆ℓ = 15 for ℓ ≥ 29 is employed for the spectra computed in the LSPE-SWIPE patch, while a constant ∆ℓ = 15 for those estimated in the SO footprint.The error-bars highlight the uncertainty on the mean of the simulations at 1σ.

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 data-sets in this analysis: LSPE-SWIPE complemented with Planck frequencies (SWIPE+PLANCK) and SO-SAT.For SWIPE+Planck it suffices to bring the full-sky Planck maps to the lower SWIPE angular resolution of FWHM=85 ′ .For SO SAT, instead, where frequency channels have different angular resolutions, such an operation is not trivial.The convolution for a beam of partial-sky maps can lead to an error in the CMB reconstruction, as in the case of needlet filtering, because the values of the pixels close to the border of the patch may be influenced by the null values outside the observed region.To assess the relevance of this effect, we generate 200 CMB B-mode simulated maps with lensing+r=0.01at the highest angular resolution of SO (FWHM =9 ′ ) and compare the average angular power spectrum, computed in the SO SAT patch, of these maps when they are brought to the lowest resolution of SO (FWHM=91 ′ ) full-sky (D BB ℓ,in ) or after having been masked (D BB ℓ,smooth ).The comparison between these spectra is shown in the bottom panels of Fig. 10 together with an assessment of the residuals through an effective tensor-to-scalar ratio: It is possible to observe that the most aggressive convolution, which has to be applied to SO maps, does not considerably affect the B-mode CMB reconstruction with values of r smooth of the order of or below 10 −4 at all the angular scales of interest for B-mode science.

NILC performance: analysis setup
After having applied a leakage correction method and needlet filtering to input partial-sky B-mode multifrequency maps, the NILC application leads to the best blind estimate of the CMB B-mode signal in the considered patch.Such a solution, as shown in Eq. 5, is contaminated by residuals of the polarised foregrounds emission, especially close to the Galactic plane.For this reason, when computing the angular power spectrum for a cosmological analysis, the expected most contaminated regions are masked, lowering the effective sky fraction.The masking strategies adopted in this work are described in Sect.4.1.
In this work, where we analyse simulated data, we can precisely determine the expected level of foregrounds and noise residual contamination in the NILC CMB solution by combining the input foregrounds and noise simulations at different frequencies with the NILC weights.Therefore, it is possible to easily assess the impact of these residuals on the estimate of the tensor-to-scalar ratio following the approach detailed in Sect.4.2.
In addition to the residual foregrounds contamination, the other main issue of the NILC CMB reconstruction is the loss of modes induced by the empirical correlation of the reconstructed CMB with other components, especially on large scales, which might lead to a negative bias in the estimated power spectrum.However, we have verified that this effect is negligible in our analysis, given the adopted configuration of needlet bands (see Appendix A).

Foregrounds masks
In order to perform a cosmological analysis of a NILC CMB solution, masking the most foregrounds contaminated regions is usually needed before estimating the angular power spectrum.We have adopted two distinct and alternative methods to generate such masks: 1. computing the polarised intensity P 2 = Q 2 + U 2 of simulated foregrounds maps at a CMB frequency channel of the considered data-set, smoothing it with a Gaussian beam of FWHM=12.5 • , thresholding it to obtain masks with decreasing f sky ; 2. considering directly the average foreground residuals of the method under study, thresholding this map to obtain binary masks with decreasing f sky , smoothing them with a Gaussian beam of FWHM=8 • , finally thresholding the smoothed masks to get binary masks with f sky equal to the original ones.
In both cases, the considered SO SAT patch is shown in Fig. 2, and corresponds, as already discussed in Sect.2, to the region where the ILC methods are applied.The rationale of the first masking strategy is to have common masks for all the considered methods (NILC, ILC and PILC).In such a way, we can assess and compare their foregrounds residuals in the same patch of the sky.These masks are obtained by assuming to have perfect knowledge of the polarised Galactic emission across the sky at a CMB channel.This is not possible when dealing with real data.However, similar masks can be obtained from foreground templates built by re-scaling with appropriate SEDs the Q and U maps of the considered data-set at the lowest (for synchrotron) and highest (for thermal dust) frequencies.The smoothing procedure is performed on the map to be thresholded in order to attenuate the small-scale fluctuations, which in real data are mainly induced by noise and CMB, and thus to avoid masks with a patchy structure.Instead, the second masking approach is more 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 ill-behaved 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 SO SAT) that we consider for the computation of the angular power spectra.
In the case of the NILC application on SWIPE+Planck Bmode maps corrected with the recycling method, we ex-clude the 4% of the pixels closest to the patch border, where some residual leakage could affect the estimate of the power spectra.This step, instead, is not needed for SO SAT, because the reduced patch of Fig. 2 already excludes the borders of the original footprint.
In this second masking strategy, the sky fractions for each case have been determined by comparing the level of the foregrounds residuals with the primordial tensor signal targeted by the experiment.

Estimates of the residuals contamination
When NILC is applied on simulated data-sets, we have access to its foregrounds and noise residuals in the B-mode solution.Therefore, it is possible to compare their angular power spectra with a primordial BB power spectrum of tensor modes with r equal to the target of the experiment under consideration.
If cross correlations among components are negligible (as proved in the Appendix A), the power spectrum of the cleaned map (C out ℓ ) can be written as follows: where C f gds ℓ and C noi ℓ are the angular power spectra of residual foregrounds and noise.The noise bias, C noi ℓ , can be subtracted at the power spectrum level by estimating its contribution through Monte Carlo simulations or by computing cross-spectra between uncorrelated splits of the data (e.g.maps from subsets of detectors or from partial-mission observations).Hence, the main systematic bias to the final angular power spectrum will be given by the residual Galactic contamination.
To assess the impact of this residual term on the estimation of r from C out ℓ , we fit r on binned C f gds ℓ b assuming a Gaussian likelihood (Hamimeche & Lewis 2008;Gerbino et al. 2020): where C r=1 ℓ b is the primordial tensor CMB B-mode angular power spectrum for a tensor-to-scalar ratio r = 1, and M ℓ b ℓ ′ b is the covariance matrix for a fiducial cosmological model with r = 0.An aggressive binning makes the angular power spectrum gaussianly distributed and mitigates correlations among different modes (induced by the use of a mask).In this work, we adopt a constant binning of ∆ ℓ = 15.We have tested the robustness of the obtained constraints by both varying the binning scheme and even considering an inverse-Wishart likelihood function.The covariance matrix is estimated from the angular power spectra of 200 NILC simulated solutions as: where A L quantifies our ability to de-lens the output map with A L = 1 for no de-lensing and A L = 0 for full delensing.Eq. 27 accounts for the cosmic variance associated to the residual lensing signal, the sample variance of the residual foreground and noise angular power spectra and all the cross-terms.In this work, we assume A L = 0.5 (which is the level that can be achieved with multi-tracers techniques; see, e.g., Yu et al. 2017).or Simons Observatory (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.

Results
In the processing of future ground-based and balloon-borne CMB polarisation experiments, NILC represents an effective alternative to the commonly used parametric methods, given the very low number of assumptions on which the algorithm is based.In Sect.5.1, we assess the ability of NILC to recover the input CMB B-mode signal in an ideal case where we assume to be able to perfectly reconstruct the needlet and pixel coefficients of B-mode data even for partial-sky observations.We compare these results with those obtained with the application of ILC and PILC.Once we have verified that it is worth applying NILC on B modes, we proceed to report in Sect.5.2 the results of the application of NILC to a more realistic case with input maps leakage-corrected with the recycling or ZB method.In this case, also the needlet filtering and the beam convolution are performed on partial-sky maps.In all the following analyses, the input CMB B-mode signal has no tensor modes and only lensing.We have verified that the same conclusions hold when injecting a primordial inflationary tensor signal.We consider two different data-sets: the three channels of SWIPE, complemented with the seven polarised Planck frequency maps (SWIPE+Planck) the six simulated frequency maps of the Small Aperture Telescope of the Simons Observatory (SO SAT) The elements of the covariance matrix of Eq. 8 associated to the Planck channels are computed considering only the observed SWIPE region (shown in Fig. 1).For SO, the considered patch is shown in Fig. 2 and represents a faithful approximation of the SO SAT sky coverage under the assumption of homogeneous noise.
For the application of NILC, the used configuration of needlet bands in all cases is the same of the test performed in Sect.3.2: mexican needlets with B = 1.5 and the first eleven bands added together.

Do we need NILC for B-modes?
At the moment, PILC (see Sect. 3) is the only blind method available to easily subtract foregrounds emission from ground-based and balloon-borne CMB polarisation data.However, it works with polarisation intensity data, which account for both E and B modes; therefore, the application of ILC and NILC algorithms could represent a valid alternative, since they are sensitive only to Galactic contamination in B modes.
We compare the performance of PILC, ILC and NILC by looking at the angular power spectra of their foregrounds and noise residuals in an ideal case, where we assume to be able to perfectly recover the B-mode signal in needlet and pixel space for partial-sky observations.This is explicitly obtained by reconstructing from simulated full-sky Q and U maps full-sky B-mode maps, which are then smoothed to a common angular resolution and needlet filtered, and only at the end they are masked according to the patch of the considered experiments.In all cases, a model of foregrounds with spatially varying spectral indices is assumed for polarised dust and synchrotron emissions (d1s1).This ideal case with a full-sky reconstruction of B-mode maps is considered just to compare the performance of the different blind methods and to assess the masking strategy needed to reach the scientific targets of the experiments under study.Then, in Sects.5.2.1 and 5.2.2 we report the results for a realistic approach with a reconstruction of Bmode maps from partial-sky Q and U observations.In Fig. 11, the angular power spectra of foregrounds and noise residuals from the application of PILC, ILC and NILC are shown.They have been computed using Galactic masks retaining the 15% (for SWIPE) and 7.5% (for SO) of the sky and generated with the first procedure described in Sect.4.1.We start by considering the SWIPE data-set complemented with Planck frequencies (SWIPE+Planck) as introduced in Sect. 2. A method which minimises the variance directly of B modes (ILC) instead of polarised intensity (PILC) leads to lower residuals of Galactic emission and instrumental noise at all angular scales.The CMB, indeed, is much brighter in P = Q 2 + U 2 (due to the presence of the E-mode signal) than in B and therefore the PILC weights are less capable of tracing and subtracting the frequency-dependent contaminants in the B-mode sky.Furthermore, extending the procedure of variance minimisation to the needlet domain (NILC) yields significantly lower foregrounds residuals, especially on large scales, with respect to ILC and PILC at the price of a mild increase of the noise residuals.This justifies the implementation of the NILC algorithm for cut-sky observations given the sensitiv-  ity to the tensor-to-scalar ratio in that multipole range.The NILC effectiveness is particularly relevant for SWIPE because we have access to information at very low multipoles.For SO, which will observe only modes at ℓ ≥ 30, indeed, NILC foregrounds residuals are comparable with those of ILC and PILC.Indeed, the range of multipoles where NILC performs better in subtracting Galactic contamination is highly dependent on the sensitivity and angular resolution of the considered data-set.However, as expected, the power of total residuals is lower for NILC than for ILC and PILC across all multipoles, because the output variance is separately minimised at different needlet scales.The above results support the effort to extend the application of NILC on B-mode data for ground-based and balloon-borne experiments.

NILC results for realistic SWIPE and SO case-studies
A realistic application of NILC on partial-sky B-mode maps can be summarised with the following steps: -Correct the E-B leakage effect in the multi-frequency data-set, which include CMB, noise and Galactic foregrounds where needed, bring all the frequency maps to the same resolution filter the maps with needlet bands perform the variance minimisation independently on the different needlet scales reconstruct a cleaned CMB map in real space.
Depending on the technique employed to correct the E-B leakage, we can identify several distinct pipelines: • r-NILC, for maps corrected with the recycling method, • rit-NILC, when the recycling method and iterative Bdecomposition are applied, • rin-NILC, for maps corrected with the recycling method and diffusive inpainting, • ZB-NILC, when the ZB technique is performed.
We have verified that all the leakage-correction methods presented in Sect 3.1 perform analogously in the correction of the E-B leakage even in the presence of instrumental noise and Galactic foregrounds.
We will report separately the NILC results obtained when recycling (and its extensions) and ZB methods are applied, for both SWIPE+Planck and SO SAT, in Sect.5.2.1 and 5.2.2, respectively.Figs. 13, 14 and 15 provide a detailed comparison of the two approaches.We analyse simulated data of: -SWIPE+Planck with d0s0 Galactic model -SWIPE+Planck with d1s1 Galactic model -SO SAT with d1s1 Galactic model The case with constant spectral indices for polarised thermal dust and synchrotron emission is considered to have a direct comparison with the recent forecast of the sensitivity on r for the LSPE experiment published in Addamo et al. (2021).Considering the results obtained in Sect.3.1 regarding the effectiveness of the different methods in correcting E-B leakage in the two footprints, we apply rit-NILC, rin-NILC, and ZB-NILC to the SWIPE+Planck data-set, while r-NILC and ZB-NILC to SO SAT simulated maps.In the case of the application of rit-NILC and rin-NILC to SWIPE+Planck, the minimum considered multipole is ℓ = 5 both in the plots and for the cosmological analysis, given the significant loss of power on larger angular scales due to B-E leakage, as found in Sect.3.1.1.In the SO SAT analysis, the leakage correction methods are performed within the full SO footprint (shown in the right panel of Fig. 1), since CMB data will be available in that region.The NILC pipelines, instead, are applied considering the patch of Fig. 2, which retains 10% of the sky and takes into account the highly inhomogeneous scanning strategy of the telescope.The reference sky fraction for the cosmological analysis in each case is established by looking at the results in Fig. 12.In this analysis, we perform a fit of r on the power spectra of NILC foregrounds residuals (see Sect. 4.2) in the ideal case considered in Sect.5.1 adopting the second masking approach described in Sect.4.1 and varying the sky fraction.It is possible to observe that when NILC is applied on d0s0 simulations of the SWIPE+Planck data-set, no masking of the foregrounds residuals is needed to achieve the scientific goal of the mission in terms of sensitivity on the tensor-toscalar ratio r.When a more realistic Galactic model (d1s1) (solid) and noise (dashed) residuals over 200 simulations.On the right: the posterior distributions of an effective tensor-to-scalar ratio fitted on the foregrounds residuals in the case of half de-lensing (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 (rit-NILC and rin-NILC).On the bottom: results when the ZB-(blue) and hybrid ZB-NILC (red) methods are applied.The adopted binning scheme is ∆ℓ = 6 for 5 ≤ ℓ ≤ 28 and ∆ℓ = 15 for ℓ ≥ 29 on the top left, while ∆ℓ = 6 for 2 ≤ ℓ ≤ 19 and ∆ℓ = 15 for ℓ ≥ 20 on the bottom left.
In both cases the chosen binning strategy is just for visualisation purposes.Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.015, 0.03].The lower bound represents the targeted upper limit at 95% CL by LSPE in case of no detection, the upper one the targeted detection at 99.7% CL.The spectra have been estimated without any masking of the regions most contaminated by foregrounds.The final 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.
is instead considered for the same data-set, 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 B-mode data of SO SAT, with a sky fraction of 8%, the foregrounds contamination is sensibly lower than the primordial tensor signal targeted by the experiment, even assuming a Galactic foreground model with varying spectral indices across the sky.
In the following, for all cases, we show the posterior distributions of an effective tensor-to-scalar ratio fitted on the foregrounds residuals average power spectrum, as described in Sect.4.2.When r f gds = 0 is within 2σ, we report upper bounds on r f gds 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 andAde et al. 2019)..

Recycling NILC (r-NILC)
In this section we report the results obtained when the NILC algorithm is applied to multi-frequency B-mode maps that include CMB, noise and foregrounds, and have been corrected with the recycling method (r-NILC).We note that, in the case of the LSPE+Planck data-set, the maps have also been post-processed either with the iterative Bdecomposition with three iterations (rit-NILC) or with the diffusive inpainting (rin-NILC).Recently, updated forecasts on r for the LSPE experiment have been published in Addamo et al. (2021) considering a Galactic model with constant spectral indices for polarised thermal dust and synchrotron.Following the same approach, we first test the recycling NILC on the SWIPE+Planck data-set 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.The plot on the upper right-hand side of Fig. 13 shows that the effective tensor-to-scalar ratio associated to foregrounds residuals has a posterior distribution with an upper limit at 95% CL of ∼ 1.5 • 10 −2 , which is fully in accordance with   (solid) and noise (dashed) residuals over 200 simulations.On the right: the posterior distributions of an effective tensor-to-scalar ratio fitted to the foregrounds residuals in the case of half de-lensing (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 4.2).On the top: blue and red lines represent the results when E-B leakage is corrected with, respectively, the iterative and inpainting recycling methods (rit-NILC and rin-NILC).
On the bottom: results when the ZB-(blue) and hybrid ZB-NILC (red) methods are applied.The adopted binning scheme is ∆ℓ = 6 for 5 ≤ ℓ ≤ 28 and ∆ℓ = 15 for ℓ ≥ 29 on the top left, while ∆ℓ = 6 for 2 ≤ ℓ ≤ 19 and ∆ℓ = 15 for ℓ ≥ 20 on the bottom left.In both cases the chosen binning strategy is just for visualisation purposes.Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.015, 0.03].The lower bound represents the targeted upper limit at 95% CL by LSPE in case of no detection, the upper one the targeted detection at 99.7% CL.The spectra have been estimated masking the regions most contaminated by foregrounds for a final f sky of 12%.
the LSPE target.This constraint has been obtained by excluding the multipoles ℓ ≤ 4 where a significant negative bias is observed in the reconstructed CMB B-mode power spectrum due to the B-E leakage effect (see Sect. 3.1.1).When a more complicated sky model is considered (d1s1 with anisotropic spectral indices for dust and synchrotron), the NILC weights try to trace locally several different SEDs at the same time, leading to a more noisy and less effective removal of the Galactic emission in each pixel.Therefore, in this case, a more aggressive masking strategy (based on the second approach described in Sect.4.1) is needed to exclude the most contaminated regions.We have found in the analysis of Fig. 12 that a sky fraction of f sky = 12% is enough to avoid biases in the estimate of the tensor-toscalar ratio.
The angular power spectra of the residuals are shown in the upper left panel of Fig. 14.Both the rit-NILC and rin-NILC cleaning methods still allow one to possibly detect both peaks in the B-mode primordial spectrum.
The posterior distribution of an effective tensor-to-scalar ratio fitted on the foregrounds residuals has an upper limit of around r ≲ 2.5 • 10 −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.
The very same pipeline can be applied to SO SAT simulations.In this case, the standard recycling method is adopted, given its effectiveness in removing E-B leakage contamination at the angular scales of interest for SO (see Sect. 3.1.1).
The target of SO is the observation of only the recombination bump, as the experiment is not sensitive to multipoles ℓ < 30 due to atmospheric contamination.Therefore, the needlet bands and the cleaning algorithm do not take into account modes on angular scales larger than almost 6 Fig. 15.SO SAT data-set with the d1s1 Galactic model.On the left: the mean angular power spectra of NILC foregrounds (solid) and noise (dashed) residuals over 200 simulations.On the right: the posterior distributions of an effective tensor-to-scalar ratio fitted to the foregrounds residuals in the case of half de-lensing (A L = 0.5).Both for the plots of the angular power spectra and the of the posterior a binning scheme of ∆ℓ = 15 has used.On the top: results when E-B leakage is corrected with the recycling method (r-NILC).On the bottom: results when the ZB-(blue) and hybrid ZB-NILC (red) methods are applied.Everywhere the grey region represents the primordial tensor BB angular power spectrum for r ∈ [0.003, 0.01].The lower bound represents the targeted upper limit at 68% CL by SO in case of no detection, the upper one the targeted detection at few σ significance.The spectra have been estimated masking the regions most contaminated by foregrounds for a final f sky of 8%.
sidering the cleanest 8% of the sky, as suggested by the analysis reported in Fig. 12.

ZB-NILC
ZB-NILC consists of implementing the variance minimisation in needlet domain after having applied the ZB leakage correction method on multi-frequency B-mode maps that include CMB, noise and foregrounds.Therefore, the cleaning algorithm is applied to B maps of Eq. 16 instead of the usual B field.The standard B-mode signal is recovered only at the end when the angular power spectra are computed.We consider the same data-sets and cases of the r-NILC analysis to closely compare the results of the two methodologies by looking at Figs. 13, 14 and 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 data-set with the d0s0 Galactic model are reported.It is possible to observe that the ZB-NILC foregrounds residuals are much higher than those of the recycling NILC.This is caused by the fact that, in ZB-NILC, the input maps are convolved for the harmonic filter N ℓ = √ (ℓ + 2)!/(ℓ − 2)! (see Sect. 3.1 for details), which boosts the power on smaller scales.Therefore, even the weights of the first needlet band are sensitive much more to contamination on intermediate scales (where noise has a not negligible contribution) than to that on the largest ones (where diffused foregrounds dominate).As a consequence, ZB-NILC weights are much less capable of removing Galactic polarised emission in B 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 ZB-NILC (HZB-NILC).HZB-NILC consists of correcting the E-B leakage effect in frequency maps with the ZB method, but then combining them using the NILC weights estimated from the B-mode maps processed with the standard recycling method.The foregrounds and noise residuals of the HZB-NILC are shown with solid and dashed red lines in the bottom panels of Figs. 13, 14 and 15.With such an implementation, we can recover foregrounds and noise residuals at a level comparable to recycling NILC, but at the same time we can exploit the greater leakage correction capabilities of the ZB technique on the largest angular scales (ℓ ≤ 4).In this case, for SWIPE+Planck with the d0s0 model, we obtain an upper bound on an effective tensor-to-scalar ratio fitted on the foregrounds residuals of 5 • 10 −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 r-NILC is motivated by the inclusion in this analysis of the lowest multipoles ℓ ≤ 4. The upper limit on the tensor-to-scalar ratio obtained with HZB-NILC is sensibly lower than that obtained in Addamo et al. (2021) even though the same Galactic model is assumed.Such a difference can be explained taking into account three key factors: i) in this work we do not take into account the scanning strategy and, thus, the anisotropic dis-tribution of the instrumental noise, which may play a role in the performance of the component separation; ii) we assume to be able to de-lens our B-mode power spectrum at a 50% level (Yu et al. 2017), thus reducing the impact of lensing B modes to the total uncertainty in the r reconstruction; iii) for this simplistic Galactic model the minimisation of the variance across the different needlet scales can in principle outperform the parametric methods, which fit spectral parameters in different regions of the sky.The ZB and hybrid ZB-NILC methods are then applied to the LSPE+Planck data-set assuming the more realistic d1s1 model of the Galaxy.The corresponding results are shown in the lower panels of Fig. 14.We again observe different capabilities of the two methods in subtracting Galactic contamination (in favour of the HZB-NILC).However, in this case, the differences between the power spectra are less evident because the considered sky components are intrinsically more difficult to subtract.With a more aggressive Galactic mask with respect to the d0s0 analysis, which retains the 12% of the sky, both the reionisation and recombination peaks with r = 0.015 would be observed.The upper bound on an effective tensor-to-scalar ratio fitted on the foregrounds residuals is r < 0.016 (95% CL) for the Hybrid ZB technique, which is in accordance with the LSPE target in case of no detection.As above, the lower constraint with respect to r-NILC pipelines is given by the possibility of including the first three multipoles in the likelihood.The same trends are obtained when ZB-NILC and Hybrid ZB-NILC have been applied to the simulated SO SAT data-set with Galactic emission simulated with the d1s1 model (see the bottom panel of Fig. 15).The angular power spectrum of the HZB-NILC Galactic residuals is below the curve of the primordial tensor signal targeted by SO (r = 0.01) for ℓ > 30, if we observe the cleanest 8% of the sky.The HZB-NILC foregrounds residuals lead to an upper bound in the posterior distribution of the effective tensorto-scalar ratio of r < 0.0028 (68% CL) which is compatible with the SO target (r ≤ 3 • 10 −3 at 1σ).We observe that for SO SAT the difference in the amplitudes of foregrounds residuals on large scales between ZB-and HZB-NILC is much less evident than the one for SWIPE+Planck.This is due to the fact that, without modes with ℓ ≤ 30, the calculation of HZB-NILC weights in the first needlet band is more affected by noise contamination on intermediate scales.This condition is close to that experienced by the ZB-NILC pipeline.

Conclusions
The main goal of future CMB experiments will be the detection of polarisation B modes generated by primordial tensor perturbations, to definitely confirm the inflationary scenario.Given its low number of assumptions, NILC represents an effective alternative for the analysis of future B-mode data to the commonly used parametric component separation methods, which are more subject to systematic biases in the case of mis-modelling of the foreground properties.However, most future surveys will be balloon-borne or groundbased and will observe only a portion of the sky, whereas, at present, NILC has been successfully applied to full-sky data from either Planck or WMAP.In this work, we explore the possibility to extend the NILC formalism to future B-mode partial-sky observations, specifically addressing the complications that such an application yields: 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 Carlo simulations for two complementary experiments: LSPE-SWIPE and the Small Telescope of Simons Observatory (SO-SAT).The former aims at a detection with high significance of both reionisation and recombination peaks with r = 0.03 or to set an upper bound of r = 0.015 at 95% confidence level; the latter, instead, is designed to be able to measure a signal at the r = 0.01 level at a few σ significance, or exclude it at similar significance, using the amplitude of the B modes around the recombination bump.In the case of LSPE-SWIPE, realistic simulated Planck maps have been included in the component separation pipeline to have a broader frequency coverage and better trace the spectral properties of B-mode foregrounds.When dealing with partial-sky observations, the estimation of B modes from CMB polarisation measurements is challenging due to the mixing of E and B modes.In order to correct for this effect, we implement two different techniques to be applied in pixel space: the recycling (Liu et al. 2019) and the ZB (Zhao & Baskaran 2010) method.We have tested their ability to recover the CMB B-mode angular power spectrum on Monte Carlo simulations, also quantifying a possible bias on the estimate of the tensor-to-scalar ratio.We find that: • the recycling method is able to reduce the E-B leakage contamination in the B-mode maps at a negligible level for SO at all the angular scales of interest (ℓ ≥ 30), while for LSPE-SWIPE only at ℓ ≳ 15 (see Fig. 4).• the ZB method, instead, allows us to recover the input CMB signal at all angular scales for both the considered experiments (see Fig. 8).However, this technique yields some complications: the need to work with the B field of Eq. 16 and to apodise the mask.
• the E-B leakage residuals of the recycling method on large angular scales (ℓ ≲ 15) in the SWIPE patch can be further reduced for ℓ ≥ 5 by applying either of the following two post-processing prescriptions: 1. the diffusive inpainting (Liu et al. 2019) 2. the iterative B-decomposition, which we introduce for the first time in this paper.Ambiguous modes, still present in the map, are 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 3. Upper bounds and constraints of an effective tensor-to-scalar ratio fitted on the average angular power spectrum of foregrounds residuals given by the application of NILC for the listed cases.
and Galactic foreground emission.For LSPE-SWIPE we have considered two different foreground models: one in which the spectral indices of the synchrotron and dust emissions are assumed to be position independent (d0s0), and another in which they vary across the sky (d1s1).For SO-SAT, instead, we generate 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 to be applied on partial-sky B-mode data from future experiments.A detailed forecast of the performance of the considered CMB experiments should take into account some additional realistic effects, 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: recycling + iterative decomposition (rit-NILC) for SWIPE+Planck recycling + diffusive inpainting (rin-NILC) for SWIPE+Planck recycling (r-NILC) for SO-SAT -ZB (ZB-NILC) for both SO-SAT and SWIPE+Planck.
For all cases, we quantify the ability of the pipeline to mitigate foreground contamination in terms of an effective tensor-to-scalar ratio, which is fitted on the angular power spectrum of foregrounds 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, respectively, at 95% and 68% CL for LSPE-SWIPE and SO, as in the official forecasts of the two experiments (see Addamo et al. 2021 andAde et al. 2019).We find that SO-SAT bounds are within the goal of the experiment regardless of the adopted E-B leakage correction approach.In this work, as mentioned above, we did not incorporate a 1/f noise component in the SO SAT simulations, which is necessary for an accurate description of the experiment.However, this component is expected to only mildly affect the performance of the method and the obtained constraints, as also proved in Ade et al. (2019).For SWIPE+Planck with the d1s1 foregrounds model, when the recycling method is applied to correct the E-B leakage, the obtained upper bound on r is greater than the target of the experiment.This result is caused by the need to exclude the first three multipoles where the B-E leakage has a significant impact on the CMB power spectrum reconstruction.For all other cases, instead, the constraints are in agreement with the expected sensitivity also for LSPE-SWIPE.When the ZB-NILC pipeline is applied (see Sect. 5.2.2), the amplitude of the foregrounds residuals is much greater than that obtained with r-NILC.This can be explained by the fact that in this case the input maps are convolved with the harmonic filter N ℓ = √ (ℓ + 2)!/(ℓ − 2)! (see Sect. 3.1.2for 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 Hybrid ZB-NILC (HZB-NILC), where the E-B leakage effect is corrected with the ZB method, while the NILC weights are estimated on the B-mode maps corrected with the recycling technique.We note that the upper limit on the tensor-to-scalar ratio obtained with the HZB-NILC pipeline and assuming the Galactic model d0s0 is sensibly lower than the one obtained in Addamo et al. (2021).This difference can be motivated by taking into account several factors: i) in this work, we do not consider the scanning strategy and, thus, the anisotropic distribution of the instrumental noise, which may play a role in the performance of the component separation; ii) we assume to be able to de-lens our B-mode power spectrum at a 50% level (Yu et al. 2017), thus reducing the impact of lensing B modes to the total uncertainty in the r reconstruction; iii) for the simplistic d0s0 Galactic model, the minimisation of the variance in needlet domain can in principle outperform parametric methods, which fit spectral parameters in different regions of the sky.In this paper, we presented the development and validation of a pipeline, based on the NILC component separation method, for the analysis of future CMB B-mode data collected by experiments that will observe the microwave sky from the ground and balloons.Taking into account realworld issues due to the incomplete sky coverage and stateof-the-art simulations of Galactic foregrounds, our results demonstrate the effectiveness of NILC, which emerges as a valid alternative to parametric component separation techniques for these kinds of experiments.Furthermore, this study permits to easily extend the application of other ILC based techniques, e.g.cMILC (Remazeilles et al. 2021) and MC-NILC (Carones et al. 2022), to multi-frequency B-mode data from ground-based and balloon-borne experiments.

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

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

Fig. 4 .
Fig. 4. On the left: mean angular power spectrum over N sims = 200 CMB-only simulations that include lensing and a primordial tensor

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

Fig. 11 .
Fig. 11.Mean angular power spectra (over 200 simulations) of foregrounds (solid lines) and noise (dashed lines) residuals when PILC (green), ILC (blue) or NILC (red) are applied on B-mode data (reconstructed on the full-sky and then masked) of SWIPE+Planck (left) or Simons Observatory (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.

Fig. 12 .
Fig. 12. Trend of the tensor-to-scalar ratio r f gds fitted on the NILC foregrounds residuals angular power spectra with the likelihood of Eq. 26 when the sky fraction of the employed mask varies.The adopted masking strategy is the second one described in Sect.4.1: we mask the most contaminated regions by foregrounds residuals.In this case, the input multi-frequency B-mode maps are constructed full-sky and then masked.The considered cases are SWIPE+Planck (left) with both d0s0 (red) and d1s1 (green) Galactic models and SO SAT data-set (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.

Fig. 13 .
Fig. 13.SWIPE+Planck data-set with the d0s0 Galactic model.On the left: the mean angular power spectra of NILC foregrounds

Fig. 14 .
Fig. 14.SWIPE+Planck data-set with the d1s1 Galactic model.On the left: the mean angular power spectra of NILC foregrounds

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 SO-SAT).We report the results for the different E-B leakage correction methods, for needlet filtering and beam convolution.The latter is relevant only for SO-SAT for which maps of the different frequency channels need to be brought at the same resolution.HZB-NILC r < 0.005 (95% CL) r < 0.016 (95% CL) r < 0.0028 (68% CL)