Optimization of foreground moment deprojection for semi-blind CMB polarization reconstruction

Upcoming Cosmic Microwave Background (CMB) experiments, aimed at measuring primordial CMB polarization B-modes, require exquisite control of instrumental systematics and Galactic foreground contamination. Blind minimum-variance techniques, like the Needlet Internal Linear Combination (NILC), have proven effective in reconstructing the CMB polarization signal and mitigating foregrounds and systematics across diverse sky models without suffering from foreground mismodelling errors. Still, residual foreground contamination from NILC may bias the recovered CMB polarization at large angular scales when confronted with the most complex foreground scenarios. By adding constraints to NILC to deproject statistical moments of the Galactic emission, the Constrained Moment ILC (cMILC) method has been demonstrated to further enhance foreground subtraction, albeit with an associated increase in overall noise variance. Faced with this trade-off between foreground bias reduction and overall variance minimization, there is still no recipe on which moments to deproject and which are better suited for blind variance minimization. To address this, we introduce the optimized cMILC (ocMILC) pipeline, which performs full automated optimization of the required number and set of foreground moments to deproject, pivot parameter values, and deprojection coefficients across the sky and angular scales, depending on the actual sky complexity, available frequency coverage, and experiment sensitivity. The optimal number of moments for deprojection, before paying significant noise penalty, is determined through a data diagnosis inspired by the Generalized NILC (GNILC) method. Validated on B-mode simulations of the PICO space mission concept with four challenging foreground models, ocMILC exhibits lower Galactic foreground contamination compared to NILC and cMILC at all angular scales, with limited noise penalty. This multi-layer optimization enables the ocMILC pipeline to achieve unbiased posteriors of the tensor-to-scalar ratio, regardless of foreground complexity.


Introduction
Accurate measurements of the temperature anisotropies of the Cosmic Microwave Background (CMB) have yielded stringent constraints on viable cosmological models [1,2].Future CMB experiments, featuring enhanced sensitivity, will provide additional insights into CMB polarization anisotropies [3][4][5][6].On the one hand, an accurate reconstruction of the CMB E-mode polarization anisotropies will lead to a cosmic-variance limited measurement of the optical depth to last scattering surface τ [5,6], providing new constraints on the cosmic reionization history [7,8] and the sum of neutrino masses [3,9,10].On the other hand, a detection of the primordial CMB B-mode polarization may open a new window on the physics of the early Universe by setting constraints on the tensor-to-scalar ratio r, thereby shedding light on the cosmic inflation scenario [11].
To achieve these goals, future CMB experiments require exquisite control of the contamination from Galactic foregrounds and instrumental systematics.In recent years, several component separation pipelines, both parametric [12][13][14][15][16][17] and blind [18][19][20][21][22][23][24], have been developed or customized for this purpose.Given the huge amplitude discrepancy between Galactic foregrounds and CMB B-modes, these pipelines need to demonstrate robustness against the complex and unknown properties of Galactic foreground emissions, including spatial variations in their spectral energy distribution (SED) along and across the lines of sight [25], as well as resulting foreground spectral distortions [26] and spectral decorrelation [27].
Parametric component separation methods are engineered to yield minimal uncertainty in the retrieval of CMB B-modes, as long as the assumed parametric foreground model aligns with the actual foregrounds in the data.However, given the aforementioned intricacies in the foreground properties, parametric methods are also prone to potential biases in the recovery of the CMB B-mode signal when confronted with incorrect or incomplete modelling of poorly understood foregrounds [28].In contrast, blind minimum-variance approaches to component separation, like the Needlet Internal Linear Combination (NILC) method [21], offer the possibility of extracting the CMB polarization E-and B-mode fields from multifrequency data without making any assumption about the foregrounds, thus avoiding mismodelling errors.However, if the polarized Galactic foregrounds turn out to be highly complex, as in the case of multi-layer contributions to thermal dust emission along the line-of-sight and spectral decorrelation [29], minimum-variance methods like NILC may return CMB polarization maps that are still significantly contaminated by Galactic foreground residuals, as highlighted in [23,24,30].
To overcome these issues, intermediate approaches at the intersection of parametric and blind methods can be pursued, using either partial information on Galactic foregrounds or some effective modelling.For instance, information on the spatial variations of the foregrounds across the sky can be inferred from the data and exploited to optimize the domains over which the NILC variance is minimized, as implemented in the Multi-Clustering NILC method [24].As an alternative semiblind approach, the Constrained Moment ILC (cMILC) method [23] relies on effective foreground modelling through statistical moment expansion around fixed pivot spectral parameters to project out some of these moments through nulling constraints added to NILC, while leaving higher-order moments and unmodelled foregrounds for blind variance minimization.Through nulling constraints for moment deprojection, the cMILC method facilitates the handling of foreground moments and SED distortions resulting from line-of-sight and beam averaging of spatial variations in foreground spectral properties by CMB instruments.In addition, it allows for imperfect parameterization of the foregrounds by leaving any higher-order correction to the moment expansion or any omitted foreground to blind variance minimization.This semi-blind technique has proven more effective than NILC in reducing Galactic contamination in the recovered CMB B-mode map when applied to complex foreground simulations of future CMB satellite experiments such as LiteBIRD and PICO [23,31].
However, the reduction of the residual foreground contamination bias in the cMILC CMB solution is counterbalanced by an increase of the overall variance associated with instrumental noise and unconstrained foreground moments, since each additional nulling constraint in cMILC effectively reduces the number of frequency channels available for variance minimization.Faced with this trade-off between bias reduction and variance increase, certain foreground moments may be more conducive to variance minimization than deprojection, particularly when they exert a dominant influence on the overall variance in the data, while other subdominant moments may be better addressed through deprojection.However, no recipe for determining which moments to deproject and which to retain for variance minimization has been provided thus far.In the baseline cMILC approach, the deprojection of zeroth-order moments is prioritized, anticipating their significant contribution to the bulk of Galactic foreground emission.Nonetheless, the optimal solution might favour the deprojection of higher-order components, relegating dominant moments to blind variance minimization.Moreover, if the actual SED of polarized thermal dust and synchrotron emissions significantly deviates from the prevalent parameterization of a single modified blackbody (MBB) and a power-law, respectively, then the deprojection of SED moments with fixed values of pivot spectral parameters across the sky and angular scales may prove less effective.
In this work, we introduce the optimized cMILC (ocMILC) pipeline, where the choice of moments to deproject, the values of pivot parameters, and the amount of deprojection of each moment is optimized through a fully data-driven approach.Specifically, the optimal number of moments to deproject before paying significant noise penalty is determined depending on the experimental configuration by applying the strategy of the Generalized Needlet ILC (GNILC) method [32].This method provides a blind diagnosis of the actual complexity of the Galactic foreground emission across the sky and angular scales, yielding maps of the local number of independent foreground degrees of freedom as byproducts.Such a diagnostic tool is valuable for informing the foreground model parameterization and tuning component separation pipelines according to the observed sky complexity and experimental configuration.The new ocMILC framework is tested and validated on a set of complex foreground scenarios using comprehensive simulations of the Probe of Inflation and Cosmic Origins (PICO) mission concept [6].We choose this experimental configuration as a proof of concept, given its high sensitivity and broad frequency coverage, allowing us to thoroughly assess the robustness of the proposed approach against various thermal dust and synchrotron models.We note that the validation of the ocMILC pipeline is performed under the assumption of perfectly reconstructed frequency maps from time-ordered data.However, as shown by the latest analyses of the Planck data (NPIPE [33], SRoll [34,35], Cosmoglobe [36]), a parametric model of the Galactic emission may still be needed for the low-level data processing steps, as calibration, bandpass corrections, and multi-component destriping, involved in producing high-fidelity frequency maps.
The paper is organized as follows: In Section 2, we provide a description of the PICO sky simulations, considering four foreground scenarios with different dust models.Section 3 outlines our component separation methodology, first reviewing the basics of the cMILC approach (Section 3.1), then describing the proposed ocMILC pipeline (Section 3.2), which involves a data-driven diagnosis of foreground complexity inspired by GNILC (Section 3.2.1)and subsequent optimization of foreground moment deprojection (Section 3.2.2).Section 4 presents our results for the different foreground scenarios.Our conclusions are summarized in Section 5.

