A Partial-Sky Gibbs ILC Approach for the Estimation of CMB Posterior over Large Angular Scales of the Sky

In this article we present a formalism to incorporate the partial-sky maps to the Gibbs ILC algorithm to estimate the joint posterior density of the Cosmic Microwave Background (CMB) signal and the theoretical CMB angular power spectrum given the observed CMB maps. In order to generate the partial-sky maps we mask all the observed CMB maps provided by WMAP and Planck satellite mission using a Gaussian smoothed mask formed on the basis of thermal dust emissions in Planck 353 GHz map. The central galactic region from all the input maps is removed after the application of the smoothed mask. While implementing the Gibbs ILC method on the partial-sky maps, we convert the partial-sky cleaned angular power spectrum to the full-sky angular power spectrum using the mode-mode coupling matrix estimated from the smoothed mask. The main products of our analysis are partial-sky cleaned best-fit CMB map and an estimate of the underlying full-sky theoretical CMB angular power spectrum along with their error estimates. We validate the methodology by performing detailed Monte Carlo simulations after using realistic models of foregrounds and detector noise consistent with the WMAP and Planck frequency channels used in our analysis. We can estimate the posterior density and full-sky theoretical CMB angular power spectrum, without any need to explicitly model the foreground components, from partial-sky maps using our method. Another important feature of this method is that the power spectrum results along with the error estimates can be directly used for cosmological parameter estimations.


INTRODUCTION
Accurate measurement of the fluctuations present in the temperature and polarization fields of CMB anisotropies provides us with a wealth of knowledge regarding the origin, geometry and composition of our Universe (e.g., see Aghanim et al. (2020c)). However, strong emissions in the microwave region of the frequency spectrum by various astrophysical sources, the foregrounds, hinders us from directly observing the underlying true CMB signal. Therefore, in order to disinter the physics encoded in the CMB anisotropies the foremost challenge is to properly remove the foregrounds present in the CMB observations and thereby to recover the underlying CMB signal.
One of the important methods to minimize the foregrounds in a (foreground) model-independent manner is the so-called the internal-linear-combination (ILC) method (Tegmark & Efstathiou 1996;Tegmark et al. 2003;Bennett et al. 2003;Eriksen et al. 2004;Hinshaw et al. 2007;Saha et al. 2006;Saha et al. 2008;Saha 2011;Saha & Aluri 2016;Sudevan et al. 2017). A unique feature of the ILC method is that it does not require an explicit modeling of foregrounds in the form of templates or the spectral index information in order to remove them. This method purely relies on the assumption that the CMB photons follow a blackbody distribution along every direction of the sky so that its temperature fluctuations in thermodynamic units is independent of frequency. In the ILC method, foreground minimization is achieved by linearly superposing all the available foreground contaminated observed CMB maps with some amplitude terms, the weight factors. These weights can be computed analytically by minimizing the variance of the cleaned map. Instead of using the usual unweighted variance, some of the authors devel-oped a new method, the global ILC method (Sudevan & Saha 2018), where the weights are estimated by minimizing a theoretical CMB covariance weighted variance. This method outperforms the usual ILC over large angular scales of the sky. Unlike the usual ILC, the new global ILC weights could effectively minimize the contributions from chance-correlations between CMB and foregrounds over large angular scales of the sky. Some of the authors, Sudevan & Saha (2020a,b), developed a Gibbs ILC method to estimate the CMB posterior density and corresponding theoretical CMB angular power spectrum given the observed data over the large angular scales of the sky in a (foreground) model-independent manner. The Gibbs ILC method involves two important steps: • Given the observed CMB data and thoretical CMB angular power spectrum, a cleaned CMB map is sampled by minimizing the foregrounds present in the observed CMB maps using the global ILC method.
• Sample a new theoretical CMB angular power spectrum from its conditional density given the cleaned CMB map and observed CMB data.
To implement Gibbs ILC method in pixel space is a computationally exhaustive task since it involves estimation of cleaned CMB maps using global ILC method (as discussed in Sudevan & Saha (2018)) during each Gibbs iteration. In Sudevan & Saha (2020a), some of the authors implemented global ILC method in spherical harmonic basis in order to preform foreground removal. Further, the efficiency of the foreground removal using the global ILC method improves when one carries out the foreground minimization in an iterative manner (Sudevan & Saha 2018). Now taking into account that the Gibbs ILC method discussed in Sudevan & Saha (2020a,b); Purkayastha et al. (2020); Yadav & Saha (2020) uses full-sky CMB observations, it is, therefore, only natural to ask whether the Gibbs ILC method can be generalized such that it can handle maps with incomplete-sky data?
The challenges of dealing with incomplete-sky maps are the following: • Given partial-sky observed CMB maps and partial-sky cleaned CMB map, how to sample a theoretical CMB angular power spectrum?
• How to incorporate a full-sky CMB angular power spectrum into the global ILC algorithm during a partialsky analysis?
• Finally, how to address the issue of "Gibbs phenomenon" 3 in a spherical harmonic analysis involving discrete masks (i.e., mask with pixel values 1 and 0 only)?.
In order to develop a truly partial-sky Gibbs ILC method these three questions need to be answered.
In this article, we incorporate the MASTER approach (Hivon et al. 2002) to our Gibbs ILC algorithm to sample a partial-sky CMB angular spectrum given the fullsky spectrum. As for Gibbs phenomenon, there exist atleast two possible solutions for suppressing its adverse affects. One of them is to allow the spherical harmonic expansion to infinitely high multipoles (Bankman 2000;Allen & Mills 2004;Carslaw 1921). Another option is to use a mask smoothed by a Gaussian beam to minimize the jump discontinuities at the region boundaries. Since there exists a maximum multipoles max upto which spherical harmonic coefficients can be expanded for a given (finite) pixel resolution of the HEALPix 4 map, utilizing a smoothed mask is, therefore, a more viable option. Hence a work-around to suppress the Gibbs phenomenon effectively in our Gibbs ILC method is to use fullsky observed CMB data multiplied with a smoothed mask. We organize the paper as follows. We discuss in Section 2 the basic formalism of our approach. We describe the input maps, masks that we use in this analysis in Section 3. In Section 4 we lay out our procedure and discuss how to estimate the posterior density in a partial-sky analysis. We present our results of analysis of WMAP and Planck frequency maps at low resolution in Section 5. In Section 6 we validate our method for estimating the joint posterior density of CMB signal and its theoretical angular power spectrum by performing detailed Monte Carlo simulations. Finally, we conclude in Section 7.

