We implement a simple main beam correction in the maximumlikelihood parametric component separation approach which allows on accounting for different beamwidths of input maps at different frequencies without any preprocessing. We validate the approach on fullsky and cutsky simulations and discuss the importance and impact of the assumptions and simplifications. We find that in the cases when the underlying sky model is indeed parametric the method successfully recovers component spectral parameters and component maps at the predefined resolution. The improvement on the precision of the estimated spectral parameters is found to be minor due to the redness of the foreground angular spectra however the method is potentially more accurate in particular if the foreground properties display strong spatial variability as it does not assume commutation of the beam smoothing and mixing matrix operators. The method permits a reconstruction of the CMB map with a resolution significantly superior to that of the lowest resolution map used in the analysis and with the nearly optimal noise level facilitating exploitation of the cosmological information contained on angular scales which would be otherwise inaccessible. The method preserves all the advantages of a pixeldomain implementation of the parametric approach and as it deals with the beams in the harmonic domain it can also straightforwardly account for spatially stationary mapdomain noise correlations.
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Main beam treatment in parametric pixelbased component separation
{'Main beam treatment in parametric, pixel-based component separation': 'Arianna Rizzieri 1 , 2 , ∗ Josquin Errard 1 , and Radek Stompor 1 , 2  \n1 Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France and 2 CNRS-UCB International Research Laboratory, Centre Pierre Binétruy, IRL2007, CPB-IN2P3, Berkeley, US (Dated: September 20, 2024)  \nWe implement a simple, main beam correction in the maximum-likelihood, parametric component separation approach, which allows on accounting for different beamwidths of input maps at different frequencies without any preprocessing. We validate the approach on full-sky and cut-sky simulations and discuss the importance and impact of the assumptions and simplifications. We find that, in the cases when the underlying sky model is indeed parametric, the method successfully recovers component spectral parameters and component maps at the pre-defined resolution. The improvement on the precision of the estimated spectral parameters is found to be minor due to the redness of the foreground angular spectra, however the method is potentially more accurate, in particular if the foreground properties display strong, spatial variability, as it does not assume commutation of the beam smoothing and mixing matrix operators. The method permits a reconstruction of the CMB map with a resolution significantly superior to that of the lowest resolution map used in the analysis and with the nearly optimal noise level, facilitating exploitation of the cosmological information contained on angular scales, which would be otherwise inaccessible. The method preserves all the advantages of a pixel-domain implementation of the parametric approach, and, as it deals with the beams in the harmonic domain, it can also straightforwardly account for spatially stationary map-domain noise correlations.', 'I. Introduction': 'Detection and characterization of the primordial B -mode signatures in the polarization of the cosmic microwave background anisotropies is one of the major goals of current research in cosmology as it offers potentially groundbreaking insights into cosmology and fundamental physics. A major effort is therefore on-going in the field aimed at achieving this goal. In this context the large angular scales are both particularly interesting and challenging. The interest is due to the fact that some of the most popular models of the origins of the observable Universe tend to predict large angular scale B -mode polarization due to primordial gravitational waves they generically produce. The challenges are related to the smallness of the expected signals of interest, as compared to signals of other origins. Among those astrophysical contaminations, affecting polarization observations throughout the microwave band, stand out as one of the major obstacles. These are generically called foregrounds and are mainly due to galactic emissions either by dust particles aligned with the magnetic field of our Galaxy or by free electrons spiraling in that same field, referred to hereafter as dust and synchrotron signals, respectively. Removal of these contributions constitutes one of the major steps down the data analysis pipeline of any modern CMB experiment. The corresponding procedures are referred to as component separation. A number of different component separation approaches have been proposed and validated on simulated and actual data, e.g. [57, 10, 12-14, 19-21]. The different methods make different assumptions concerning the properties of the underlying signals, and include approaches which make no explicit assumptions about the frequency dependence of the foreground components, generically called non-parametric methods, as well as methods where the frequency scaling of all components is assumed to be known down to a few parameters, as in the parametric approaches. A successful approach has  \nto be also sufficiently flexible to include instrumental effects, which while originated from the instrument may affect the outcome of the component separation procedures potentially biasing the obtained results [4, 18, 24]. In this work we focus on accounting on the frequency dependent experimental beams in the context of the parametric component separation technique as introduced in the next section. Specifically, we propose a simple approach accounting for main beam differences between different single frequency maps, constituting the input data for the component separation step, directly in the component separation procedure. We demonstrate the approach on nearly full-sky and cut sky cases with characteristics motivated by forthcoming LiteBIRD [9] and Simons Observatory ( SO ) surveys [2], respectively. We show that the proposed technique can result in a significantly better resolution of the derived maps of the components, while achieving nearly optimal noise levels and being more robust with respect to systematic biases due to spatial variability of the foreground properties. These improvements come however at increased computational cost, we introduce and discuss novel, specialized numerical algorithms which make the implementation of the method feasible.', 'A. Parametric component separation formalism': 'We assume hereafter that the data collected by the experiment can be cast as a set of sky maps, each reflecting the signal as seen at different frequency and with well-characterized uncertainties. This is indeed the case of the scanning experiments, a category which includes the majority of on-going and forthcoming experimental efforts. We will also assume that the signal measured at each frequency can be represented as a linear mixture of sky templates each corresponding to a signal due to a different physical source, and referred to as component, as observed at some predefined fiducial frequency. This results in the following model of the input data,  \nwhere d is the data vector combining all the maps at each frequency band and for each Stokes parameter, A is the component mixing matrix, and s stands for the templates of the actual signals for each component at the fiducial frequency. n is statistical uncertainty assigned to the data assumed to be Gaussian-distributed and due to the instrumental noise. The objective of the component separation procedure is to provide reliable estimates of all (or at least some) sky components, s . The challenge is due to the fact that neither the mixing matrix, A , nor the component templates are known a priori and have to be both derived from the noisy data.  \nIn the parametric methods, the frequency dependence of the component amplitudes, as encoded by the columns of the mixing matrix, A , is assumed to be known down to a few unknown parameters, β i . For the parametric methods we therefore have A = A ( β i ). These parameters are called spectral parameters and are derived from the data as part of the parametric component separation procedure. Many existing parametric component separation techniques are based on the maximum likelihood principle, e.g. [10, 14, 21] .  \nThe quantities appearing in Eq. (1) can be defined in different domains, such as pixel, harmonic, or wavelet. The choice of domain will in general depend on the specific context, defined by the expected properties of the foregrounds, as well as the specificity of the chosen component separation method. The parametric approach can be implemented in any of these domains. Among these, the pixel domain offers certain advantages, as this is the domain in which accounting on the inhomogeneous noise and spatial dependence of the foreground properties is the most straightforward. The former feature is common to essentially all experiments and is related to the non-uniform scanning of the sky, and the latter is well confirmed observationally and particularly relevant for experiments targeting large patches of the sky as for example is the case of satellite missions. This choice has however a number of downsides, one of which is accounting for instrumental beams. We discuss this problem in detail below.  \nThe consideration presented in this paper are fully general and applicable to any parametric component separation method. However, for concreteness, we follow the implementation of [21] and assume that the procedure is implemented in two steps. The first, non-linear step renders the estimates of the spectral parameters, via a maximization of the socalled spectral likelihood,  \n-2ln L sp ( β i ) = const -d T N -1 A ( A T N -1 A ) -1 A T N -1 d , (2)  \nwhere N is the noise covariance matrix. The second, linear step estimates the sky components from the input data given the precomputed parameters as,  \ns = ( A T N -1 A ) -1 A T N -1 d . (3)', 'B. Accounting on beams': 'The instrumental beams generically characterize the response of the instrument to sky signals. A variety of issues related to instrumental beams are relevant in the data analysis of CMB experiments, first of which the modeling and characterization of the beams themselves. However, our goal here is to deal with beams at the map-level in the component  \nseparation step, where multiple frequency maps each with its specific, and potentially different beam, need to be processed simultaneously. We need to account for the different beams to avoid biases in the component separation procedure, and to determine the beam width of the recovered component maps.', '1. Basic data model': 'The instantaneous beams can be generally very complex, displaying both complicated spatial and frequency dependencies and resulting in complex, pixel-dependent effective beams in the recovered sky maps. While many of these effects may need to be accounted for in a complete analysis, our focus here is on the main property of the beams, namely, the dependence of their main lobe on frequency bands. We thus assume, hereafter, that such sufficiently accurate effective beams exist and that they are independent on the sky pixel location.  \nUnder these assumptions the input data set can be modeled as,  \nd = B true ( As ) + n (4)  \nwhere the operator B true represents the convolution with the true beams. This operator is a sparse matrix, which is blockdiagonal performing the convolutions for each frequency map separately. We note that while the matrices B and A appear in the data model in a similar manner, they act differently on the maps, with the mixing matrix, A , combining amplitudes of different components corresponding to the same pixel, and the beam matrix, B , combining pixels of the same map. In the case of a full sky coverage, the application of the beam operator to a sky map would correspond to a multiplication of their representations in the harmonic space. This observation is the basis of the (approximate) algorithm proposed below.', '2. Commutation of component mixing and beam smoothing': "Let us assume that the beams for all considered frequency maps are the same. In this case, we can wonder whether convolving the single frequency maps with the corresponding beam operator, B true , is equivalent to convolving first the component maps with the same beam and combining them together by applying the mixing matrix. Whenever this is the case, on defining ˆ B as the beam operator acting on the component maps, we can write,  \nBA = A ˆ B , (5)  \nas the beam smoothing and the component mixing operations can be swapped. We refer to such a case as a beam smoothing and mixing matrix commutation.  \nWhile for the commutation to be the case, the experimental beams have to be the same for all the channels, this is not always sufficient. Indeed, we have to also assume that the mixing matrix, A , does not have spatial variability, thus it is the same for all the pixels. However, even if A is spatially varying on scales comparable to the beam sizes we can still hope that a generalized commutation holds which reads,  \nBA = A ' ˆ B , (6)  \nwhere we have defined modified (effective) mixing matrix A ' accounting for the effect of the beam convolution. We note that while there are no generic prescriptions how A should be generalized to A ' , a possible approach has been suggested in [8]. In the following we do not consider this option, instead we discuss the implications of assuming the commutation as in Eq. (5) and evaluate the error this entails in cases of a spatially varying mixing matrix in cases of different scales and amplitudes of the spatial variability of the spectral parameters. We expect that the error will generally be larger the larger the beam considered and the more spatial variability is present in A . We note that the formalism presented below could be straightforwardly generalized to include the cases as in Eq. (6), if needed.", '3. Specialized data models': "There are inherently two different manners in which we can proceed starting from Eq. (4). The first one is to smooth all the single frequency maps to some common resolution prior to the component separation, as done in e.g. [9, 26]. In the following we refer to this method as common resolution approach and use it as a reference method to compare other methods against. On denoting the corresponding beam operator as B final , we transform the data set modeled as in Eq. (4), as follows,  \nd ' ≡ B final B -1 true d = B final As + B final B -1 true n = As final + B final B -1 true n , (7)  \nwhere the last equality assumes the commutation of the effective beam matrix and the mixing matrix from left to right and s final stands for the component templates all smoothed with the same final beams. The components in the new data set, d ' , are now combined together as defined by the mixing matrix without involving any beam operators, as assumed in the standard component separation techniques. The beam operations are now either hidden in the component templates, which are now smoothed with the effective beam, as well as in the properties of the map-level noise, which is now both enhanced and correlated. Both these effects would have to be taken into account on the first step of the component separation, as any mischaracterization of the noise may lead to significant biases in the estimated spectral parameters, Eq. (2). In the common beam approaches these are typically mitigated by first smoothing all the input maps to the common resolution as given by the largest experimental beam, i.e., usually that of the lowest frequency band, and subsequently downgrading the map resolution via increasing the pixel size. The pixel-pixel noise correlation induced by the beam smoothing are consequently suppressed and the effective noise can be considered (nearly) uncorrelated, if the final pixel size is sufficiently large. The downgraded data can be then described with the data model,  \nd ' = As final + n eff , (8)  \nwhere n eff differs from n as its level reflects the downgrading procedure and which, by construction, can be presumed to be uncorrelated. We emphasize that the smoothing/downgrading procedure has to be applied to all input,  \nsingle frequency maps consistently. We also assume implicitly that it commutes with the mixing matrix. Indeed, given that the downgrade has to be performed on scales larger than the final beam, this may set requirements on spatial variability of the mixing more demanding than the beam commutation itself as already assumed earlier. We also note that while it is possible to use these maps to perform the second step of the component separation procedure, Eq. (3), this may not be desirable as the resulting component maps, will have complex smoothing determined by the pixelization. Indeed, typically the mixing matrix estimates derived from the smoothed and underpixelized maps are instead used in Eq. (3) together with the smoothed but not underpixelized input maps. The resolution of the resulting component maps is then determined by the adopted common beam. This procedure assumes usually simple diagonal noises. This is justifiable given the linearity of this step ensuring that the penalty for using inappropriate noise characteristic is not a bias but rather sub-optimal statistical uncertainty. This however also implies again some implicit assumptions about the mixing matrix variability as we use its lower resolution estimates for the higher resolution maps. If this is justifiable then the map resolution used on this step could be even higher than that of the largest input beam and as high as the one for which this assumption holds. We discuss some of the aspects of this choice in the following noting here that in practice making the right choice can be rather difficult and objective-dependent. Indeed, it may well be that using the underpixelized, and oversmoothed component maps is the most robust way forward for cosmological parameter estimation, while higher resolution maps with more complex noise and potentially higher residuals can still be a useful source of information for other science goals.  \nAll of the above issues can be sidestepped if we explicitly account for the beams in the component separation problem. Indeed, let us start from the data model in Eq. (4) and introduce the target beam of the output component maps, hereafter denoted as B final . Hereafter, we implicitly assume that the pixel size is set as dictated by the smallest beam out of the ones considered in the problem. Algebraically, the final beam enters the data model as,  \nd = B true A ˆ B -1 final ˆ B final s + n = B true A ˆ B -1 final s final + n , (9)  \nwhere s final are the beam-smoothed component templates. While any value of B final can in principle be assumed, some values results in better numerical stability. We discuss this later on. We note that this data model makes no assumptions about the commutation between (any of) the beam operators and the mixing matrix.  \nIn the cases where the commutation of the final beam and the mixing matrix can be assumed, i.e.,  \nB -1 final A = A ˆ B -1 final (10)  \nthe data model can be further simplified to read,  \nd = B true B -1 final As final + n . (11)  \nWhile the validity of the commutation has to be evaluated for each experimental setup and spatial variability allowed for in the mixing matrix, see Section I B 2 for more  \ndetails, the advantage of this data model is that it permits on numerically more efficient implementation resulting in reduced computational cost with respect to the previous case. We also note that this data model is equivalent to the data model in Eq. (7), where the mixing matrix is kept in its original form but the noise covariance is made more complex, as opposed to the data model in Eq. (11), where it is the mixing matrix which is modified. This is however only strictly true if the commutation relation in Eq. (10) is fulfilled. We emphasize that whichever method is used the (de)convolution of the beam will always lead to some spurious effects whenever only partial sky is observed. Indeed, for any pixel located within the beam size from the boundary any such operation is bound to fail. In the following we discuss the relevance of this effect in the context of the proposed method. We also note that the beam correction as discussed here is similar to the treatment used in some sampling based component separation codes, e.g., COMMANDER [14]. These methods perform the component separation while characterizing simultaneously statistical properties of the components, such as their power spectra. Such applications can in principle avoid the boundary problem mentioned above as they estimate full sky samples conditioned on the observed, cut-sky data. The price to pay are the assumptions concerning the statistical properties of the sky components, such as Gaussianity and stationarity.", 'A. Framework': 'We employ the two-step parametric component separation method, Eqs. (2 and (3), [21], assuming that the input data model is as given in Eq. (9) or in Eq. (11). The formalism therefore remains the same and the only change is that of the mixing matrix, A , which is now replaced with its generalized version denoted as ˜ A which reads either,  \n˜ A ≡ B true AB -1 final , (12)  \nfor the data model in Eq. (9) or,  \n˜ A ≡ B true B -1 final A , (13)  \nfor the data model in Eq. (11).  \nThe component separation problem, Eqs. (2) and (3), is then solved by maximizing,  \n-2ln L sp ( β i ) = const (14) -( d T N -1 ˜ A )( ˜ A T N -1 ˜ A ) -1 ( ˜ A T N -1 d ) .  \nand estimating the component maps via,  \ns final = ( ˜ A T N -1 ˜ A ) -1 ˜ A T N -1 d . (15)  \nWhile conceptually simple, the need for accommodating the generalized mixing matrices calls for more involved numerical algorithms needed for the implementation of the method. We describe them in the following section.', '1. Algebraic operations': 'The algebraic operations involved in these two steps can be in general performed either in harmonic or pixel domain, each having its own advantages. In the harmonic domain applying the beams is more straightforward, at least for the full sky coverage and axially symmetric beams as each beam operator B is represented by a diagonal matrix. In the pixel domain instead we can easily account on the spatial variability of the foregrounds by allowing the mixing matrix to be pixel-dependent. The noise inhomogeneities are easier to treat in pixel domain, as is the general case of the non-skystationary correlated noise, the computational load needed for this latter case notwithstanding. The harmonic domain allows for easy treatment of the sky-stationary noise, what in some applications can be a sufficient approximation. Here, we choose an alternative approach which tries to capitalize on advantages of both these domains. Therefore, each matrix in Eqs. (14) and (15) is applied in its most suitable domain and the necessary back and forth transforms between the pixel to harmonic domain are performed using spherical harmonic transforms. Specifically, in our implementation the noise correlations matrices, as well as the component mixing matrix, A , are applied in the pixel domain. This allows us to deal rigorously with the cut sky, inhomogeneous noise, and the spatial variability of the foregrounds, effects that are naturally living in the pixel domain. A general (non-skystationary) correlated noise should also be treated preferentially in pixel domain though it could be computationally rather involved as it may in general require map-making type operations within the component separation code. We leave such an extension for future work.  \nAsky-stationary correlated noise can instead be accounted for most straightforwardly in the harmonic domain and treated together with the beams, as the beam operators are applied in the harmonic domain. Given that we do not have typically access to the full sky we can not recover actual harmonic representation of the sky maps. We therefore can deal instead only with the pseudo harmonic coefficients, which are formally a convolution of the observed sky mask and the true harmonic coefficients. While it is an approximation, it is usually a very good one as the modes of interest, i.e., those affected by the beam smoothing, are the ones much shorter than the observed patch size. In pixel domain, for well-behaved beams, we expect the failures to be well localized and limited to a narrow strip of pixels along the patch boundary and that they quickly vanish when we move away from it towards the patch center.  \nThere are three basic operations we need to be able to perform in order to calculate the likelihood and component signal estimates in Eqs. (14 and (15). These are:  \n- · noise weighting of frequency domain maps , N -1 m : in the current implementation this is performed directly in the pixel domain owing to the assumption of diagonal noise covariances. If needed the skystationary, correlated noise component can be easily added as appropriate noise weighting applied in the harmonic domain.\n- · generalized component mixing , ˜ As : this is performed  \nas a sequence of operations with the operators applied to pixel domain objects from right to left. All the beam operations pass through harmonic domain, while the actual mixing of the components, i.e., A , is applied directly in the pixel domain. This requires back and forth transformations between the pixel and harmonic domains. The number of the necessary transformations depends on the adopted data model, with the generalized mixing matrix in Eq (12) requiring twice as many of those as the mixing matrix in Eq. (13).', '· transposed generalized component mixing , ˜ A T m :': 'this operation is also implemented as a sequence of operations, analogous to the case above but performed  \nin the reverse order.  \nThese operations suffice for the computation of ˜ A T N -1 d . The computation of the full covariance matrix, ˜ A T N -1 ˜ A is more involved as in the case at hand it does not have simple, block diagonal structure with the blocks corresponding to different observed sky pixels. Thus in general the construction and the inversion of the full pixel-domain covariance is therefore only possible for the lowest resolution cases with a limited number of sky pixels. These problems can be overcome in a specific case of a pixel-independent mixing matrix, A , and homogeneous noise. We can then apply both the mixing matrix and the noise covariance matrix directly in the harmonic domain and represent all the map objects via their spherical harmonic representations. Then matrix ˜ A T N -1 ˜ A , which would then operate on harmonic coefficients, would be block-diagonal and thus simply invertible block-by-block.  \nIn cases of more general component mixing or noise patterns, such approaches are not applicable, and instead we develop an iterative solver to compute directly a product of the inverse noise covariance and a pixel-domain map-like vector.  \nIn our approach we use a preconditioned conjugate gradient (PCG) solver to solve the system,  \n˜ A T N -1 ˜ Ax = ˜ A T N -1 d , (16)  \nobtaining the current component map estimates, x , which we then insert back in the spectral likelihood,  \n-2ln L sp ( β i ) = const -( d T N -1 ˜ A ) x . (17)  \nThe implementation relies on the PCG solver of scipy [25]. The operations discussed earlier are sufficient for the implementation of all the computations required by the PCG solver. To accelerate it, we employ a preconditioner constructed, and applied, in pixel domain. It looks like the system matrix but with the actual beam free mixing matrix A , with the spectral parameters fixed at their current estimate: A T N -1 A . This preconditioner is able to account for the spatial variability of the foregrounds, but not the beams. Consequently this PCG solver converges rapidly to a solution at the largest angular scales, but needs many more iterations to converge the smaller angular scales, where the beams impact is more important and the preconditioner inefficient.  \nWe have also explored an alternative preconditioner which is expected to perform well if the mixing matrix variability across the sky is sufficiently small. The preconditioner is constructed as the inverse noise covariance as described above in  \nthe case of homogeneous noise and pixel-independent mixing matrix. It is thus built and applied directly in the harmonic domain where it has a block-diagonal structure and, by construction, can account nearly strictly on the beam effects. Due to sky variability of the mixing matrices considered in the examples studied in this work, this preconditioner was found to not be competitive with the one described earlier. Nevertheless, some hybrid constructions retaining the advantages of both could outperform any of the two. As could some more generic preconditioners, i.e., such as those based on two-level constructions [11, 22]. We leave such investigations to future work and use the first preconditioner as proposed here throughout this work.  \nOur implementation is based on the fgbuster code [16], but it accounts explicitly for the beams in the mixing matrices and uses the PCG solver to evaluate the spectral likelihood. The effective mixing matrix ˜ A is implemented as an operator, which applies it to any arbitrary pixel-domain vector without however ever constructing it. In the cases where the sky-stationary noise correlations are present, and described by the noise power spectrum, N /lscript , we first modify the beam function, b true ,/lscript → b true ,/lscript N -α /lscript , α = 0 . 5 or = 1 for the left, respectively right, hand side operators in Eq. (16) and then call the operator as usual. This operator is then used in all the generalized spectral likelihood calculations as well as in the PCG solver. The code is released as part of the fgbuster suite of codes, at [17].', '2. Likelihood optimization': 'The extended fgbuster package includes also all the functionality necessary for the maximization of the spectra likelihood what constitutes the first step of the component separation approach. It has been implemented and validated as part of the described work. Consequently, though this functionality has no direct bearing on the results presented hereafter, for completeness, we provide below a brief description of details of the numerical implementation.  \nThe likelihood maximization is performed with help of the optimize.minimize module of the scipy package [25]. This package offers various methods to reach the convergence. As long as we fit for a limited number of spectral parameters we were able to reach the convergence with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm computing the gradient of the spectral likelihood numerically. Still, this implementation was quite slow to converge, due to the fact that each evaluation of the spectral likelihood takes order of minutes to converge even for lower resolution maps, mostly driven by the number of PCG iterations required and the spherical harmonic transforms cost when considering higher resolution cases. To minimize the number of the spectral likelihood evaluations, we have implemented an analytical expression for the gradient of the spectral likelihood, given by,  \n-2 ∇ β ln L sp = -2 d T N -1 ˜ A ( ˜ A T N -1 ˜ A ) -1 ˜ A , T β Pd , (18)  \nwhere,  \nP ≡ N -1 -N -1 ˜ A ( ˜ A T N -1 ˜ A ) -1 ˜ A T N -1 . (19)  \nWe note that the only additional operations required in this computation involve an application of the derivative of the generalized mixing matrix with respect to the spectral parameters to a frequency map vector. This is straightforward and computationally cheap as it can be done without any additional spherical harmonic transforms, for instance, via generalizing the general mixing matrix operator so it performs simultaneously the product of the generalized mixing matrix and its derivative by an arbitrary (but common) pixeldomain vector. All the other operations, in particular those included in the PCG iterations, have to be computed anyway for the spectral likelihood evaluations and reused here. This allows us to use more efficient optimization algorithms, which require fewer computations of the likelihood function. This optimization is particularly important when we fit for many spectral parameters as in the cases with spatial variable foreground scalings.', '3. Numerical efficiency': 'As mentioned in the previous paragraph, the computational time of a single likelihood evaluation is essentially proportional to the number of iterations required for the PCG solver to converge. This emphasizes the need for very efficient preconditioners. In the next section we discuss the performance of the preconditioner proposed here.  \nFor a given PCG iteration, the computation time increases with the resolution of the maps, as it is dominated by the time needed to compute the spherical harmonic transforms. These need to be performed every time we apply the generalized mixing matrix or its transpose. At present this is done with help of the healpy package [15, 27], however, we currently explore a potential use of a JAX -based package (and GPUs). We anticipate that this will significantly alleviate the computational limitations of the current code and will be made available in the forthcoming release. All the results shown in this work have been produced with healpy -based code.', 'III. Results': 'We present in this section some results of the method presented above with the goal of comparing the validity and performances of the different data models given in Eqs. (4), (9), (11) among themselves and with the reference common resolution approach, Eq. (8), in realistic experimental scenarios. We consider in the following only beams given by a main lobe and with circular Gaussian profile, however any axially symmetric beam profile could be assumed instead if needed.  \nThe impact of the beam treatment in the component separation can be expected to depend on the specifics of the instrumental setup assumed, in particular on the range of resolution of the input frequency bands of a given instrument: the broader this range the more significant the impact will be. We therefore show the results for two different experimental setups, spanning the range of possible choices. We emphasize that the cases studied hereafter are highly idealized as our goal is to validate and evaluate the new beam  \ntreatment and not to explore yet all possible complications. We note, however, that many of the relevant effects, such as bandpasses, beam chromaticity, etc., which without a doubt will be important in actual applications, could be included in our formalism rather straightforwardly. We leave that to future work.', 'A. Full-sky case': "The first setup we consider is that of a full sky mission, with resolution of the input frequency bands going from 17 . 9 ' to 70 . 5 ' , as it would be the case for the LiteBIRD instrument [9]. In Table I we list all the specifics of this setup, which we will refer to as LiteBIRD -like setup.  \nTABLE I: Frequency, sensitivity and FWHM of Gaussian beams for the LiteBIRD -like setup. Taken from [9].  \nThe frequency maps are simulated at nside =512, and the component separation is performed on the full sky, with a galactic plane mask at f sky = 60% from [1] applied on the recovered components when necessary. To create the simulations we perform the following steps: first the sky signal is simulated through the PySM [23] package, the true beams of Table I are then applied, finally the noise is simulated as a Gaussian realisation and added to the convolved sky signal. As sky model we take the PySM c 1 model as a CMB realisation with lensing and no primordial B modes, to which we add thermal dust and synchrotron emissions. For these latter we consider in the following two options, first the PySM model d 0 s 0, which simulates the dust frequency scaling with a modified black body and the synchrotron one with a power law, assuming as spectral parameters the spectral index of the modified black body and its temperature, and the spectral index of the power law, respectively with values, β dust = 1 . 54, T dust = 20 K and β synch = -3 across the full sky. We then consider the d 1 s 1 model, again given by a modified black body and a power law but with spectral parameter values varying across the sky, with most of the variations along the Galactic plane. This latter model is adequate for our scope of exploring the effect of dealing with the beams in a scenario with spatially varying scaling laws, and a more complex sky model with a higher level of spatial variability would still suffer the complications of those in an enhanced way with respect to what presented below, according to its level of", '1. Without spatial variability': "We apply here the proposed approach to the simulations without spatial variability of the foreground properties and compare it with the common resolution approach. As the assumed scaling laws are pixel independent the beam and mixing operators commute. The expected differences in performance of different methods is therefore due to treatment of the map resolution and the noise. For the common resolution approach we perform the preprocessing by convolving each input frequency map with the effective beam B final B -1 true , where the first term is a Gaussian beam with FWHM given by the largest instrumental beam, 70 . 5 ' . In addition, for the sake of the first step of the component separation, we downgrade to the pixel resolution of nside =32, in order to suppress the noise correlations the additional smoothing unavoidably introduces. The value of B final in general is chosen in a way that the enhancement of the noise power due to the true beam deconvolution, B -1 true , is suppressed. As this is followed up by the underpixelization of the maps, there is some freedom in selecting the final beam size. Downgrading to such a low nside dramatically decreases the angular scales contributing to determining the spectral parameters. In the case of the pixel-independent foreground properties this may only lead to some loss of precision. While we could try using larger nside to compensate for this, this poses a potential risk due to a bias ensued by the remaining correlations in the noise, which are not easily incorporated in the noise model. While there may be some potential tradeoff to explore, we use, hereafter, nside =32 as the resolution of the downgraded map. The spectral parameter estimates recovered then on the first step of the component separation are subsequently used to produce the component maps on the second step. In this latter case, we can in principle keep the original pixel resolution of nside =512 and choose the final beam size accordingly. As the second step is linear, this choice affects the noise level of the produced maps, but does not result in biases, and the beam size can be then optimized accordingly. For simplicity, unless mentioned otherwise, we adapt here the value of 70 . 5 ' (FWHM) for B final , corresponding the largest experimental beam, to efficiently compensate for B -1 true . As far as the method proposed here is concerned, given that the beam and mixing operators commute the data models in Eqs. (9) and (11) are equivalent. Moreover, as the noise covariance is also pixel independent we can perform all the operations in the harmonic domain (including the spectral likelihood minimization), and come back in the pixel domain only once we have the recovered components.  \nThis simplifies the comparison as we have only two methods to compare, the reference, common resolution method and one of the variants of the beam-treating approach. We first compare the estimation of spectral parameters. We do so by simulating 1 , 000 input data sets with different instrumental noise realizations. The histograms are shown in Figure 1. We see that both methods perform similarly. This may seem surprising given the fact that the new method should be able to exploit all the small scale information not accessible to the oversmoothed, common resolution technique.  \nWhile in principle this is indeed so, we note that due to the redness of the power spectra of all non-CMB components the information concerning their spectral parameters residing at small-scales is essentially negligible. Instead the differences can be seen on the second step of the component separation, when the component maps are estimated. As it shown in Figure 2 the impact of accounting for the beams in the component separation can be indeed significant. Given the very simple foreground models assumed in this case, the CMBresiduals are dominated by the instrumental noise after component separation. We point out the different recovered components start to be dominated by the noise at different scales: at about /lscript = 150 for synchrotron, at about /lscript = 350 for dust in the common resolution approach and /lscript = 800 for dust in the beam-treating approach.  \nThe noise spectra after component separation are shown in Figure 3. In addition the different recovered components have the noise exploding at different scales, which depends on the beam convolution remaining in the noise after component separation and which is given by a combination of all the B -1 true , which varies for the different components, but is steeper in the common resolution approach then in the approach dealing with the beams in the component separation. The gain in resolution in the latter approach is thus due to fact that the noise starts to dominate later and when it does it goes up in a less steep way. The wider range of accessible scales with this method is clearly visible in the plot. Also, the gain in resolution is more important for the dust component, less for the synchrotron component. This is because the additional resolution that we can now exploit comes from the higher frequency bands with smaller beam FWHM, and thus affecting mainly the dust component. With the proposed approach about a factor 1.5 is gained in the recovered CMB residuals before the noise stops on being white and more than a factor two for the dust component. We show a patch of the maps of the recovered CMB residuals for one realization of the sky, thermal dust emission and synchrotron emission on the Galactic plane in Figure 4 and in a cleaner region of the sky in Figure 5. In both cases the improvement in resolution is visible in the approach dealing with the beams in the component separation, in particular for the dust component. The maps shown here are smoothed with a Gaussian beam with the width as recovered from the noise power spectra shown in Figure 3 suppressing the excess noise on the smallest scales and defining the resolution. While better resolution is in principle possible it will come at the cost of significant noise excess. The main goal of the proposed procedure is to facilitate access to the information on the smaller angular scales, what is perhaps the most important for Galactic science goals. In terms of constraining cosmological parameters the most significant impact is expected for parameters relevant for small scales. Indeed with this particular setup we do not get a significant improvement on the uncertainty on r . A mild improvement could instead be reached with other experimental setups, with a broader range of resolution of the input frequency bands or for experiments aiming at constraining r from medium scales, which is for example the case of SO or any other ground based experiment for which accessing the reionization peak in the B -mode power spectrum is excluded.  \nFIG. 1: Recovered spectral parameters from component separation on 1000 simulations with different noise realizations. The horizontal axis give the difference with the true spectral parameter values. In blue : common resolution approach, dealing with the beams before the component separation with B final = 70 . 5 ' and downgrading to nside =32; in red : proposed approach, dealing with the beams in the component separation. The two distributions are both peaked on the true spectral parameter values and have very close width, given the number of samples they are statistically consistent.  \n<!-- image -->", '2. With spatial variability': "For the foreground simulations we now consider the PySM d 1 s 1 model with spectral parameter variations across the sky. To have low CMB residuals from the component separation we are now obliged to account for the spatial variability of the mixing matrix.  \nWe first explore the effect of performing the commutation of the mixing matrix and the beam operator, as in Eq. (5), in the data preprocessing of Eq. (7). We do this by first computing the difference,  \nBAs -ABs , (20)  \nwhere as A we use the true d 1 s 1 spectral parameter values, s are the component maps at nside =512, and B a Gaussian beam with FWHM 70 . 5 ' . We find this difference to be non-negligible for all the frequency bands, except for the one at which we have normalized the mixing matrix, which is spatially constant in that band, Figure 6. The difference is particularly noticeable along the Galactic plane, as this corresponds to the sky region where the spectral parameters considered are more spatially varying, Figure 7. While at its face value this difference seems clearly important, of real interest is how it may affect the recovered components, and how much this depends on whether the Galactic plane is masked or not. We thus compute the residuals that we would have because of the beams in the case of noiseless input simulations and with the mixing matrix with one patch per pixel, at nside =64, and with B final = 70 . 5 ' . Doing so we ensure that there is no bias due to the choice of the scaling laws in the mixing matrix, and the residuals that we find are entirely due to the effect of the beam commutation, the same setup without beams would indeed give numerically zero residuals. Two effects due to the beams contribute to the residuals: first, the fact that when smoothing with  \nB final , because of the non-linearity of the spectral parameters and their spatial variability, the scaling laws are not correct anymore; second the fact that B final and A do not commute. However we do not have a way to distinguish between these two contributions in the residuals. We find that these latter are comparable with the target signal in the fullsky case, while they decrease in case we mask the Galactic plane, Figure 8.  \nTo account for the spatial variability of the mixing matrix in the following we consider a simplified version of what was done in [9], where they allow for the spatial variability of the mixing matrix for each spectral parameter on different sky pixels, corresponding to different values of nside . As our goal here is not to perform the best component separation possible for the given experimental setup, but rather study the effect of the beams, we do not explore all possible options but, instead, choose one, which is reasonable given what is shown in [9]. We thus consider the spatial variability values of nside =64 for β d , nside =8 for T d , and nside =4 for β s .  \nWe perform here only the second step of the component separation with the beams, assuming the spectral parameter values are given. First we solve the PCG in the noiseless case, with the true spectral parameter values in one patch per pixel, to verify if we are able to recover the components as expected. For the data model Eq. (9) we can recover the correct components with numerically zero residuals, this is true for any FWHM value of B final , even if we neglect it altogether as in the data model in Eq. (4), and we do it assuming a /lscript max value of about 300 as this is where the residuals will start to increase significantly due to the convolved noise. In contrast, this is not the case for the data model of Eq. (11), see Figure 9, in which the residuals due to the commutation, for different fraction of the sky are shown. The residuals that we get in this latter case give us the approximate measure of the error introduced assuming the commutation of  \nFIG. 2: Average of the recovered foreground components and CMB residuals for 1000 noise realizations, for the LiteBIRD -like setup with input full sky c1d0s0 . The components are shown at the frequency of 100 GHz. With the proposed approach dealing with the beams in the component separation we can access smaller angular scales as compared to the approach dealing with the beams before the component separation.  \n<!-- image -->  \nBB  \n<!-- image -->  \nCMB residuals (common reso approach)  \ndust (common reso approach)  \nsynch (common reso approach)  \nCMB residuals (beams in comp sep)  \ndust (beams in comp sep)  \nsynch (beams in comp sep)  \nD reference (r=0.001)  \nFIG. 3: Noise spectra after component separation for a given simulation, for the two approaches component separation approaches.  \n<!-- image -->  \nA and B final with this particular experimental setup. We note however that the level of the residuals is much smaller than the targeted signal, which means that with this particular instrumental setup and amount of spatial variability performing the commutation is an acceptable approximation.  \nThe value of B final has been chosen to be 17 ' , hence as the smallest of the true beams. In fact we find that the choice of the final beam size is largely irrelevant in the cases studied in this work, to the extent that we can neglect it altogether. Instead, a care has to taken about the band-limiting of all the pixel domain objects. This may or may not be the case for more complex simulations and real data, where subpixel power aliasing may become important. We therefore believe that using the explicit final beam is always prudent. We also find that it has an effect on the performances of the PCG.  \nWe have studied how the number of required PCG steps varied to reach the tolerance of 10 -7 as a function of the B final chosen. The qualitative behaviour that we have found does not change with the chosen nside or beam value: we see that the minimum of iterations is needed for a beam width in between the values of the true beams, closer to the lower bound. For smaller values the number of iterations needed significantly increases, the extreme case is the one of no B final , equivalent to the data model of Eq. (4). While on the other hand for larger beams the increase in number of iterations is even steeper. We explain this behaviour as a trade-off between having badly constrained modes (too small  \nCMB residuals Q  \nFIG. 4: LiteBIRD -like experiment, input sky c1d0s0 full sky, component separation with common resolution approach, first row, versus accounting for the beams in the component separation, second row. We show a patch of the sky for the recovered CMB residuals, thermal dust and synchrotron maps of about 30 · × 30 · centered on the galactic plane. We plot the Q Stokes parameter. The improvement in resolution in the proposed approach is appreciable even by eye and is more important for the dust component.  \n<!-- image -->  \n<!-- image -->  \nDust Q  \n<!-- image -->", 'Synchrotron Q': "beam size) or strongly correlated modes (too large beams) both negatively affecting the convergence. Hence, having both B true and B final with similar values allows to have better numerical stability. We show an example of this behaviour for a case at nside =64 and true beams ranging from 150 to 200 ' . We also notice that, given a particular B final , the data model with commutation, Eq. (11), converges faster than the one without commutation, Eq. (9), requiring about a factor 5 less iterations.  \nWe now perform the second step of the component separation assuming the spectral parameter values as recovered from a run with the common resolution approach . We compare the recovered components with the two data models (9) and (11), see Figure 11. We can here make similar considerations to what is discussed in the case without spatial variability. In particular with the data models with the beams we can reach smaller scales before the recovered residuals increase because dominated by the noise convolved to a combination of the B tue values, with respect to the common resolution approach. Again, a much wider range of /lscript is now gained for the thermal dust component, while less for the synchrotron components, according to the frequency bands where the gain on the beams resolution is more important. We also notice that the residual is now not fully dominated by the noise, as it was in the case without spatial variability, but are instead mainly dominated by systematic and statistical errors on spectral parameter estimates. The noise starts dominating only at smaller angular scale, where the residuals  \ngo up. Those scales, from about /lscript > 100, as are the ones affected by the beams are the ones which requires more PCG iterations to converge. Indeed, on large scales, thanks to the preconditioner used, good convergence is reached after few iterations ( O (1)), while all the following iterations are needed to bring the final scales up as due to the beams. The number of iterations needed strongly depends on the values used in the mixing matrix. When using the true spectral parameter values used also to create the simulations in this experimental setup with nside =512, about 40 iterations are needed, while this more and more iterations are required as the spectral parameters differ from these values. In particular with the spectral parameter used here form the common resolution run in the multiresolution configuration described above, with relative small nside the number of iterations required is of order of tens, while it significantly increases with nside =512, up to order of thousands to reach a PCG tolerance of 10 -6 . To reduce the number of iteration, making the approach more feasible we reduce the scale to which we recover the components. This is because the smaller angular scales are the ones mostly affected by the beam and thus requiring more iterations to converge. Typically it is enough to keep the scales up to where the noise starts dominating and exploding. Limiting the /lscript range we consider also has the benefit of reducing the computational cost of the SHTs.  \nOn a practical note, to allow for the PCG convergence one should ensure that the RHS and LHS of Eq. (16) have the same band limits. This is particularly the case with  \n<!-- image -->  \nDust Q  \n<!-- image -->  \nFIG. 5: LiteBIRD -like experiment, input sky c1d0s0 full sky, component separation with common resolution approach, first row, versus accounting for the beams in the component separation, second row. We show a patch of the sky for the recovered CMB residuals, thermal dust and synchrotron maps of about 30 · × 30 · centered in a clean region of the sky at longitude=60 · and latitude=40 · . We plot the Q Stokes parameter. The improvement in resolution in the proposed approach is appreciable even by eye and is more important for the dust component.  \n<!-- image -->  \nFIG. 6: Commutation of A and B final in the full-sky, true spectral parameter case, as it would occur in the common resolution approach. We use maps at nside =512 and FWHM of B final = 70 . 5 ' . The commutation introduces a non-negligible error except for the band at which the mixing matrix is normalized and thus spatially constant. The lensing B-modes curve is plotted for reference.  \n<!-- image -->  \nFIG. 7: Commutation of A and B final for the 235 GHz band. The plot shows ABs -BAs with and without masking the Galactic plane. The mask used is the f sky = 60% mask of [1]. Note that the values are driven by the fact that the mixing matrix is normalized at the frequency of 100 GHz. Bands closer to 100 GHz will thus have smaller amplitudes of these residuals, while it will be larger for bands further from 100 GHz.  \n<!-- image -->  \nSynchrotron Q  \nFIG. 8: Residual due to the beam in the common resolution approach , arising both from the commutation of the mixing matrix with the beam operator and the smoothing of the original maps with the beam. The residual is comparable with the targeted signal for the full-sky case, and is progressively reduced by masking the Galactic plane where the spectral parameters present the most spatial variability.  \n<!-- image -->  \nFIG. 9: CMB residuals due to the commutation of the mixing matrix and the beam operator B final in the proposed approach, for different Galactic plane masks. The level of the residuals due to the beam commutation is here lower than how it is for the approach with the common resolution beam, thanks to the smaller beam that we can use in this approach.  \n<!-- image -->  \nFIG. 10: Number of PCG iterations needed to recover the components with different choices of B final . Example for maps at nside =64 and with B true FWHM in the range from 150 to 200 ' . The minimum number of iterations required is for a B final FWHM value inside the range of values for B true and closer to its lower bound.  \n<!-- image -->  \nthe generalized mixing matrix, Eq. (12), as it requires the application of the beam operator as the last operation, which imposes a band limit when going back from harmonic to pixel domain, while it may not be the case for the data model of Eq. (13), for which then an additional, explicit bandlimiting of the RHS and LHS after the application of A is required.", 'B. Cut-sky case': "As a study case here we consider a Simons Observatory -like setup, and in particular its small aperture telescopes, with the specifics given in Table II [2]. We perform the sec-  \nTABLE II: Frequency, sensitivity and FWHM of Gaussian beams for the SO -like setup [2].  \nond step of the component separation with the signal only coming from about 10% of the sky. Given the partial sky our method is now only approximately valid. In the following we want to assess the impact of this approximation.  \nWe perform the second step of the component separation procedure on d 1 s 1 noiseless simulations ( nside =512 and /lscript max for the SHTs in the PCG 1024), assuming their true spectral parameters. This choice allows us to assign the full residuals that we get to the error introduced by the SHTs with the partial sky, as on full sky we would get numerically  \n<!-- image -->  \nFIG. 11: Recovered components (100 GHz) for sky d 1 s 1, LiteBIRD -like, full-sky and partial-sky (Galactic plane masked with f sky = 0 . 4). The first step of the component separation is performed with the common resolution approach, while the second step with the beams in the component separation with both data model proposed, Eq. (9) and Eq. (11). In the two cases the recovered components are in very good agreement and are plotted here up to the scales at which the noise starts dominating. The smaller scales would be out of reach if performing the second step with the common resolution approach. The CMB residuals could be further lowered by optimizing both the patches to deal with the foreground scaling spatial variability and the Galactic plane mask.  \n<!-- image -->  \nzero residuals. For numerical efficiency we set the beamsize of the beams smaller than 17 ' to 17 ' .  \nThe recovered components are shown in Figure 12. The main bias is visible along the patch border due to the signal being zero inside the mask. We can easily remove those biased pixels at the cost of reducing the sky fraction. We also notice that for an apodized initial mask the effect of the border is faded. We show the recovered B -mode residual power spectrum for this latter case, Figure 13. Given the error due to the mask is limited and not impacting the scientific target the proposed method can be used in this study case too, with analogous gain in resolution to what describe in the previous section. At the cost of loosing numerical efficiency one could keep the lower resolution channels to their native values and performing the computation at nside =1024, which would lead to even further gain in resolution.", 'IV. Conclusion and prospects': 'In this work we have discussed a treatment of main beam in the context of the parametric component separation and the input single frequency maps characterized by different observational beams. We have considered the standard approach, referred to as a common resolution approach, where all the input maps are smoothed to the common resolution prior to the component separation. We have found out that this approach not only tends to produce results with lower resolution but also that it implicitly assumes commutations of the common beam operator and mixing matrix. We have shown that this is the case even for mildly spatially-varying foreground properties, and this unavoidably leads to a bias, which can, and has to, be controlled with appropriate masking of the high contrast regions, where the variability of the foreground properties is particularly important. While we have found the bias to be manageable in the examples studied here, it depends strongly on the scale of the smoothing and thus it increases for large beams, becoming progressively more important and/or requiring more aggressive masking. We have then proposed two alternative approaches where the beams are explicitly treated as part of the component separation procedure. They either do not assume the beam and mixing matrix commutation at all, Eq. (9), or assume but only for beams much finer that the ones typically required for the common resolution approach, Eq. (11). We have developed numerically viable implementations of both these methods, which are algorithmically much more involved than the common resolution approach. We have found out that indeed these new methods allow us to recover sky signals over a significantly broader range of angular scales than the common resolution approach and can do so while controlling much better the level of beam-induced residuals. We have demonstrated both these methods on nearly full-sky simulations as expected from the next generation satellite missions and on the limited sky observations as could be observed from the ground. We have found out that the boundary effects due to the cut sky can be successfully dealt with. From the two proposed methods we have found that the method which does not assume any beam-mixing matrix commutation tends indeed to produce the most accurate results, however, this comes at the cost of significant numerical load, as the PCG solver in this case converges very slowly. Instead, the simpler method, which assumes the commutation of a high resolution beam and the mixing matrix, seems to offer an interesting trade-off between the computational efficiency and the level of bias, ensuring that the latter is subdominant for the sensitivity levels expected for the next generation of experiments and that the solution is produced in reasonable compute time. In both these cases further algorithmic improvements are however still needed to allow for a full exploration of the relevant likelihoods. This is a target of the on-going work.  \nThe obtained results are very promising and provide a proof-of-concept for this type of approaches, which will be instrumental in facilitating the full exploitation of the scientific potential of the future CMB data sets. The current implementation is however limited to axially symmetric beams and either spatially correlated stationary noise or inhomogeneous, uncorrelated noise. All these limitations can, and should, be bypassed with new, more efficient numerical  \nFIG. 12: Recovered CMB residuals for an SO -like patch, d 1 s 1 noiseless sky and with true spectral parameter values, the recovered CMB residuals are thus due to the error introduced by performing the SHTs on the cut-sky. In the first column the majority of the residual is in the pixels along the patch border, those have been masked in the central column, while on the right column the PCG has been solved on the apodized patch (with effective f sky ≈ 0 . 09). Both these two latter cases shows residuals by about an order of weaker than the former case.  \n<!-- image -->  \nalgorithms and their implementations. Furthermore, a full treatment of the problem will also require an inclusion of instrumental systematic effects. We leave such extensions for future work noting only here that all these effects are also relevant for the common resolution approach, where their impact also has to be yet assessed and fully understood.', 'Acknowledgments': "The authors thank Dominic Beck and Hamza El Bouhargani for useful discussions in the initial phases of the project. AR also thanks Magdy Morshed for various feedback during the project. Some of the results in this paper have been derived using the healpy [27] and HEALPix package [15] and with the NaMaster package [3]. The authors acknowledge the SCIPOL project funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 101044073). This work was granted access to the HPC resources of IDRIS under the allocation 2024-AD010414161R1 made by GENCI. Numerical computations were partly performed on the DANTE platform, APC, France.  \nFIG. 13: BB power spectrum from the components reconstructed with the data model of Eq. (9) and the PCG with tolerance 10 -6 , with the apodized mask with effective f sky ≈ 0 . 09. The components are plotted for the frequency of 100 GHz. The input sky is d 1 s 1 noiseless and the spectral parameters used are the true ones, the recovered CMB residuals are thus due to the error introduced by performing the SHTs on the cut-sky. They are lower than the targeted signal and they could be further reduced by removing the pixels on the border of the patch.  \n<!-- image -->  \n- L. Colombo, H. Eriksen, et al. Beyondplanck-i. global bayesian analysis of the planck low frequency instrument data. Astronomy & Astrophysics , 675:A1, 2023.\n- [5] C. Baccigalupi, L. Bedini, C. Burigana, G. De Zotti, A. Farusi, D. Maino, M. Maris, F. Perrotta, E. Salerno, L. Toffolatti, et al. Neural networks and the separation of cosmic microwave background and astrophysical signals in sky maps. Monthly Notices of the Royal Astronomical Society , 318(3):769-780, 2000.\n- [6] W. Brandt, C. Lawrence, A. Readhead, J. Pakianathan, and T. Fiola. Separation of foreground radiation from cosmic microwave background anisotropy using multifrequency measurements. Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 424, no. 1, p. 1-21 , 424:1-21, 1994.\n- [7] J.-F. Cardoso, M. Le Jeune, J. Delabrouille, M. Betoule, and G. Patanchon. Component separation with flexible models-application to multichannel astrophysical observations. IEEE Journal of Selected Topics in Signal Processing , 2(5):735-746, 2008.\n- [8] J. Chluba, J. C. Hill, and M. H. Abitbol. Rethinking cmb foregrounds: systematic extension of foreground parametrizations. Monthly Notices of the Royal Astronomical Society , 472(1):1195-1213, 2017.\n- [9] L. Collaboration, E. Allys, K. Arnold, J. Aumont, R. Aurlien, S. Azzoni, C. Baccigalupi, A. Banday, R. Banerji, R. Barreiro, et al. Probing cosmic inflation with the litebird cosmic microwave background polarization survey. Progress of Theoretical and Experimental Physics , 2023(4):042F01, 2023.\n- [10] E. De la Hoz, P. Vielva, R. B. Barreiro, and E. MartínezGonzález. On the detection of cmb b-modes from ground at low frequency. Journal of Cosmology and Astroparticle Physics , 2020(06):006, 2020.\n- [11] H. El Bouhargani, A. Jamal, D. Beck, J. Errard, L. Grigori, and R. Stompor. Mappraiser: A massively parallel mapmaking framework for multi-kilo pixel cmb experiments. Astronomy and Computing , 39:100576, 2022.  \n- [12] H. Eriksen, J. Jewell, C. Dickinson, A. Banday, K. Górski, and C. Lawrence. Joint bayesian component separation and cmb power spectrum estimation. The Astrophysical Journal , 676(1):10, 2008.\n- [13] H. K. Eriksen, A. Banday, K. Górski, and P. Lilje. On foreground removal from the wilkinson microwave anisotropy probe data by an internal linear combination method: limitations and implications. The Astrophysical Journal , 612(2):633, 2004.\n- [14] H. K. Eriksen, C. Dickinson, C. Lawrence, C. Baccigalupi, A. Banday, K. Górski, F. Hansen, P. Lilje, E. Pierpaoli, M. Seiffert, et al. Cosmic microwave background component separation by parameter estimation. The Astrophysical Journal , 641(2):665, 2006.\n- [15] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. ApJ, 622:759771, Apr. 2005.\n- [16] https://fgbuster.github.io/fgbuster/index.html.\n- [17] https://github.com/a-rizzieri/fgbuster-beams.\n- [18] B. Jost, J. Errard, and R. Stompor. Characterizing cosmic birefringence in the presence of galactic foregrounds and instrumental systematic effects. Physical Review D , 108(8):082005, 2023.\n- [19] N. Katayama and E. Komatsu. Simple foreground cleaning algorithm for detecting primordial b-mode polarization of the cosmic microwave background. The Astrophysical Journal , 737(2):78, 2011.\n- [20] R. Stompor, J. Errard, and D. Poletti. Forecasting performance of cmb experiments in the presence of complex foreground contaminations. Physical Review D , 94(8):083526, 2016.\n- [21] R. Stompor, S. Leach, F. Stivoli, and C. Baccigalupi. Maximum likelihood algorithm for parametric component separation in cosmic microwave background experiments. Monthly Notices of the Royal Astronomical Society , 392(1):216-232, 2009.\n- [22] M. Szydlarski, L. Grigori, and R. Stompor. Accelerating the cosmic microwave background map-making procedure through preconditioning. Astronomy & Astrophysics , 572:A39, 2014.\n- [23] B. Thorne, J. Dunkley, D. Alonso, and S. Naess. The python sky model: software for simulating the galactic microwave sky. Monthly Notices of the Royal Astronomical Society , 469(3):2821-2833, 2017.\n- [24] C. Vergès, J. Errard, and R. Stompor. Framework for analysis of next generation, polarized cmb data sets in the presence of galactic foregrounds and systematic effects. Physical Review D , 103(6):063507, 2021.\n- [25] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods , 17:261-272, 2020.\n- [26] K. Wolz, S. Azzoni, C. Hervías-Caimapo, J. Errard, N. Krachmalnicoff, D. Alonso, C. Baccigalupi, A. B. Lizancos, M. L. Brown, E. Calabrese, et al. The simons observatory: pipeline comparison and validation for large-scale bmodes. arXiv preprint arXiv:2302.04276 , 2023.\n- [27] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon, and K. Gorski. healpy: equal area pixelization  \npython. Journal of Open Source Software , 4(35):1298, Mar. 2019."}