Sky simulations
As a proof of concept, we perform simulations of multifrequency sky maps based on the experimental setup of the PICO mission.PICO, a NASA probe-scale space mission concept, is anticipated to be the most sensitive among the proposed next-generation CMB space missions, targeting a sensitivity on the tensor-to-scalar ratio of r ≤ 0.0002 at a 95% confidence level (CL) and aiming for a cosmicvariance limited measurement of the optical depth to last scattering, σ(τ) = 0.002 [6].To meet these objectives, PICO will conduct full-sky observations in 21 frequency channels spanning from 21 to 799 GHz with high signal-to-noise to disentangle Galactic foreground emission from the CMB polarization signal.Given its extensive frequency coverage and high sensitivity, PICO serves as an ideal experimental framework to assess the robustness and flexibility of the proposed optimization of moment deprojection across various thermal dust and synchrotron models.
The PICO sky maps are generated by co-adding simulated maps of CMB, Galactic foregrounds, and instrumental noise across the 21 frequency bands, according to the instrument specifications reported in [6].
CMB polarization anisotropies are generated using the HEALPix Python package1 [37,38], based on the CMB angular power spectra of the Planck 2018 best-fit parameters [39] with a tensorto-scalar ratio r = 0, and include the contribution from the gravitational lensing effect.
The instrumental noise is simulated as white and isotropic Gaussian realizations, with standard deviations in each pixel corresponding to the polarization sensitivities of the baseline PICO configuration as reported in [6].The combined noise level for this configuration is equal to 0.87 µK • arcmin.In [6,30], it has been demonstrated that the targeted sensitivity r < 0.0002 (95% CL) can be achieved for most foreground scenarios with the so-called Current Best Estimate (CBE) PICO configuration, featuring a combined noise level of 0.61 µK •arcmin.However, for the purpose of our study relying on the more conservative baseline PICO configuration, with a combined noise level of 0.87 µK • arcmin, we estimate a target sensitivity of about r < (0.87) 2 /(0.61) 2 × 0.0002 = 0.0004 (95% CL).Here, the noise levels are squared, as the tensor-to-scalar ratio parameter is associated with power spectrum rather than map noise.In Section 4.2, we show that, by adopting the CBE PICO configuration, we can indeed achieve r < 0.0002 (95% CL) with the ocMILC method.
Maps of polarized Galactic emission are generated using the PySM Python package [40,41].Given the current lack of knowledge about the polarized Galactic emission at the targeted sensitivity, we explore a range of foreground sky models that remain consistent with existing observations.Four distinct models of Galactic thermal dust emission are examined in this study and are described below, using the terminology of the PySM: • d1: The thermal dust emission at each line-of-sight n and each frequency ν is characterized by a modified blackbody (MBB) SED in intensity units where X = {Q, U} are Stokes parameters, B ν (T ) is the blackbody spectrum, β d is the dust spectral index, T d is the dust temperature, and X d (n, ν d ) represents the dust template at a reference frequency ν d .Both β d and T d vary across the sky, as estimated from the Planck data analysis using the Commander algorithm [42].The thermal dust template is given by the Commander dust Q and U maps at 353 GHz of the Planck 2015 data release, smoothed with a Gaussian beam of full width at half maximum (FWHM) 2.6 • , with stochastic small-scale fluctuations added via the procedure described in [40].
• d4: In this model, the thermal dust emission is composed of two distinct MBB components, as fitted to Planck and DIRBE/IRAS data in [43].The temperatures of these two components vary across the sky, while the spectral indices remain constant with values β (1) d = 1.63 and β (2)  d = 2.82.The Q and U templates at the reference frequency are obtained by scaling the dust intensity templates of [43] to 353 GHz and applying the polarization fraction and polarization angles of the d1 model.Gaussian small-scale fluctuations are generated similarly to the d1 model.
• d7: Simulations of thermal dust emission are based on a physical model of dust grains with a distribution of sizes, composition, and temperatures.Additionally, an environmental radiation strength field is assumed across the sky based on the dust temperature map obtained in [42].
The model also accounts for magnetic dipole radiation from iron grains [44], resulting in a SED that deviates from a simple MBB spectrum.
• d12: In this model, the dust emission is composed of six distinct emitting layers [29], so that the effective SED varies not only across the sky but also along the line-of-sight, leading to spectral decorrelation across frequencies.It represents the most complex dust model among the four considered in this study.Specifically, the maps of Stokes parameters at a frequency ν are given by the superposition of six templates scaled with a MBB law, each associated with a different map of spectral indices and temperatures: such that it matches the total dust emission of the Planck GNILC intensity and polarization maps at 353 GHz [45].
Galactic synchrotron emission is simulated using the PySM s1 model, which assumes a powerlaw SED in brightness temperature units [46] where β s is the synchrotron spectral index, varying with the position n in the sky.X s (n, ν s ) represents the synchrotron template at a reference frequency ν s .The chosen templates are the WMAP 9-year 23-GHz Q and U maps [47], smoothed with a Gaussian kernel of FWHM 5 • , with small-scale fluctuations added using the procedure outlined in [40].The synchrotron spectral index map is derived through a spectral fit to the Haslam 408-MHz data and the WMAP 23-GHz 7-year data [48].In summary, four distinct PICO data sets, hereafter referred to as d1, d4, d7, and d12, are created by combining the respective dust model maps with the PySM s1 synchrotron maps and random realizations of CMB.At each frequency, these components are first smoothed using a Gaussian beam with the corresponding FWHM of that PICO channel and then co-added to instrumental noise.Such resulting PICO frequency maps are ultimately brought to a common resolution of FWHM = 38 ′ and downgraded to a HEALPix resolution of N side = 512 before undergoing component separation analysis.For each data set, 50 simulations are generated, each featuring distinct realizations of CMB and instrumental noise.

Component separation methodology
In polarization, we can assume that data at different frequencies, d ν , are the superposition of three distinct types of components: where (assuming CMB temperature units) c(n) denotes frequency-independent CMB anisotropies, f ν (n) is the Galactic foreground emission, n ν (n) the instrumental noise, n the sky position, and ν = 1, . . ., N ν , with N ν the number of available frequency channels.An estimate ĉ(n) of CMB anisotropies can be formed via a weighted linear combination of the data: where the specific weights w ν assigned to data d ν across frequency channels are sought to minimize the variance due to foreground and noise contamination while preserving the CMB signal c(n).This is the strategy of ILC methods such as NILC [21].In this section, we seek to optimize these weights to minimize foreground variance rather than global variance as a function of sky regions and angular scales with our proposed ocMILC (optimized cMILC) pipeline.We first review the basics of the cMILC formalism (Section 3.1) which forms the skeleton of the proposed ocMILC pipeline, the specific properties of which are described afterward (Section 3.2).