The Gibbs ILC Method
In the Gibbs ILC method, to estimate the joint CMB posterior density, P(S,C |D), we draw samples of S and C from their respective conditional densities using the Gibbs sampling (Eriksen et al. 2008(Eriksen et al. , 2007Geman & Geman 1984) technique. Here S, C and D represent the true CMB signal, theoretical CMB angular power spectrum, and the observed CMB data respectively. At the beginning of some iteration (i + 1) in a Gibbs sampling procedure, a cleaned CMB signal S i+1 and a theoretical CMB angular power spectrum C i+1 is sampled from their respective conditional densities P 1 (S|D,C ) 3 Refers to the peculiar manner in which the Fourier representation behaves at a jump discontinuity Gibbs (1898); Bankman (2000); Allen & Mills (2004); Carslaw (1921). 4 Hierarchical Equal Area Iso-Latitude Pixelation of Sphere which was developed by Górski et al. (2005). and P 2 (C |D, S) as follows, (1) Using the pair of samples, S i+1 and C i+1 generated at the end of (i + 1) th iteration, we repeat the above two sampling steps (Eqns. 1 and 2) for a large number of iterations till convergence is achieved. Ignoring few samples generated during the initial (burn-in) phase, rest of the samples of S and C appear as if they are sampled from the joint posterior density P(S,C |D).
In the Gibbs ILC method, in order to estimate P(S,C |D) in a (foreground) model-independent manner, we minimize the foregrounds present in the observed CMB maps using the global ILC algorithm at each Gibbs step. Using the global ILC method, with CMB maps d i observed at n different frequencies, an estimateŜ of the underlying true S is obtained asŜ where w i is the weight corresponding to the map from i th frequency channel. These weights are subject to a constraint that they should sum to unity, i.e., n i=1 w i = 1. This constrain ensures that the CMB signal will not undergo any undesired modification during the entire foreground removal procedure. By maximizing the likelihood of the model given the full-sky observed CMB data and the theoretical CMB angular power spectrum in the spherical-harmonic domain, the conditional density P 1 (S|D,C ) is obtained as follows: σ i j in the above equation is the cross power spectrum between the observed CMB maps d i and d j . Form the Eqn. 5, we define an estimator σ 2 r given by which is then minimized in order to obtain the weights. C (= C B 2 P 2 ) in Eqn. 6 is the beam and pixel smoothed theoretical CMB angular power spectrum (C ). The choice of weights which minimizes σ 2 r is obtained by following a Lagrange's multiplier approach, where,Â W is a (n × 1) weight vector and e is a (n × 1) shape vector of CMB in thermodynamic units.
To draw the samples of C given S and D in the Gibbs ILC method, we first define a variable z =Ĉ (2 + 1)/C , where C is estimated from the cleaned CMB map. We obtain the conditional density P 2 (C |D, S) as From Eqn. 9, one can infer that the variable z follows a χ 2 distribution with 2 − 1 degrees of freedom. We then sample a theoretical CMB angular power spectrum as follows where z is drawn from the χ 2 distribution of 2 − 1 degrees of freedom.
2.2. Partial-Sky Analysis The Gibbs ILC method (Sudevan & Saha 2020a,b), discussed in Section 2.1, uses full-sky observed CMB maps provided by WMAP (Jarosik et al. 2011) and Planck (Aghanim et al. 2020a) CMB observations. Oftentimes in the point of view of the foreground removal problem it is advantageous to mask certain highly contaminated regions of the sky like the Galactic plane, strong point sources etc. These maskedobserved CMB maps are then fed into the foreground minimization pipeline. The advantage of implementing this approach is that the ILC (and the global ILC) weights can be tuned in such a way that they are estimated based on the strength of the foreground contaminations across the sky. A similar partial-sky approach is not direct in the case of Gibbs ILC method. Even though one can sample a partial-sky CMB map by removing the foregrounds using global ILC method given the partial-sky observed maps in pixel space. However, a similar implementation of the global ILC method in spherical harmonic space for the partial-sky maps is not straightforward. Another challenge is to understand how to sample a theoretical CMB angular power spectrum given a partial-sky cleaned CMB map and partial-sky CMB data.
In the current analysis, we formalize a technique such that the Gibbs ILC method can incorporate only certain regions of the sky (preferably low foreground contaminated regions) from the full-sky observations to estimate the joint CMB posterior density and full-sky cleaned theoretical CMB angular power spectrum.

Sampling a Theoretical CMB angular power spectrum
In this section, we briefly review the MASTER algorithm (Hivon et al. 2002) implemented in our framework to sample a theoretical CMB angular power spectrum from its conditional density given a partial-sky cleaned CMB map and partial-sky CMB data. We begin with a full-sky CMB anisotropic map which can be expanded in terms of spherical harmonics as follows with a m = dΩ S(n)Y * m (n) and angular power spectrum defined as, A partial-sky map (S) used in this analysis is a full-sky map multiplied with a Gaussian smoothed mask (M) i.e., Details about the smoothed mask is given in the Section 3.2. Substituting Eqn. 11 in Eqn. 13, we can expand the masked map in spherical harmonic representation as follows where K m m is the coupling kernel defined as The mask can be expanded in spherical harmonics with the coefficients Now substituting Eqn. 17 in Eqn. 16, the coupling kernel reads as m m m is Wigner 3-j symbol also known as Clebsch-Gordan coefficient which describes the coupling of the three angular momentum vectors. After using the orthogonality condition of Wigner Symbol, where, substitute the kernel expansion Eqn. 18 and Eqn. 15 in Eqn. 12 we getC where, C M is the smoothed mask (M) power spectrum. Defining the mode-mode coupling matrix corresponding to the mask as, FIG. 1.-We present various masks that is used in this analysis. Left panel shows the mask based on the thermal dust emission, the ThDust5000 mask. The central galactic region of the ThDust5000 mask (shown in violet) corresponds to very intense thermal dust emissions. In the middle panel we show the mask obtained after smoothing the ThDust5000 mask, the ThDust5000_Sm, using a Gaussian beam of FWHM = 9 • . In the last panel we show the effective mask, i.e., ThDust5000_Eff, which consists of all pixels set to 1 whose corresponding value in ThDust5000_Sm mask is ≥ 0.95 and rest 0.
we finally arrive atC = M C .
The importance of using Eqn. 26 is that we can estimate angular power spectrum corresponding to the unmasked region by multiplying the full-sky power spectrum with the mode-mode coupling matrix. Similarly, the converse is also true i.e., We can use Eqn. 27 to estimate the full-sky angular power spectrum from the partial-sky angular power spectrum given the mode-mode coupling matrix. During the implementation of the Gibbs ILC method, at any given Gibbs step i, we use Eqn. 27 to estimate the full-sky CMB angular power spectrum from the partial-sky CMB map. A new full-sky theoretical CMB angular power spectrum is then sampled from its conditional density as discussed in Section 2.1 using the full-sky CMB angular power spectrum. Afterwards, we employ Eqn. 26 to estimate the partial-sky CMB angular power spectrum from the full-sky sampled theoretical CMB angular power spectrum. This partial-sky CMB angular power spectrum is then used to sample a partial-sky cleaned CMB map. The details of sampling a partial-sky CMB map is presented in the following sub-section.

Sampling a CMB map in a Partial-sky Analysis
In Sudevan & Saha (2018), we discussed how to implement global ILC method in an iterative manner in pixel space. Given that during a partial-sky analysis, the observed CMB maps at any given N side consist ofÑ pix (where,Ñ pix < N pix 5 ) number of pixels, the theoretical CMB covariance weighted variance in this case in pixel space is as follows: HereS is the a partial-sky cleaned CMB map withÑ pix of pixels. Similarly,C is a (Ñ pix ×Ñ pix ) matrix consisting of theoretical CMB covariance between only the surviving pixels. But in the Gibbs ILC algorithm using global ILC in pixel space is a computationally expensive procedure since at each Gibbs step we need to evaluate the theoretical CMB covariance matrix using the sampled theoretical CMB angular power spectrum. Hence, instead of using pixel space global ILC we implement the global ILC method in harmonic space (Eqn. 8) in the full-sky observations. For a partial-sky Gibbs 5 where Npix is the total number of pixels in a HEALPIX full-sky map at a given pixel resoultion Nside. The Npix and Nside is related as follows: Npix = 12 × N 2 side ILC pipeline, we propose the following estimator based on the Eqn. 5 for estimating the weights: whereσ i j is the partial-sky cross-power spectrum corresponding to the (i, j) th input CMB maps. The partial-sky weights are calculated by minimizing Eqn. 29. We use these new weights to estimate a partial-sky cleaned CMB map as follows, 3. INPUT DATA 3.1. Input Maps In our analysis, we use all of the WMAP nine-year difference assembly (DA) maps (Jarosik et al. 2011) and seven Planck-2015 maps -three of them correspond to LFI frequencies (30, 40 and 70 GHz) (Ade et al. 2016) and the rest are four HFI frequencies (100, 143, 217 and 353 GHz) (Adam et al. 2016). The processing of all the input maps remains identical to Sudevan & Saha (2018) and results in a total of 12 input maps (five at WMAP and seven at Planck frequencies). All our input maps are at a pixel resolution N side = 16 and beam smoothed by a Gaussian beam of FWHM = 9 • after properly taking taking care of the native beam resolutions. At this pixel and beam resolution we can ignore the contributions of noise. The monopole is removed from all the full-sky input maps before the analysis.

Masks
In the current analysis we use three different types of mask. The details of the various masks are given below.

ThDust5000
In order to generate a mask based on the thermal dust emissions we follow the below steps.
1. We downgrade Planck 70 and 353 GHz maps to N side = 256 and both at same effective beam resolution of FWHM = 6 • .
2. We take the difference between Planck 353 and 70 GHz maps, which will give a map dominated by the thermal dust emissions.  3.-The Best-Fit partial-sky cleaned CMB map obtained following our Gibbs ILC procedure at N side = 16 and smoothed by a a Gaussian beam of FWHM = 9 • is shown in the top panel. While estimating the final best-fit map we masked the galactic region using our ThDust5000_Eff mask. In the middle left and middle right panel we show the difference maps comparing our best-fit map with Planck COMMANDER and NILC cleaned maps respectively. In the bottom panel from left to right we compare our best-fit map with SMICA and SEVEM cleaned maps respectively. All the difference maps are masked using ThDust5000_Eff mask. All the color scales are in µK thermodynamic temperature unit.
3. We set all the pixels in the difference map with pixel value ≤ 5000 µK to 1. All the remaining pixels are set to 0.
We refer to the mask obtained by following the above procedure as "ThDust5000" and is shown in the left panel of Fig. 1. Out of a total of 3072 pixels, the ThDust5000 mask consists of 679 pixels with pixel value '0' and rest of the pixels are '1'.

ThDust5000_Sm & and ThDust5000_Eff
As mentioned in the Introduction, using a binary mask such as ThDust5000 in the Gibbs algorithm results in Gibbs phenomenon at jump discontinuities across the mask boundaries. In order to avoid this situation, one feasible solution is to smooth the mask. In this work, we use the Gaussian smooth-ing kernel for smoothing the mask where θ is a separation angle. The Gaussian smoothing kernel can be considered as a low-pass filter with window function Here, σ = FWHM/ √ 8ln(2) and FWHM of the smoothing kernel is set to 9 • so that in our Gibbs ILC method, where we perform a spherical harmonic expansion up to = 32, all the the multipoles > 32 are sufficiently suppressed. The smoothed mask ThDust5000_Sm obtained after smoothing the ThDust5000 mask is shown in the middle panel of smoothed mask might take small positive values and similarly the pixels with initial pixel value 1 might no longer be 1 after the smoothing operation. As we move away from the region boundaries the pixel values will no longer be affected significantly by the smoothing process.
While estimating the final cleaned map we need to consider a new mask since some of the '1' valued pixels in ThDust5000 are no longer '1' in ThDust5000_Sm mask. Hence CMB estimated in those pixels will have strength lower than the actual value. In order to avoid the under estimation of CMB, we construct a new mask "ThDust5000_Eff" by setting all the pixels in ThDust5000_Sm mask whose pixel value is ≥ 0.95 as 1 and rest as 0. The ThDust5000_Eff consists of 1003 pixels with pixel value 0 out of total 3072 pixels.
Unlike the present case, where full-sky CMB maps are available from WMAP and Planck satellite missions, if we consider a CMB mission (e.g., a ground based observatory) where only part of the entire sky can be observed, then the smoothing of the mask should be modified as follows.
1. Let Mask1 selects the portion of the sky observed by the experiment. Use Mask1 to construct a new mask Mask2 such that the new mask lies well inside of Mask1.
2. The area of the new mask, Mask2, should be selected in such a way that after Gaussian smoothing the smoothed Mask2 pixel values in the unobserved region is = 0.

METHODOLOGY
We multiply all the input maps at a pixel resolution N side = 16 and smoothed by a Gaussian beam of FWHM = 9 • with the ThDust5000_Sm mask (discussed in Section 3.2) thus mimicking partial-sky observed maps. Our Gibbs ILC method has ten different chains, each with 10000 Gibbs steps. The initial choice of C for each chain is made by drawing them uniformly within ±3∆C around the Planck best-fit theoretical power spectrum (Aghanim et al. 2020b), where ∆C denotes error due to cosmic variance alone. Since our input maps are all partial-sky maps, we cannot use the above sampled C directly for computing the weights. For each chain, we multiply the respective initial C with the modemode coupling (M ) matrix to obtain the corresponding partial-sky C . The M is estimated corresponding to the ThDust5000_Sm mask, which relates a full-sky CMB angular power spectrum of a given map to the corresponding partial-sky angular power spectrum estimated from the ThDust5000_Sm masked map. The partial-sky angular power spectrum is then used to sampleS using Eqn. 30 where the weights are estimated by minimizing Eqn. 29. After samplingS we estimate the partial-sky cleaned CMB power spectrumC . In order to sample a full-sky theoretical CMB angular power spectrum, we use Eqn. 27 to estimate the full-sky angular power spectrumĈ usingC . This full-sky power spectrum is then mulitplied with ( + 1)/2π, beam and pixel window functions and later binned with a bin width of 3 mulitpoles. The bin-middle values are at = 3, 6, 9, 12, 15, 18, 21, 24, 27 and 30. After binning we divide the binned full-sky power spectrum by ( + 1)/2π, beam and pixel window functions. Now using this binned power spectrum we sample a new theoretical CMB angular power spectrum by following Eqn. 10. The new sampled full-sky theoretical CMB angular power spectrum is then multiplied by the M to obtain a sample of partial-sky power spectrum which is then used in the next Gibbs step.
At the end of all 10 chains, we are left with a total of 100000 joint samples of cleaned map and theoretical CMB power spectrum to sample the posterior density P(S,C |D). The burn in phase in each chain is very brief. Visually, this phase does not appear to contain more than a few Gibbs iterations. We, however, remove the initial 50 Gibbs iterations from each chain as a conservative estimate of burn-in period. After the burn-in rejection, we end up with a total of 99500 samples from all chains.

Cleaned Maps
In order to generate the best-fit cleaned CMB map we follow the procedure as outlined in Sudevan & Saha (2020a). We use all samples of cleaned CMB maps after rejecting those generated during the initial burn-in phase from each chain. We estimate the marginalized probability density functions of CMB temperature corresponding to the pixels not masked out by the ThDust5000_Eff mask. We show the normalized probability density functions corresponding to some random pixels in Fig. 2. The normalization for each probability density function is such that the peak corresponds to a value of unity. Taking the mode-value of the density function corresponding to each surviving pixel we form the best-  fit partial-sky cleaned CMB map which is shown in the top panel of Fig. 3. We compare our best-fit cleaned CMB map with Planck Commander, NILC (needlet space ILC), SMICA and SEVEM cleaned maps (Akrami et al. 2020). We downgrade the Planck cleaned maps to N side = 16 and smooth with a Gaussian beam of FWHM = 9 • after properly taking care of their respective individual beams. We show the difference maps obtained after taking the difference between our best-fit cleaned map and the cleaned maps provided by Planck COM-MANDER, NILC, SMICA and SEVEM in the middle and bottom panels of the Fig. 3. We remove the central galactic region by multiplying the maps with the ThDust5000_Eff mask. Clearly, from the difference map plots we can see that our best-fit map matches with rest of the cleaned CMB maps quite well.
Using all the 99500 partial-sky cleaned CMB maps from our Gibbs ILC method, we estimate a mean cleaned CMB map, the mean map. We show the mean map after multiplying with the ThDust5000_Eff mask in the top panel of Fig. 4. In the middle panel of Fig. 4 we show the difference between best-fit CMB map and mean CMB map. We see that the mean map matches very well with the best-fit CMB map within an absolute difference of < 1µK. In order to quantify the reconstruction error while estimating a partial-sky cleaned CMB using Gibbs ILC method given partial-sky input CMB data and theoretical CMB angular power spectrum, we generate a standard deviation map using all 99500 partial-sky cleaned maps. We show this map in the bottom panel of Fig. 4. We can infer from this figure that the reconstruction error is very small all over the region of the sky selected by ThDust5000_Eff mask.

Angular Power Spectrum
We estimate the marginalized probability density of CMB theoretical angular power spectrum for multipoles corresponding to the bin-middle values and show the respective density functions in Fig. 5. We normalize these density functions to a value of unity at their peaks. The horizontal axis of each plot of this figure represents ( + 1)C /(2π) in the unit of 1000 µK 2 . The density functions show long asymmetric tails for low multipoles (e.g., = 3, 6, 9). For large multipoles ( ≥ 21) the asymmetry of the densities become gradually reduced. The region within the two vertical lines in each plot corresponds to the 1σ (68.27%) confidence interval for the CMB theoretical power spectrum for each multipoles.
In the top panel of Fig. 6 we show the binned best-fit theoretical CMB angular power spectrum, with black colored points, defined by the positions of peaks of marginalized angular power spectrum density functions. The asymmetric error bars at each binned values show the 1σ confidence interval for the theoretical angular power spectrum. The light green line shows the theoretical power spectrum consistent with Planck 2018 results (Aghanim et al. 2020b). The best-fit theoretical angular power spectrum agrees well with the binned spectra estimated from Commander, NILC and SMICA cleaned maps, which are shown by red, blue and violet points respectively. In the bottom panel of Fig. 6, we show the difference of the best fit and Commander NILC, SMICA and SEVEM angular power spectra. We emphasize here that the all the Planck cleaned maps angular power spec- best-fit cleaned CMB map obtained from the first Monte Carlo simulation, followed by the difference map obtained after taking the difference between the best-fit cleaned CMB map the input CMB map used in that particular simulation. In the top right panel we show mean of all 1000 such difference maps estimated from each Monte Carlo simulation. In the bottom right pane, we show the standard deviation map from the first simulation. The standard deviation map signifies the error while reconstructing a cleaned CMB map following our Gibbs ILC algorithm. In the bottom middle and bottom right panels we show the mean and standard deviation of all such 1000 standard deviation maps. All the maps that we show in this figure are multiplied by ThDust5000_Eff mask. All the color scales are in µK thermodynamic temperature unit. tra are estimated from the respective full-sky maps.

MONTE CARLO SIMULATIONS
For our detailed Monte Carlo simulations, we generate 1000 sets of input maps at the 12 different WMAP and Planck frequencies at a Gaussian beam resolution 9 • and N side = 16 following the procedure as described in Sudevan & Saha (2018) and Sudevan & Saha (2020a). Our foreground model consists of synchrotron, free-free and thermal dust emissions. We use Planck 2018 foreground templates for generating foreground maps at all WMAP and Planck frequencies. The random CMB realizations that we use in the input frequency maps are generated using the CMB theoretical angular power spectrum consistent with Planck 2018 results. Once we simulate the foreground contaminated CMB maps, we then proceed to mask all the 1000 sets of 12 input maps using the ThDust5000_Sm mask. The masked input maps are then used in our Gibbs ILC code which consists of a total of 10 Gibbs chains, each with 10000 Gibbs steps. We follow the same methodology as adopted in the real data.
As was the case with real data, the initial burn-in period ends quickly. We continue keeping a conservative limit of initial 50 samples corresponding to the burn-in phase which are then removed from each Gibbs chain. After initial burnin rejection, we are left with a total of 9950 joint samples of cleaned map and theoretical angular power spectrum from each Gibbs chain.
Using all cleaned map samples from all chains after burnin rejection we form marginalized density functions corresponding to each surviving pixels in a N side = 16 map. We show the normalized density functions corresponding to some random pixels from a representative simulation (here for instance, the Monte Carlo Simulation No. 1) in Fig.7. A CMB map formed from the pixel temperatures corresponding to the modes of these density functions define the bestfit partial-sky CMB cleaned map obtained from the simulation. We only consider those pixels which have pixel value 1 in the ThDust5000_Eff while estimating the best-fit cleaned map. We show our partial-sky best-fit cleaned map from the above simulation in the top-left panel of Fig. 8. The difference of the best-fit and input CMB realization after masking with the ThDust5000_Eff mask is shown in the top middle panel of the same figure. Clearly, the best-fit CMB map matches very well with the input CMB realization used in the simulation. The maximum difference (∼ 10 µK) between the two maps is observed along the edges of the ThDust5000_Eff mask. This shows that our method removes foreground reliably in the outer regions of the Galactic plane. The top right panel of Fig. 8 shows the mean of all the 1000 difference maps estimated by taking the mean of difference between best-fit and input CMB maps obtained from all Monte Carlo simulations. The bottom right panel of Fig. 8 shows the error map computed from the Gibbs samples. The error map is a standard devaition map computed using all the cleand CMB samples from Monte Carlo simulation No. 1. The maximum error of 1 µK is observed along the edges of the ThDust5000_Eff masked region. We show the mean and standard deviation of 1000 such error maps from each Monte Carlo simulations in the bottom middle and bottom left panels of Fig. 8. We see that the error while reconstruct-  ing a cleaned CMB pixel following our partial-sky Gibbs ILC method is very small. In summary, using the Monte-Carlo simulations, we see that our Gibbs ILC method reliably minimizes the foreground across all the 1000 simulations.
We show normalized density functions corresponding to the bin-middle multipoles from some representative simulation in Fig.9. The normalization is performed such that the peak value of the density functions are set to unity. The horizontal axis for each sub-plot represents ( + 1)C /(2π) in the unit 1000 µK 2 . The red line in each sub-plot corresponds to the mode-value of the corresponding probability density functions. The blue line represents the mean of all 1000 mode values corresponding to each density function and the black dashed line is the value of Planck-2018 theoretical CMB angular power spectrum. We show the mean of the best-fit estimate of binned underlying CMB theoretical angular power spectrum from all the 1000 Monte Carlo simulations in Fig. 10 in blue points. The multipole values corresponding to each point is the bin-middle value. The underlying theoretical angular power spectrum is shown in light green. The blue vertical lines shows the standard deviation corresponding to each multipoles.

CONCLUSIONS & DISCUSSIONS
In this article, we formalize a foreground model independent approach to estimate the posterior density of CMB map and the corresponding theoretical angular power spectrum following Gibbs-ILC method using sky region outside the galactic plane. We use all the 10 detector set maps provided by WMAP and 7 frequency maps from Planck mission 44,143,271,353 GHz). Before foreground removal all the input maps are preprocessed as discussed in Section 3 which results in a total of 12 input maps at Healpix pixel resolution N side = 16 and smoothed by a Gaussian beam of FWHM = 9 • .
We exclude the galactic region of the sky by multiplying the input maps with a smoothed mask. Using smoothed mask mitigates the adverse effects due to improper behavior of spherical harmonic transforms at jump discontinuities around the mask boundaries. There are two major outcomes of generalizing the full-sky Gibbs ILC method over the partial-sky. First, the weights that are the primary variables to remove the foregrounds, can now be more effectively tuned to minimize the foregrounds based upon the local information of foreground strength or depending upon the nature of foreground contaminations in input maps. Secondly, the method can now be applied on the CMB observations which can produce maps only over a fraction of region of the entire sky (e.g., as expected from ground based CMB observations). Finally, the final CMB products of this new work are independent of detailed foreground modelling.
We perform detailed Monte Carlo simulations to validate the Gibbs ILC method over partial-sky. We simulate 1000 noise and foreground contaminated CMB maps (based on the foreground model provided by Planck 2018 results) at the same WMAP and Planck frequencies which are included in the case of foreground removal from real data. The simulation results show that our method effectively and accurately removes the foregrounds present in the unmasked region. The mean of all binned best-fit theoretical angular power spectra from all simulations agree excellently with the input Planck-2018 theoretical CMB angular power spectrum used to generate the random realizations of CMB maps.
For the case of real data, using all the 99500 samples of partial-sky cleaned CMB maps and binned estimates of full-sky CMB angular power spectra we estimate the best-fit (partial-sky) cleaned map and best-fit theoretical CMB angular power spectrum respectively. We see a nice agreement between our results and those obtained from the Planck COM-MANDER, NILC, SMICA cleaned CMB maps. This work is based on observations obtained with Planck (http://www.esa.int/Planck) and WMAP (https://map.gsfc.nasa.gov/). Planck is an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We acknowledge the use of Planck Legacy Archive (PLA) and the Legacy Archive for Microwave Background Data Analysis (LAMBDA). LAMBDA is a part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is supported by the Astrophysics Science Division at the NASA Goddard Space Flight Center. We use publicly available HEALPix Górski et al. (2005) package (http://healpix.sourceforge.net) for the analysis of this work.