Baseline cMILC approach
As demonstrated in [23], although the blind NILC method provides the CMB solution with minimum overall variance, it does not achieve minimum foreground variance.This minimum foreground variance can be attained through the more constrained semi-blind cMILC method, which involves the deprojection of certain foreground components.Achieving minimum foreground variance is essential to alleviate systematic biases in the recovered tensor-to-scalar ratio, especially in cases where residual foreground contamination after component separation surpasses the CMB and noise power at large angular scales.Such unfavourable outcome was observed in the case of NILC for the multi-layer dust scenario with PICO [30] and for less complex foreground skies in experimental configurations with narrower frequency coverage [31].Galactic emission exhibits three-dimensional characteristics [25], extending along the line of sight, and its spectral parameters vary spatially both along and across the lines of sight.When integrated along the line of sight and within the instrumental beam by CMB instruments, the superposition of emitting elements with varying spectral parameters results in spectral distortions of the baseline SEDs of the Galactic components.As demonstrated in the seminal work of [26], these foreground distortions can be accurately captured through a Taylor moment expansion of the baseline SEDs around pivot values of the spectral parameters.Foreground subtraction can then be enhanced by deprojecting certain moments of the Galactic foreground emission from the output CMB solution, a technique introduced by the cMILC method [23].
Referring to Equation 2.3, the Galactic synchrotron emission can be expanded around a baseline power-law SED up to second order in β s (n) as: where is a power-law with fixed pivot spectral index βs .Equation 3.3 highlights moments of the synchrotron spectral index, (β s (n) − βs ) k , with uniform SEDs ∂ k β s f sync (ν, βs ).Similarly, from Equation 2.1, the moment expansion of the thermal dust emission around a baseline MBB SED up to second order reads as: where is the MBB function for data in brightness temperature units, with fixed pivot values of the spectral index and temperature βd and Td .We assume the observed averages of spectral parameters across the sky as pivot values to ensure that most of the Galactic emission projects onto the zeroth-, first-, and second-order moments of the expansion: βs ≃ −3 [45,[48][49][50], βd ≃ 1.5, and Td ≃ 20 K [45].However, further refinement of the most suitable pivot values for spectral parameters could be achieved by considering the observed average in the sky area outside the Galactic mask used in the analysis or by requiring the cancellation of first-order moments in the average [51].
The cMILC weights w are built to deproject some of these foreground moments by adding to NILC nulling constraints against their homogeneous SEDs, while the remaining unconstrained contamination is left to blind variance minimization: where w = {w ν (n)} represents the vector of frequency-dependent, pixel-dependent cMILC weights, ⟩} denotes the data covariance matrix with the same dependencies as the weights, )} are the frequency-dependent but pixel-independent dust and synchrotron moment SEDs, and a CMB = {a CMB,ν } stands for the CMB SED, which is independent of frequency in CMB temperature units.The number of nulling constraints on moments in Equation 3.7 is only limited by the number of available frequency channels and the sensitivity of the experiment.
By collecting the SED vectors of the CMB and selected foreground moments into the mixing matrix and the deprojection coefficients into the vector: the cMILC weights can be expressed in the compact form: The cMILC method, initially introduced in [23], has since been applied in various subsequent works [31,[52][53][54][55][56].In such a framework, NILC represents simply a subcase of cMILC where only the first element of e and first column of A are employed in Equation 3.9.
Adding constraints to null some of the foreground moments in Equation 3.7 and, thereby, remove systematic biases generally results in an increase of the variance of the noise and unconstrained foreground components, such as non-included higher-order moments.There is thus a trade-off between deprojection and blind variance minimization, where components that significantly contribute to the overall variance in the data may be better suited for variance minimization rather than deprojection.However, in the baseline cMILC approach, there is no inherent procedure for selecting which moments should be prioritized for deprojection and which should be left for variance minimization, nor is there a method for assessing the optimal number of moments to deproject before incurring significant noise penalties.Moreover, if the effective SED of foreground emission turns out to highly deviate from the adopted baseline SEDs of Equations 3.3 and 3.5, moment expansion may not fully capture such distortions, in which case partial deprojection or variance minimization could be more effective.Therefore, there is room for optimization of moment deprojection in cMILC, as anticipated and preliminarily explored in [23].

Optimized cMILC (ocMILC) method
Our proposed ocMILC (optimized cMILC) method, an enhanced version of cMILC, aims to optimize moment deprojection and variance minimization in four ways: (i) allowing the number of deprojected moments to vary across sky regions and multipole ranges via a blind diagnosis of local foreground complexity; (ii) selecting moments for deprojection independently of the natural order of the Taylor expansion; (iii) allowing pivot values to vary across sky regions and multipole ranges to reduce the significance of higher-order moments in the moment expansion; (iv) partially filtering out the contribution of specific moments by modifying the deprojection vector in Equation 3.8 as follows: with |ϵ i | < 1.This last adjustment aims to strike a balance between reducing bias and increasing noise variance, drawing inspiration from [57].The method is fully data-driven.The first stage in the ocMILC pipeline aims at assessing the optimal number of moments to deproject across the sky and angular scales.This is achieved by leveraging information from byproducts of the Generalized Needlet ILC (GNILC) method [32,58], which allows for a blind diagnosis of the effective foreground complexity in the data by determining how many independent modes of the Galactic emission exhibit greater power than noise across the sky and angular scales.The GNILC diagnosis tool which is implemented in ocMILC is detailed in Section 3.2.1, while subsequent optimization layers of the ocMILC pipeline are described in Section 3.2.2.

Diagnosis of foreground complexity
The complex properties of the polarized Galactic foregrounds at microwave frequencies are not yet fully understood at the anticipated sensitivity levels of future CMB experiments.Furthermore, the integration of foreground emission along the line-of-sight and within the beam by CMB instruments leads to spatial averaging of each foreground spectral parameter, which results in spectral distortions of the expected foreground SEDs [26,59] and spectral decorrelation of the foreground emission across frequencies [25,27].Additional spectral averaging occurs during data processing steps such as spherical harmonic transforms.All these processes effectively augment the list of expected foreground parameters, adding further complexity to the modelling of the foregrounds.To avoid making incorrect assumptions about the number of foregrounds and their spectral properties, it is useful to adopt a data-driven approach which can inform us on the actual complexity of foregrounds in sky observations.The GNILC method [32,58] offers an avenue, by providing maps of the effective number of independent foreground components across the sky and different angular scales as byproducts, based on the experiment's sensitivity and frequency coverage.These GNILC byproduct maps, revealing the effective dimension of the foreground subspace, can in turn inform the cMILC method [23] about the optimal number of nulling constraints or moments to deproject without incurring significant noise penalty, depending on the sky area and the range of angular scales.Besides, parametric component separation methods may also benefit from these diagnostic maps to gauge the expected volume of the foreground parameter space depending on the sky regions.
Instead of making spectral assumptions, GNILC relies on statistical decorrelation to estimate the number of independent foreground degrees of freedom and perform component separation.It capitalizes on the fact that while the noise is largely uncorrelated across frequencies, the foreground emission, f ν (n), remains strongly correlated, such that it can be decomposed, at any frequency ν, onto a common basis of a few independent (i.e.uncorrelated, not physical) templates t where the number m of independent templates is typically smaller than the number of frequency channels.The focus of ocMILC is on mapping this number m rather than the templates themselves to provide insight into the complexity of the foreground sky.From Equation 3.1, the elements of the N ν × N ν data covariance matrix for frequency pairs (ν, ν ′ ) read as where C CMB , F, and N denote respectively the CMB, foreground, and noise covariance matrices.In Equation 3.12 we fairly assume that the three components are statistically uncorrelated.In practice, the ensemble average of Equation 3.12 is substituted with a sample average: where W (n, n′ ) are the weights associated to a Gaussian convolution kernel centred at n, with a standard deviation large enough to sample an adequate number of modes.
Considering that our focus is on Galactic foregrounds, we rewrite Equation 3.12 as where C N is the nuisance covariance matrix collecting the contribution from CMB and instrumental noise covariance matrices.Unlike the foreground covariance matrix, the nuisance covariance matrix can be estimated from simulations based on the knowledge of the instrument and the CMB statistics.The estimated nuisance covariance matrix, ĈN (n), is then used as an input of the GNILC pipeline to compute the whitened data covariance matrix: where the last equality follows from Equation 3.14 and is exact if the estimated nuisance covariance matrix perfectly captures the statistics of the CMB and noise components, and 1 denotes the identity matrix.Therefore, the eigendecomposition and diagonalization of the whitened data covariance matrix must have the following structure: where all eigenvalues significantly greater than one correspond to the independent eigenmodes of the Galactic emission which locally dominate over the variance of the noise and CMB components.Therefore, in each point of the sky n, we are able to estimate the number of independent Galactic foreground modes, m(n), above the noise, given the sensitivity of the considered experiment, through this model-independent approach where we only need to assume the two-point statistics of the noise and CMB components in the polarization data.The m independent foreground eigenmodes are simply collected in the N ν × m matrix U F .
In the language of principal component analysis (PCA), the template maps of Equation 3.11 are given by t and correspond to the m principal components detected above the noise along the line-of-sight n.They form uncorrelated, not physical, foreground templates, as indicated by their diagonal covariance matrix.If the foreground emission was completely uncorrelated across frequencies, as instrumental noise is, the method would return m = N ν , i.e. as many independent templates as the number of available frequency channels.In contrast, the greater the correlation among physical foregrounds across frequencies and with each other, the lower the value of m will be.In addition, when noise dominates the local data, the effective number m of independent foreground components being detected is also reduced.
However, foreground emission models with same overall amplitude, i.e. equal 'foregroundto-noise' ratio, yet different levels of decorrelation, can lead to a different number of independent components, thus holding different degrees of complexity and producing different GNILC outcomes.
We note that, in the above treatment, we considered the CMB as a nuisance whose covariance matrix can be predicted.This may not be exactly the case at large angular scales for CMB B-modes since the value of the tensor-to-scalar ratio r is unknown.However, this issue can be addressed in ocMILC by considering only the noise in the nuisance covariance matrix (C N = N), in which case our method outputs the number of independent modes from both foregrounds and CMB that dominate over the noise: with m CMB (n) = 0 or 1, as the CMB is fully correlated between frequencies.Given the sensitivity of future CMB experiments, we expect the lensing CMB B-mode signal to dominate over instrumental noise across the range of observed angular scales (see e.g.figure 1 of [30]).Therefore, The GNILC approach implemented in ocMILC thus allows one to perform a blind diagnosis of the complexity of the Galactic foreground emission in the observed sky and to possibly compare it with that expected from various foreground models.
We expect the 'foregrounds-to-noise' ratio, as traced by the whitened covariance matrix C(n) (Equations 3.15 and 3.16), to vary not only across the sky but also across angular scales.Therefore, most likely, m fgds will be higher at large angular scales where diffuse Galactic emission is brighter than CMB and noise by orders of magnitude, while it will be lower at small angular scales where the noise dominates the covariance budget.This reasoning supports an implementation of our diagnostic tool in wavelet domain, using a specific frame of spherical wavelets known as needlets [60,61].In such a framework, the data d ν (n) are convolved by a set of needlet bandpass filters in harmonic domain in order to form needlet coefficient maps, β j ν (n), including sky emission from specific ranges of angular scales: where a ℓm,ν are the spherical harmonic coefficients of the frequency maps d ν (n), b j ℓ is the j th needlet bandpass, and Y ℓm (n) are spherical harmonics.By construction, the needlet bandpasses satisfy j (b j ℓ ) 2 = 1 for all multipoles, so that the inverse needlet transform retrieves the initial data.
The methodology presented above is applied to the multifrequency needlet coefficient maps, β j ν (n), instead of the original sky maps, d ν (n).The data and nuisance covariance matrices, C j (n) and C j N (n) of Equation 3.14, are thus computed separately for each range of angular scales selected by the needlet bands.This implementation allows us to produce maps of the number of independent foreground modes, m j fgds (n), across the sky and in different ranges of angular scales.By diagnosing the number of independent modes of the Galactic foreground emission above the noise level, the blind approach outlined in this section can be leveraged to enhance the optimization of CMB component separation.In ocMILC, we exploit this information to optimize the selection of deprojected moments of the Galactic foreground emission within the framework of the semi-blind cMILC component separation method.

Optimization of moment deprojection
For each range of angular scales (needlet scale j), the GNILC diagnosis conducted in ocMILC (Section 3.2.1)provides a partitioning of the sky into clusters of uniform dimension, m j fgds , of the expected foreground subspace.Each cluster of the partition, as illustrated later in Section 4.2.1, offers constraints on the optimal number of moments to deproject in that region of the sky and range of angular scales without incurring significant noise penalty.In ocMILC, this information is fully leveraged, and only m j fgds (n) moments are filtered out at pixel n and needlet scale j.
Once the optimal number of moments to deproject is determined across the sky and needlet scales, ocMILC selects the optimal configuration of m j fgds moments, pivot values, and deprojection coefficients among the following sets of 1. moment SEDs: A distinct deprojection coefficient is assigned to each selected moment.The optimal combination of moments, pivots, and deprojection coefficients is determined so as to minimize the cMILC output variance corrected for residual noise variance in each GNILC cluster: where C represents the data covariance matrix and N denotes the instrumental noise covariance matrix.
The GNILC partitioning of the sky into clusters of pixels characterized by uniform foreground complexity is used to recalculate the local covariance matrix of the data and ocMILC weights in each cluster independently, drawing inspiration from the multiclustering strategy of [24].The covariance matrices in Equation 3.24 and the optimization procedure are thus computed within each independent GNILC region, representing a connected set of pixels sharing the same value of m j fgds .Therefore, the ocMILC CMB map is the solution with minimum local foreground variance obtained by partially (or fully) deprojecting m j fgds moments with varying pivot parameters across the sky and angular scales.

Data analysis and results
In this section, we analyse the simulated PICO d1, d4, d7, and d12 data sets described in Section 2, representing four different foreground scenarios, using the ocMILC pipeline.The performance of ocMILC is compared against that of the baseline cMILC approach and the fully blind NILC approach.
The data are decomposed into needlet coefficient maps (Equation 3.19) for analysis by the NILC, cMILC, and ocMILC pipelines.In this work, we employ cosine-shaped needlet bandpasses defined as where ℓ peak = 0, 50, 100, 200, 300 for j = 1, . . ., 5.These needlet bandpasses in harmonic domain are plotted in Figure 1, with each of them covering a specific range of multipoles.Since the needlet bands have compact support in harmonic domain, the needlet coefficient maps, β j ν (n), are obtained at a HEALPix N side equal to the smallest power of 2 larger than ℓ j max /2, where ℓ j max is the maximum multipole of the j th needlet band.
The outcomes from the pipelines are first inspected by means of the angular power spectra of foreground and noise residuals, calculated by applying the weights derived from the specified data set to the foreground and noise inputs of each simulation.Such spectra are compared against the primordial CMB B-mode signal across a range of tensor-to-scalar ratio values, r ∈ [0.4 × 10 −3 , 4 × 10 −3 ], where r = 0.4 × 10 −3 represents the actual upper bound at 95% CL targeted by PICO for its baseline instrumental configuration, as adopted in this work (refer to Section 2 for details on rescaling PICO CBE to baseline sensitivity).Throughout the paper, to alleviate residual foreground contamination along the Galactic plane, the angular power spectra for d1, d4, and d7 are computed over f sky = 70% of the sky using the publicly available Planck GAL70 mask applied to the maps.For d12, a more conservative masking strategy is employed, by extending the Planck GAL70 mask through smoothing the average d12 foreground residuals with a 5 • beam and thresholding them to achieve a final sky fraction of f sky = 40%.Although the d12 mask is constructed from foreground residuals that may not be accessible in real data analysis, it is employed solely for comparing outcomes from different pipelines within the same clean 40% region of the sky.The angular power spectra are computed using the pymaster Python package2 of the MASTER estimator [62,63], which facilitates proper deconvolution of masking and beam smoothing effects [64].
To further assess the quality of the CMB B-mode reconstruction, we examine the impact of residuals on the estimation of the tensor-to-scalar ratio, r.We estimate the tensor-to-scalar ratio from the recovered B-mode power spectra using the inverse-Wishart log-likelihood [65]: is the primordial CMB B-mode power spectrum for tensor perturbations with r = 1.We adopt A lens = 0.27 to align with the expected 73% delensing capability of PICO [30].As our simulations do not include primordial tensor perturbations, the best-fit value of r is only driven by residual foreground contamination, while its uncertainty is influenced by residual lensing, noise, and foregrounds.While the likelihood function in Equation 4.2 is not exact for a masked sky as it neglects correlations between different multipoles induced by masking, we ensure that these offdiagonal terms of the covariance matrix are negligible by binning the angular power spectrum with a constant binning scheme of ∆ℓ = 10.
We begin by discussing some limitations in the performance of the baseline cMILC approach (Section 4.1), then present the results obtained using the ocMILC method for all foreground scenarios (Section 4.2).This includes showcasing diagnostic maps illustrating the actual foreground complexity in each dataset (Section 4.2.1) and then presenting ocMILC outcomes related to CMB B-mode reconstruction and tensor-to-scalar ratio estimation (Section 4.2.2).

cMILC results and limitations
To illustrate some limitations of the baseline cMILC approach, we applied both NILC and cMILC pipelines to different PICO data sets as discussed in Section 2, using the needlet frame implementation presented above.We implemented three versions of cMILC: 1. cMILC-0, where all zeroth-order moments are deprojected.
In all cases, the pivot values of the spectral parameters were held fixed to their average across the sky in the PySM d1s1 model: βs = −3, βd = 1.54, and Td = 19.6K.The results for PICO data sets d1 and d7 are reported in Figure 2, showing the angular power spectra of foreground and noise residuals.Residuals from the most constrained version cMILC-012 are exclusively shown for the d1 data set.The results from Figure 2 lead to several conclusions.For the PICO d1 data set, the cMILC foreground residuals are lower than those of NILC at intermediate and small angular scales.However, on the largest angular scales, the averaging processes significantly distort the foreground SEDs, and the rigid cMILC deprojection fails to adequately suppress these low-ℓ distortions of the Galactic emission compared to simple variance minimization.As previously demonstrated in [23], imposing the deprojection of more moments in Equation 3.7 results in a penalty to the residual noise budget.This is particularly noticeable when examining the residuals from the application of cMILC-012 (where all zeroth-, first-, and second-order moments are deprojected) to the PICO d1 data set: the CMB solution, derived from the linear combination of 21 frequency channels with nine fully deprojected moments, suffers from a significant increase in noise contamination.Therefore, a compromise in the optimal number of moments to deproject must be found.Moreover, in this baseline cMILC approach, the moments have been projected out in the natural order of the expansion.However, there is no a priori reason to filter out lower-order moments first in the context of minimum-variance methodologies.It may be more effective to leave some major low-order contributions for blind variance minimization, as they are anyway prevailing in the budget of the observed data covariance matrix of Equation 3.12.When considering PICO maps with different dust models, such as d7 (right panel of Figure 2), the cMILC foreground residuals might even appear brighter compared to NILC in some cases (cMILC-01), indicating that the employed moment deprojection may not entirely capture the actual deviations of the dust model from the MBB SED.The application of ocMILC's optimized moment deprojection effectively overcomes all these limitations of the baseline cMILC approach, as demonstrated in the following sections.

ocMILC results
As described in Section 3.2, the ocMILC pipeline consists of two main steps: a blind GNILC diagnosis of the effective foreground complexity to determine the optimal number of moment deprojection constraints across the sky and angular scales, and the subsequent optimized moment deprojection.The outcomes from these two steps are presented in Sections 4.2.1 and 4.2.2, respectively.

Diagnostic maps of foreground complexity
We perform our GNILC diagnosis on the simulated PICO d1, d4, d7, and d12 data sets, each representing different foreground scenarios detailed in Section 2. In this analysis, the CMB signal is considered part of the nuisance budget of GNILC, allowing us to estimate the number, m fgds , of independent degrees of freedom associated with Galactic emission detectable by a telescope like PICO for each foreground scenario.We perform this diagnosis for both E-and B-modes, presenting results for the baseline realization n sim = 0, although similar conclusions can be drawn for other realizations.
Maps of m j fgds (n) for the fifth needlet scale ( j = 5) are shown in Figure 3 for the four distinct data sets.They provide the effective number of independent foreground components detected across the sky in both in E-and B-mode fields at this scale, thereby partitioning the sky into clusters of uniform foreground complexity.These diagnostic maps reveal different degrees of complexity across the sky depending on the underlying foreground scenario, even though all considered foreground models are consistent with Planck observations.In particular, d12 and d1 models exhibit a larger number of independent foreground modes detected across more extensive sky regions compared to d4 and d7 models, with indications of higher complexity for d12.These results hold for both E-and B-modes.The effective lower complexity of d4 with two MBB components for dust compared to the single-MBB d1 model, is attributed to the absence of spatial variations of each MBB spectral index in d4.Indeed, unlike the d1 model, where spatial variations of the spectral index induce higher-order moments, the constant spectral indices of the two MBB components in d4 results in the absence of moments beyond zeroth-order in Equation 3.5, thereby reducing the list of effective components contributing to the overall foreground emission.The highest level of complexity detected for d12 can be attributed to the multi-layer structure of the dust emission in this model.Also in this case, not all the coadded seven physical components are detected by the GNILC diagnosis, as some may be correlated or below the noise level.In summary, with high polarization sensitivity, as in the case of PICO, our approach offers a reliable tool for discriminating between different foreground models.Moreover, the observed decrease in the effective number of independent foreground modes towards high Galactic latitudes aligns with the spatial extent of Galactic emission in the sky and the prevalence of noise in the data at these latitudes.In E-modes, we also observe larger values of m fgds than in B-modes over larger regions of the sky.This discrepancy arises from the foreground E/B ratio, which exceeds 2. However, this trend is partially compensated by the even higher CMB E/B ratio, mitigating the discrepancy between E-and B-modes concerning the number of independent foreground modes.
In Table 1 we report the excursion set of m fgds observed in the four PICO data sets at the five different needlet scales.For each value, the corresponding sky fraction is also reported in parentheses.As expected, the effective number of independent degrees of freedom of the Galactic emission decreases towards smaller angular scales, where the foreground power decreases while the CMB and noise have relatively larger power.Differences among the four data sets are noticeable.The d1 and d12 data sets exhibit larger values of m fgds compared to d4 and d7 at the largest angular scales and over larger portions of the sky at intermediate and small angular scales, both in E-and B-modes.Compared to other data sets, the d12 data set, which incorporates multilayer dust emission, exhibits higher complexity in broader regions of the sky for all needlet scales.Again, in general, higher values of m fgds are found on larger portions of the sky for E-modes.At j = 2, corresponding to angular scales of ∼ 3 • , we find no sensible differences among the first three foreground models, as we reach the resolution limit of original Planck reference templates used in these PySM models.Conversely, at smaller angular scales, artificial fluctuations have been added to these templates through different power spectra extrapolations [40], thus returning again visible differences across the sky among the three models as diagnosed by GNILC.).In cases where multiple values are observed, the corresponding sky fraction (in percentage) associated with each value of m fgds is provided in parentheses.

CMB B-mode reconstruction
The clusters of pixels depicted in the diagnostic maps of Figure 3, each featuring distinct dimensions of the foreground subspace, guide ocMILC in determining the appropriate number of constraints for moment deprojection in each cluster at each needlet scale.This constitutes the first layer of optimization in ocMILC.The next three layers of optimization in ocMILC consist in selecting, for each specified number of constraints, the optimal combination of moments, pivot values, and deprojection coefficients to minimize residual foreground contamination in CMB B-mode reconstruction, as detailed in Section 3.2.2.The improvement induced by combining the various layers of optimization in ocMILC is demonstrated by comparing the power spectrum of ocMILC residuals with those obtained when optimization is exclusively applied to specific layers.This comparison is performed using the example of the PICO d7 data set (similar results are obtained for other models) and considering the following intermediate ocMILC pipelines: • ocMILC-v0, where optimization focuses on the number of moments, based on information from GNILC diagnostic maps, deprojecting the first m fgds moments in the natural order of the Taylor expansion (Equation 3.21); • ocMILC-v1, where optimization involves both the number and the set of moments, deprojecting the set of m fgds moments that lead to the lowest denoised output variance (Equation 3.24); • ocMILC-v2, where optimization extends to the number of moments, set of moments, and pivot values.
The comparison results, shown in Figure 4, reveal that deprojecting the first m fgds moments in the natural order of the Taylor expansion (ocMILC-v0) does not improve foreground subtraction across most scales, especially when the dust model, like d7, deviates substantially from a MBB.In such cases, nulling constraints result in a notable penalty in residual foreground variance attributed to higherorder moments.This suggests that the mitigation of foregrounds could have been better handled if some of the low-order moments had been left for variance minimization, while the higher-order moments had been deprojected.Optimizing the choice of moments (ocMILC-v1) indeed brings improvement, further enhanced by allowing the values of pivot parameters to vary freely (ocMILC-v2), enabling a more effective accommodation of the dust model through moment expansion.Enabling the weights to also partially deproject specific moments (ocMILC) provides additional flexibility to the blind variance minimization process, allowing for the mitigation of foreground modes not adequately captured by the moment expansion.This final refinement, in conjunction with the previous optimization steps, enables ocMILC to reconstruct the cleanest CMB solution among all considered intermediate pipelines for all data sets, achieving minimum foreground variance with limited noise penalty.This increased flexibility proves advantageous in effectively handling various thermal dust models through moment expansion.
To validate the accuracy of the GNILC diagnosis in determining the optimal number m fgds of moments to deproject before incurring significant noise penalties, we apply ocMILC to the PICO d1 data set while imposing the deprojection of m fgds − 1, m fgds + 1, and m fgds + 2 moments in each point of the sky and for each needlet scale.The resulting angular power spectra are compared with those of NILC and the proper ocMILC, with m fgds moments deprojected, in Figure 5.The outcomes clearly demonstrate that deprojecting one or two additional moments, beyond what is suggested by the GNILC diagnosis, leads to a significant increase in residual noise.Conversely, deprojecting m fgds − 1 moments results in larger Galactic contamination on intermediate and small scales.These results highlight the efficacy of the GNILC diagnosis in providing a data-driven assessment of the optimal number of moments for deprojection with ocMILC.The full ocMILC pipeline, integrating all layers of optimization, has been applied to all PICO B-mode data sets, encompassing diverse foreground scenarios outlined in Section 2. The angular power spectra of ocMILC residuals are compared with those of NILC for all foreground scenarios in Figure 6.For the d1, d4, and d7 models, ocMILC yields significantly lower foreground residuals across all angular scales compared to NILC, with only a moderate increase in noise contamination relative to NILC.These improvements in terms of residual foreground contamination, regardless of the foreground model, can be attributed to the several layers of optimization in ocMILC: multiclustering from GNILC, varying number and set of deprojected moments, varying pivot parameters, and varying deprojection coefficients, across the sky and angular scales.For the most challenging foreground model d12, which includes multi-layer dust components along the line-of-sight and spectral decorrelation [29], the reduction of residual foreground contamination is more modest than in other foreground scenarios but still outperforms NILC across angular scales overall.For all the considered data sets, the reduced residual foreground contamination of the ocMILC CMB B-mode solution is achieved with only a moderate noise penalty relative to NILC, thanks to the information provided by the GNILC diagnosis.These results showcase the flexibility of ocMILC in effectively handling deviations of thermal dust emission from a simple MBB through optimized moment deprojection.
The resulting r-posteriors for NILC and ocMILC across the four distinct foreground scenarios are displayed in the right panels of Figure 6.Biased estimates of the tensor-to-scalar ratio (r) are reported when the value r = 0 falls outside the 68% confidence level.Notably, ocMILC, in contrast to NILC, produces unbiased posteriors for r across all considered foreground scenarios even for such large fractions of sky.This bias removal is characterized by both a shift in the distribution's peak towards lower values of r due to reduced foreground contamination and an increase in the distribution's width due to increased noise.The biased r-distribution observed in NILC results from more Table 2. Optimal configuration of moments and dust pivot parameters for deprojection, as determined by ocMILC across various sky regions and needlet scales for each foreground data set.The optimal deprojection coefficient for each moment SED is indicated in parentheses.As detailed in Section 3.2.2, the optimal combination of these parameters is chosen by minimizing the denoised output variance of Equation 3.24.ter (sky region with uniform foreground complexity), the table includes information on the optimal number and combination of moments for deprojection, along with the corresponding optimal values of deprojection coefficients and dust pivot parameters.Examining these configurations allows us to draw conclusions consistent with the discussions in previous sections: • The optimal pivot values of the spectral parameters vary for different foreground models, sky regions, and ranges of angular scales.
• Higher-order moments are sometimes preferred for deprojection (e.g.d4s1), while the subtraction of lower-order components can be effectively handled by minimum variance.
• Most moments require only partial deprojection.
The ranges of m fgds values derived from the GNILC diagnosis for different Galactic models and angular scales, as reported in Table 2 and used by ocMILC, are similar to those shown in Table 1.However, minor differences exist in a few cases, particularly at large and intermediate scales and for small clusters.This discrepancy arises because ocMILC and Table 2 adopt the more conservative estimate m fgds = m GNILC − 1, where m GNILC is computed considering only instrumental noise as a nuisance component, and one mode is subtracted to account for the expected CMB contribution.In contrast, Table 1 reports m fgds = m GNILC , including CMB in the nuisance term.These minor variations result from the fact that, on large and intermediate angular scales, the CMB lensing Bmodes are not consistently detected across the entire sky as an independent component above the PICO noise level.This is expected, as lensing modes primarily contribute at small scales.However, these differences are negligible and do not impact the performance of the ocMILC technique.
Throughout this paper, we adopted the baseline instrumental configuration of PICO, consistent with required specifications detailed in the PICO report [6], which provides an overall sensitivity of 0.87 µK • arcmin.To complement our results, we present NILC and ocMILC outcomes for the PICO d1 data set, employing the alternative PICO configuration known as CBE (Current Best Estimate) [6,30], featuring an enhanced overall sensitivity of 0.61 µK • arcmin.The results are shown in Figure 8, showcasing the angular power spectra of residuals averaged over all realizations and corresponding posteriors of the estimated tensor-to-scalar ratio.Even with the increased sensitivity of the PICO CBE configuration, ocMILC still outperforms NILC, exhibiting lower foreground residuals across all angular scales with only a moderate rise in noise contamination.As anticipated, the noise residuals in both NILC and ocMILC CMB B-mode reconstructions are lower compared to those obtained with the baseline PICO configuration (displayed in the top panel of Figure 6).Consequently, this leads to lower upper bounds on the estimated tensor-to-scalar ratio (r < 0.00025 at 95% CL), as shown in the right panel of Figure 8. Furthermore, ocMILC improves upon the biased upper bound obtained with NILC for such a large sky fraction (70%).The ocMILC constraint aligns with the quoted forecast in [30] for PICO CBE sensitivities, albeit for a significantly larger sky fraction compared to [30].

Sensitivity to AME and synchrotron emission
The results in Section 4.2.2 demonstrated the capability of ocMILC to handle diverse dust models efficiently, thanks to the optimization of moment deprojection.This was achieved by integrating information from the GNILC diagnosis and optimizing the selection of moments, pivot values and deprojection coefficients across various sky regions and angular scales.In this section, our focus shifts to assessing the pipeline's robustness against varied models of Galactic synchrotron emission, while maintaining a fixed thermal dust emission model (d1) to isolate the impact of synchrotron model variations.Additionally, we investigate the pipeline's response to potential polarization of the anomalous microwave emission (AME).Specifically, in addition to s1, two alternative PySM synchrotron models are considered: 1. s5: This model features a power-law SED as defined in Equation 2.3, but the spectral index map from s1 has been rescaled to match the larger variability observed in S-PASS data [50].
In addition, small-scale non-gaussian anisotropies have been incorporated into the templates using the Logarithm of the Polarization Fraction Tensor (logpoltens) formalism.
2. s7: This model features a power-law SED with some curvature, where ν s = 23 GHz, and the X s (n, ν s ) and β s maps remain consistent with those in s5, while a spatially varying curvature term C s is introduced based on the analysis of intensity data observed in a patch of the ARCADE experiment [66].
As of now, there is no reported detection of AME polarization, and the most stringent upper bounds on the AME polarization fraction are at the level of 1% or lower [67][68][69][70], so that the actual contribution of AME to polarized emission still remains uncertain.Therefore, we explore three AME scenarios with different constant degrees of polarization across the sky: Π AME = 0.5, 1, and 2%.The AME model from PySM corresponds to the sum of two spinning dust components based on the Commander fit of Planck temperature data [42].The SED of the first component features a spatially anisotropic peak frequency, while the SED of the other component maintains a constant peak frequency across the sky.The polarization angles of AME are set to match those of thermal dust emission.
We begin by examining the outcomes of the GNILC diagnosis for the various synchrotron and AME models in Table 3, contrasting them with those of the d1s1 data set, where the synchrotron was represented by the s1 model and AME was assumed to be unpolarized.Clearly, the number of independent foreground B-modes detected by GNILC at each needlet scale and the associated sky areas tend to increase with more complex synchrotron models and higher degrees of AME polarization.The larger variability of the synchrotron spectral index in the s5 model compared to s1 is noticeable to GNILC, particularly at intermediate and small scales ( j = 2, 3 and 4), where larger regions of the , 3 (20.4), 4 (26.6)sky exhibit higher values of m fgds .In contrast, GNILC is not particularly sensitive to the spatially varying curvature of the s7 model (equation 4.3) for the frequency range and sensitivity of PICO, as the number of independent foreground modes detected remains similar for s5 and s7.Comparing the first lines of the top and bottom tables highlights that even a modest AME polarization fraction of Π AME = 0.5% is detected by GNILC in all needlet bands for the frequency range and sensitivity of PICO.This trend persists with increasing polarization fractions (Π AME = 1 and 2%), suggesting that higher degrees of AME polarization further enhance the effective foreground complexity across the sky and angular scales.
After evaluating the effects of more complex synchrotron and AME models in the GNILC diagnosis, we proceed to assess their impact on NILC and ocMILC CMB B-mode reconstruction in Figure 9.This is to determine whether there is a need to extend the optimization procedure of ocMILC to accommodate the higher complexity of synchrotron and AME polarized emissions.Indeed, despite ocMILC's optimization including synchrotron moments, the current pipeline maintains a fixed pivot value of βs = −3 for the synchrotron spectral index (Section 3.2.2).However, a higher variability of this parameter across the sky may necessitate optimizing the corresponding pivot values, similar to what was done for dust spectral parameters.Furthermore, the framework of moment expansion has not yet been extended to AME [26], so the current ocMILC pipeline does not attempt to deproject any AME moments, leaving them for blind variance minimization.
The angular power spectra of foreground and noise residuals in the NILC CMB B-mode reconstruction under various synchrotron and AME scenarios are displayed in the left panels of Figure 9. Having a more complex synchrotron emission (s5 and s7) results in a slight increase in foreground residuals across all angular scales up to ℓ ∼ 300.However, the impact of varying synchrotron complexity on NILC residuals is notably less pronounced compared to the impact of varying dust models, -25 -  as discussed in previous sections.Consequently, further optimization of ocMILC on the pivot values of the synchrotron spectral index is not expected to significantly enhance ocMILC results.Across all extensions of synchrotron emission considered, the ocMILC's CMB solution with a fixed pivot of βs = −3 indeed exhibits lower foreground contamination compared to NILC (top right panel), with percentage improvements similar to those obtained with the least complex model s1.In the case of the d1s5 data set, we even observe a 10% greater improvement of ocMILC over NILC at large angular scales compared to d1s1.Additionally, the increase in noise residuals in ocMILC remains consistent across all synchrotron scenarios.The negligible impact of synchrotron curvature is anticipated since, as demonstrated in [23,26], curvature is accounted for by the second-order term in the moment expansion (Equation 3.3), which is included in our optimized deprojection.These results show that, currently, with the available and foreseeable more complex synchrotron models, employing a fixed pivot for the synchrotron spectral index ( βs = −3) still yields highly effective moment deprojection with ocMILC.
Concerning AME scenarios, we observe a slight increase in NILC residuals for ℓ ≲ 250 as the AME polarization fraction increases (bottom left panel of Figure 9).This indicates that the presence of a polarized AME component affects the variance minimization, given the sensitivity and frequency coverage of PICO.Examining the ratio between ocMILC and NILC residuals (bottom right panel), we find that, while the improvement achieved by ocMILC at intermediate and small angular scales (ℓ ≳ 70) slightly degrades with increasing AME polarization compared to the case without polarized AME, it remains highly significant.However, for very low multipoles (ℓ ≲ 15), ocMILC foreground residuals become comparable to or larger than NILC foreground residuals depending on the degree of AME polarization.This highlights the necessity, in such scenarios and considering the frequency coverage of PICO, of incorporating AME moments into the ocMILC pipeline in the first needlet band.We leave this extension for future investigation, noting that there is currently no clear observational evidence for such levels of polarized AME emission across the sky.

Conclusions
The success of upcoming CMB experiments hinges on the ability to measure the primordial CMB Bmode polarization down to r ≲ 10 −3 without bias to validate the inflationary scenario.This requires exquisite control of instrumental systematics and Galactic foreground subtraction.While the blind NILC method enables the reconstruction of the CMB B-mode polarization field without encountering biases from foreground mismodelling that parametric methods may face, fully blind overall-variance minimization might not always be sufficient to reduce large-scale Galactic contamination below the primordial signal, especially for the most challenging foreground scenarios, depending on the available sensitivity, frequency channels, and the targeted r value.
Deprojecting moments of Galactic emission with the semi-blind cMILC method theoretically allows for minimizing foreground variance rather than overall variance.However, the associated noise penalty, including the increased variance of unconstrained foreground moments, does not always guarantee optimal outcomes for the most complex foreground scenarios.Specifically, deprojecting foreground moments in the natural order of the Taylor expansion may not be optimal, particularly for the moments that dominate the overall variance, which may be better suited for blind variance minimization.To address this trade-off, we introduced an optimized variant of cMILC, named ocMILC (Section 3.2), which refines the deprojection process through a data-driven approach by optimizing the number and selection of moments to deproject, as well as pivot values and deprojection coefficients.
In ocMILC, the optimization of the number of moments to deproject, depending on sky areas and angular scales, relies on the GNILC diagnosis of data complexity (Section 3.2.1),which maps out the number of independent observable modes of Galactic emission, m fgds , in the given data set across the sky and at different angular scales (Figure 3 and Table 1).This diagnosis informs ocMILC on the optimal number of moments to deproject locally before paying significant noise penalties.This strategy is corroborated by the results obtained in Figure 5.The specific subset of moments to deproject locally is also optimized to retain some dominant moments for blind variance minimization, while selectively deprojecting higher-order moments that are not predominantly captured by covariance estimation.The optimization of local pivot values of the dust spectral parameters across different sky regions and ranges of angular scales allows for accommodating local moment expansions associated with local averaging processes.Finally, allowing the deprojection coefficients ϵ k to be non-zero for partial deprojection of certain moments offers additional flexibility for limiting noise penalties.In addition, it allows to better handle foreground modes not properly described by the moment expansion of a single MBB SED.The optimal pivot values, subset of moments, and deprojection coefficients are determined to minimize the denoised output variance (Equation 3.24), and hence the residual foreground variance, in each GNILC cluster.
Testing on simulated PICO data demonstrated ocMILC's superiority over NILC in reducing Galactic contamination across all angular scales and all the considered thermal dust models over large sky fractions, while maintaining noise levels within the target sensitivity for the tensor-to-scalar ratio (Figures 6-7).Furthermore, ocMILC outperformed the standard cMILC, showcasing its efficacy in handling diverse dust emission models, including the multi-layer dust model d12 with decorrela-tion, while controlling noise penalties.While current synchrotron emission models do not necessitate pivot value optimization, the inclusion of mildly polarized AME emissions may require the incorporation of AME moments into ocMILC, a potential avenue for future exploration given the lack of observational evidence.Overall, ocMILC stands as a robust optimization for semi-blind CMB polarization reconstruction, offering promise for application to other satellite missions like LiteBIRD [5], with potential extensions to ground-based experiments [3,4] through appropriate E-to-B leakage corrections [71][72][73].

Figure 1 .
Figure 1.Cosine needlet bandpass windows in harmonic domain used in this work.

. 2 )
Here Ĉℓ = A lens C lens ℓ + C fgds ℓ + C noi ℓ represents the angular power spectrum of the output CMB Bmode map, averaged over all realizations, with C lens ℓ the lensing power spectrum, A lens its residual amplitude after some expected level of delensing, C fgds ℓ and C noi ℓ the foreground and noise residual power, while C ℓ (r) = r C r=1 ℓ + A lens C lens ℓ + C noi ℓ represents our theoretical power spectrum model, where C r=1 ℓ

Figure 2 .
Figure2.Angular power spectra of foreground (solid) and noise (dashed) residuals in B-modes resulting from the application of NILC (blue), cMILC-0 (orange), and cMILC-01 (green) to PICO d1 (left) and d7 (right) data sets.For the d1 case, results of cMILC-012 (red) are also presented.All the power spectra are binned with ∆ℓ = 5 and computed over 70% of the sky.The grey shaded area denotes the range of primordial tensor CMB B-mode power spectra corresponding to tensor-to-scalar ratio values r ∈ [0.0004, 0.004].

Figure 3 .
Figure 3. Maps of the detected number, m j fgds , of independent foreground E-mode components (left) and Bmode components (right) at needlet scale j = 5.The results are shown for four distinct foreground scenarios, d1 (first row), d4 (second row), d7 (third row), and d12 (fourth row), and they reveal varying degrees of complexity.These diagnostic maps, which partition the sky into clusters of uniform foreground complexity, were generated using the GNILC methodology outlined in Section 3.2.1,applied to full-sky frequency maps of PICO E-and B-modes.

Figure 4 .
Figure 4. Angular power spectra of foreground (solid) and noise (dashed) residuals over f sky = 70% of the sky resulting from the application of NILC (blue), ocMILC-v0 (cyan), ocMILC-v1 (red), ocMILC-v2 (green), and ocMILC (magenta) to the PICO d7 B-mode data set.This illustrates the progressive improvement in residuals achieved by combining several layers of optimization in ocMILC.

Figure 5 .
Figure 5. Angular power spectra of foreground (solid) and noise (dashed) residuals over f sky = 70% of the sky obtained from the application of NILC (blue) and ocMILC (others) to the PICO d1 B-mode data set, with variations in the number of moments deprojected compared to the number m fgds diagnosed by GNILC.The results of ocMILC are shown for four different cases, where, respectively, m fgds (magenta), m fgds − 1 (cyan), m fgds + 1 (red), and m fgds + 2 (green) foreground moments are deprojected across the sky and needlet scales.

Figure 6 .
Figure 6.Left: Angular power spectra of foreground (solid) and noise (dashed) B-mode residuals obtained from the application of NILC (blue) and ocMILC (magenta) to PICO d1 (top), d4 (second row), d7 (third row), and d12 (bottom) data sets.Power spectra have been computed over f sky = 70% of the sky for models d1, d4, and d7, and over f sky = 40% of the sky for model d12.Right: Posterior distributions of the recovered tensor-to-scalar ratio for both pipelines.In all plots, the grey shaded area represents the range of values of the tensor-to-scalar ratio: r ∈ [0.0004, 0.004].

Figure 8 .
Figure 8. Same as the top panels of Figure6, but for PICO CBE sensitivities.In both panels, the grey shaded area represents the range of values of the tensor-to-scalar ratio: r ∈ [0.0002, 0.004].The lower bound (r = 0.0002) corresponds to the targeted sensitivity at 95% CL for PICO CBE configuration, while the upper bound (r = 0.004) roughly corresponds to the energy scale of Starobinski's inflation model.

Figure 9 .
Figure 9. Left: Angular power spectra of foreground (solid) and noise (dashed) B-mode residuals over 70% of the sky resulting from NILC applied to the PICO d1 data set with varying synchrotron models (top) or varying AME polarization models (bottom).Right: Ratio of the power spectra between ocMILC and NILC for foreground (solid) and noise (dashed) B-mode residuals across different synchrotron and AME scenarios considered.

Table 1 .
Number m fgds of independent modes characterizing Galactic foreground emission, as detected by GNILC within each needlet band, for the PICO baseline configuration and four distinct foreground scenarios.The results are presented separately for E-modes (top table) and B-modes (bottom table

Table 3 .
Number m fgds of independent modes characterizing Galactic foreground emission, as detected by the GNILC in each needlet band for the baseline PICO configuration, as a function of the synchrotron (top panel) and AME (bottom panel) scenarios.In cases where multiple values are observed, parentheses indicate the corresponding sky fraction in percentage associated with each value of m fgds .