An Application of the Global ILC Algorithm over Large Angular Scales to Estimate the CMB Posterior Using Gibbs Sampling

Since the discovery of the cosmic microwave background (CMB; Penzias & Wilson 1965), rapid advancements in the field of its observation made it possible to map the primordial signal over the entire sky with increasingly higher resolution (Smoot et al. 1991; Bennett et al. 2013; Planck Collaboration et al. 2019b). Accurate measurement of the temperature (and polarization) anisotropy of the CMB, which arguably forms one of the cornerstones of the precision era of modern precision cosmology, provides us with a wealth of knowledge regarding the geometry, composition, and origin of the universe (e.g., see Planck Collaboration et al. 2019a and references therein). However, the observed CMB signal in the microwave region is strongly contaminated by foreground emissions that are due to different astrophysical sources within and outside our galaxy. Hence the challenge is to accurately recover the CMB signal for a cosmological analysis by minimizing contributions from various foregrounds emissions.

For a reliable estimation of cosmological parameters, a desirable property of a CMB reconstruction method is that it produces the best guesses for the signal and its angular power spectrum, along with their corresponding error (and bias, if any) estimates. Eriksen et al. (2004b, 2007, 2008a, 2008b) and Planck Collaboration et al. (2016a, 2016b) propose and implement a Gibbs sampling (Geman & Geman 1984) approach to jointly estimate the CMB map, its angular power spectrum, and all foreground components along with their error estimates using the Wilkinson Microwave Anisotropy Probe (WMAP; Hinshaw et al. 2013) and Planck (Planck Collaboration et al. 2019b) observations. Eriksen et al. (2006) and Gold et al. (2011) use a maximum likelihood approach to simultaneously reconstruct CMB and foreground components using prior information about the CMB as well as detector noise covariance matrices and foreground models. Although these methods are extremely useful for a simultaneous reconstruction of the CMB and all foreground components, an alternative approach for the CMB reconstruction alone is the so-called internal-linear-combination (ILC) method (Tegmark & Efstathiou 1996; Bennett et al. 2003; Tegmark et al. 2003; Eriksen et al. 2004a; Saha et al. 2006), which does not rely upon any explicit model of the foreground spectrum. In recent years, the method has been investigated extensively. Saha et al. (2006) use this method to estimate CMB cross-power spectra by removing detector noise bias using WMAP maps. These authors also report a possible negative bias at the low multipoles. Saha et al. (2008) perform for the first time a rigorous analytical study of negative bias at the low multipoles for a single-iteration ILC foreground removal procedure in harmonic space. Later, Sudevan et al. (2017) find and correct for a foreground leakage in the iterative ILC algorithm in harmonic space by applying their technique on high-resolution Planck and WMAP observations. Sudevan & Saha (2018) propose a global ILC method in pixel space by taking into account prior information of the CMB covariance matrix under the assumption that detector noise can be ignored over the large angular scales of the sky. The method considerably improves the usual ILC method at low resolution, where no prior information about the CMB covariance is used. In spite of these progresses, a joint analysis of the CMB signal and its angular power spectrum posterior density in a manner independet of foreground model has not yet been explored in the literature. The current article is aimed to provide a mechanism to solve just this problem. By estimating the posterior density of the CMB signal and the CMB theoretical angular power spectrum given the observed data over the large angular scales of the sky using a ILC method similar to that used by Sudevan & Saha (2018), we provide the best-fit estimates of the CMB map and of the theoretical angular power spectrum, along with their confidence interval regions. In the current article, we replace the CMB signal reconstruction technique by a faster harmonic-domain algorithm than the pixel-space algorithm of Sudevan & Saha (2018). We use the Gibbs sampling method (German & German 1984) to draw samples from the joint conditional density. There are two important advantages of our method. First, the theoretical power spectrum results can directly be integrated into the cosmological parameter estimation process. Second, the CMB posterior estimation can be achieved without any need to explicitly model the foreground components. The results therefore cannot be sensitive to foreground modeling uncertainties.

The early work of reconstructing the CMB component has been performed by Bennett et al. (1992) using a variant of the ILC algorithm, where prior information of a free–free spectral index is used. Bunn et al. (1994) and Bouchet et al. (1999) developed a Wiener filter approach. Basak & Delabrouille (2012, 2013) proposed an ILC algorithm in needlet space, which can take into account local variation of foreground spectral properties both in the pixel and needlet space. Saha & Aluri (2016) used the ILC method to jointly reconstruct the CMB Stokes Q polarization signal and other foreground components in presence of spatially varying spectral properties of polarized synchrotron emission using simulated observations of WMAP. In an interesting application of the ILC method, Saha (2011) and Purkayastha & Saha (2017) reconstructed CMB maps using the Gaussian nature of the CMB and the non-Gaussian nature of astrophysical foregrounds.

In Section 2 we discuss the basic formalism. We describe the posterior estimation method in Section 3. In Section 4 we present the results of our analysis of WMAP and Planck frequency maps at low resolution. We discuss convergence tests of the Gibbs chains in Section 5. We validate the posterior density estimation method by performing detailed Monte Carlo simulations using realistic foreground and detector noise models consistent with WMAP and Planck observations in Section 6. We discuss the impact of the uncertainty on the residual calibration of the input frequency maps on the final CMB products in Section 7. Finally, we conclude in Section 8.

2.1. Data Model

Let us assume that we have observations of foreground-contaminated CMB maps at n different frequencies. Without sacrificing any generality, we assume that each of these maps has the same beam (and pixel) resolutions.¹ The observed data set,  D, can be represented as ${\boldsymbol{D}}=\{{{\boldsymbol{X}}}_{1},{{\boldsymbol{X}}}_{2},...,{{\boldsymbol{X}}}_{n}\}$, where ${{\boldsymbol{X}}}_{i}$, i ∈ {1, 2, ..., n} is an N × 1 column vector denoting the input foreground-contaminated CMB map (in thermodynamic temperature units) at a frequency ν_i. N represents the number of pixels in each input frequency map, and  D is an N × n matrix. Assuming detector noise is negligible² , we have

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{i}={\boldsymbol{S}}+{{\boldsymbol{F}}}_{i},\end{eqnarray} \tag{ 1 }$$

where  S is an N × 1 column vector, representing the CMB signal³ , and  F_i denotes a map of the same size, representing the net foreground contamination at the frequency ν_i.

2.2. CMB Posterior and Gibbs Sampling

2.3. Density of the Pure CMB Signal

The probability density function of the pure CMB signal  S given a theoretical CMB angular power spectrum C_ℓ is given by

$$\begin{eqnarray}&&{P}_{3}({\boldsymbol{S}}| {C}_{{\ell }})=\displaystyle \frac{1}{\sqrt{{2}^{r}{\lambda }_{1}{\lambda }_{2}...{\lambda }_{r}}}\exp \left[-\displaystyle \frac{1}{2}{{\boldsymbol{S}}}^{T}{{\boldsymbol{C}}}^{\dagger }{\boldsymbol{S}}\right],\end{eqnarray} \tag{ 4 }$$

where  C denotes the N × N pixel-pixel CMB covariance matrix of  S at the chosen beam and pixel resolution. As discussed in Sudevan & Saha (2018), and as is the case in this article, the rank r of  C is lower than its size N, implying that  C is a singular matrix. ${{\boldsymbol{C}}}^{\dagger }$ therefore represents the Moore–Penrose generalized inverse (Moore 1920; Penrose 1955) of  C. The element of  C can be computed from the knowledge of the CMB theoretical angular power spectrum using Equation (6) of Sudevan & Saha (2018) by assuming that the CMB map is statistically isotropic. The set of ${\lambda }_{k},k\in \{1,2,...,r\}$ in the denominator of Equation (4) represents the nonzero eigenvalues of  C.

2.4. Drawing Samples of S

How do we draw samples of  S given  D and C_ℓ? We must do this without knowing or sampling the foreground components to keep our method independent of the foreground model. This would be possible if we could somehow remove all foregrounds without using their model, given  D and C_ℓ. The cleaned map obtained by using the global ILC method described in Sudevan & Saha (2018) can be used for exactly this purpose if we assume that the detector noise is negligible and one has a sufficient number of input frequency maps to remove all foreground components, as discussed in the current section. Let us consider the cleaned map,  Y, obtained by using linear combination of n frequency maps { X_i},

$$\begin{eqnarray}&&{\boldsymbol{Y}}=\sum _{i=1}^{n}{w}_{i}{{\boldsymbol{X}}}_{i},\end{eqnarray} \tag{ 5 }$$

where w_i represents the weight corresponding to the ith input frequency map. Clearly, we can neglect any detector noise contribution in  Y because  X_i themselves are assumed to contain negligible detector noise (e.g., see Equation (1)). Because the CMB follows a blackbody distribution, to preserve the CMB signal in the cleaned map, the weights for all frequency maps are constrained to add to unity, i.e., ${w}_{1}+{w}_{2}+{w}_{3}+...+{w}_{n}=1$. Minimizing the CMB covariance-weighted variance ${\sigma }^{2}={{\boldsymbol{Y}}}^{T}{{\boldsymbol{C}}}^{\dagger }{\boldsymbol{Y}}$ of the cleaned map  Y, subject to the above constraint on weights, as in Sudevan & Saha (2018), we obtain,

$$\begin{eqnarray}&&{\boldsymbol{W}}=\displaystyle \frac{{\hat{{\boldsymbol{A}}}}^{\dagger }{\boldsymbol{e}}}{{{\boldsymbol{e}}}^{T}{\hat{{\boldsymbol{A}}}}^{\dagger }{\boldsymbol{e}}},\end{eqnarray} \tag{ 6 }$$

where  W is an n × 1 column vector with the ith element given by the weight factor w_i.  e denotes an n × 1 column vector with all entries equal to unity, representing the frequency shape vector of the CMB component. Finally, the (i, j) element of the matrix $\hat{{\boldsymbol{A}}}$ can be computed in the pixel domain following

$$\begin{eqnarray}&&{\hat{A}}_{{ij}}={{\boldsymbol{X}}}_{i}^{T}{{\boldsymbol{C}}}^{\dagger }{{\boldsymbol{X}}}_{j}.\end{eqnarray} \tag{ 7 }$$

However, computing the above in pixel space is a numerically expensive process. One can considerably simplify Equation (7) in harmonic space. As shown in Appendix A, Equation (7) can conveniently be expressed in the multipole space following

$$\begin{eqnarray}&&{\hat{A}}_{{ij}}=\sum _{{\ell }=2}^{{{\ell }}_{\max }}\left(2{\ell }+1\right)\displaystyle \frac{{\hat{\sigma }}_{{\ell }}^{{ij}}}{{C}_{{\ell }}^{{\prime} }},\end{eqnarray} \tag{ 8 }$$

where ${C}_{{\ell }}^{{\prime} }$ represents the beam- and pixel-smoothed CMB theoretical power spectrum, ${C}_{{\ell }}^{{\prime} }={C}_{{\ell }}{B}_{{\ell }}^{2}{P}_{{\ell }}^{2}$, B_ℓ and P_ℓ being the beam and pixel window functions, respectively, and ${\hat{\sigma }}_{{\ell }}^{{ij}}$ represents the cross angular power spectrum between frequency maps  X_i and  X_j.⁵

We assume that there are n_f different foreground components, each with a constant spectral index all over the sky.⁶ We denote the shape vector of the kth foreground component by  f_k (with k ∈ {1, 2, ..., n_f}) , each one of which is an n × 1 column vector. Using Equations (1) and (6) in Equation (5), we obtain

$$\begin{eqnarray}&&{\boldsymbol{Y}}={\boldsymbol{S}}+\left[{{\boldsymbol{W}}}^{T}\sum _{k=1}^{{n}_{f}}{{\boldsymbol{f}}}_{k}\right]{{\boldsymbol{F}}}_{k}^{0},\end{eqnarray} \tag{ 9 }$$

where ${{\boldsymbol{F}}}_{k}^{0}$ is an N × 1 column vector representing an appropriately chosen template for the k foreground component. Equation (9) shows that the cleaned maps contain the pure CMB signal plus some foreground residual given by the second term. To find these residuals, introducing matrix notation, we first write Equation (8) as

$$\begin{eqnarray}&&\hat{{\boldsymbol{A}}}=\sum _{{\ell }=2}^{{{\ell }}_{\max }}\displaystyle \frac{\left(2{\ell }+1\right)}{{C}_{{\ell }}^{{\prime} }}{\hat{{\boldsymbol{\Sigma }}}}_{{\ell }},\end{eqnarray} \tag{ 10 }$$

where n × n data covariance matrix, ${\hat{{\boldsymbol{\Sigma }}}}_{{\ell }}$, in the harmonic space can be written in terms of the CMB angular power spectrum, ${\hat{C}}_{{\ell }}$, of the particular random realization under consideration and the foreground covariance matrix ${{\boldsymbol{C}}}_{{\ell }}^{f}$ as

$$\begin{eqnarray}&&{\hat{{\boldsymbol{\Sigma }}}}_{{\ell }}=\left[{\rm{\Gamma }}{{\boldsymbol{\Gamma }}}^{T}{\hat{C}}_{{\ell }}+{{\boldsymbol{C}}}_{{\ell }}^{f}\right]{B}_{{\ell }}^{2}{P}_{{\ell }}^{2}.\end{eqnarray} \tag{ 11 }$$

Using Equations (11) and (10) in Equation (6) and following a procedure similar to that in Saha & Aluri (2016), we obtain

$$\begin{eqnarray}&&{\boldsymbol{W}}=\displaystyle \frac{\left({\boldsymbol{I}}-{{\boldsymbol{C}}}_{f}{{\boldsymbol{C}}}_{f}^{\dagger }\right){\boldsymbol{e}}}{{{\boldsymbol{e}}}^{T}\left({\boldsymbol{I}}-{{\boldsymbol{C}}}_{f}{{\boldsymbol{C}}}_{f}^{\dagger }\right){\boldsymbol{e}}},\end{eqnarray} \tag{ 12 }$$

where  I denotes the n × n identity matrix and

$$\begin{eqnarray}&&{{\boldsymbol{C}}}_{f}=\sum _{{\ell }=2}^{{{\ell }}_{\max }}\displaystyle \frac{\left(2{\ell }+1\right)}{{C}_{{\ell }}}{{\boldsymbol{C}}}_{{\ell }}^{f}.\end{eqnarray} \tag{ 13 }$$

The product ${{\boldsymbol{C}}}_{f}{{\boldsymbol{C}}}_{f}^{\dagger }$ represents the projector on the column space, ${ \mathcal C }({{\boldsymbol{C}}}^{f})$, of  C_f. If we assume that n_f < n, which may be achieved by using a sufficiently large number of input frequency maps, the null space of  C_f is an nonempty set and ${\boldsymbol{I}}-{{\boldsymbol{C}}}_{f}{{\boldsymbol{C}}}_{f}^{\dagger }$ is a projector on this null space. Because the shape vector,  f_k of each foreground component with constant spectral indices completely lies on ${ \mathcal C }({{\boldsymbol{C}}}_{f})$, we must have ${{\boldsymbol{W}}}^{T}{{\boldsymbol{f}}}_{k}=0$ for all k. Therefore, from Equation (9) one finds that the foreground contamination in the final cleaned map at each pixel due to all foreground components disappears. Hence ${\boldsymbol{Y}}={\boldsymbol{S}}$.

Based upon the preceding discussions, to sample  S from ${P}_{1}({\boldsymbol{S}}| {C}_{{\ell }},{\boldsymbol{D}})$, we use Equation (5). The cleaned map in this case has the probability density as given by Equation (4), with the same covariance structure as mentioned therein.

2.5. Drawing Samples of C_ℓ

As discussed in Appendix B, the conditional density ${P}_{2}({C}_{{\ell }}| {\boldsymbol{S}},{\boldsymbol{D}})$ can be written as

$$\begin{eqnarray}&&{P}_{2}\left({C}_{{\ell }}| {\hat{C}}_{{\ell }}\right)\propto {\left(\displaystyle \frac{1}{{C}_{{\ell }}}\right)}^{\left(2{\ell }+1\right)/2}\exp \left[-\displaystyle \frac{{\hat{C}}_{{\ell }}\left(2{\ell }+1\right)}{2{C}_{{\ell }}}\right],\end{eqnarray} \tag{ 14 }$$

where the variable $z={\hat{C}}_{{\ell }}\left(2{\ell }+1\right)/{C}_{{\ell }}$ follows a χ² distribution of 2ℓ − 1 degree of freedom. To draw samples of C_ℓ from Equation (14), we first draw z from the χ² distribution of 2ℓ − 1 degrees of freedom, which is achieved by drawing 2ℓ − 1 independent standard normal deviates and forming the sum of their squares. Given the value of ${\hat{C}}_{{\ell }}$ estimated from the map, we then find C_ℓ following ${C}_{{\ell }}={\hat{C}}_{{\ell }}\left(2{\ell }+1\right)/z$.

We use 10 WMAP nine-year difference assembly (DA) maps and 7 Planck 2015 maps, three of the latter are at LFI frequencies (30, 40, and 70 GHz) and the rest are at four HFI frequencies (100, 143, 217, and 353 GHz). The processing of input maps remains identical to Sudevan & Saha (2018) and results in a total of 12 input maps, 5 at WMAP and 7 at Planck frequencies. We note that because we are interested in an analysis over large angular scales of the sky, where detector noise can be ignored, we chose a low pixel resolution defined by the HEALPix⁷ parameter N_side = 16 and a Gaussian beam-smoothing of 9° full width at half-maximum (FWHM) for each input map. We remove both monopole and dipole from all the full-sky input maps before the analysis. To sample the posterior density $P({\boldsymbol{S}},{C}_{{\ell }}| {\boldsymbol{D}})$, we simulate a total of 10 Gibbs chains, each containing 5000 joint samples of cleaned maps and theoretical power spectrum, following the three sampling steps described in Section 2. At any given chain and at any given iteration, to sample  S, we use Equation (5), where the weights are described by the vector  W defined by Equation (6). The elements of matrix  A that appears in Equation (6) are computed following Equation (8) using the last sampled C_ℓ values. After sampling  S, we estimate its full-sky power spectrum , ${\hat{C}}_{{\ell }}$, which we use to obtain a new sample of C_ℓ using the method described in Section 2.5. We emphasize that both ${\hat{C}}_{{\ell }}$ and C_ℓ contain beam- and pixel-smoothing effects absorbed in them. We note in passing that the weights as given by Equation (6) are insensitive to such smoothing, however, because in Equation (8), both C_ℓ and ${\hat{\sigma }}_{{\ell }}^{{ij}}$ contain same smoothing effects. 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 (Planck Collaboration et al. 2016c), where ΔC_ℓ denotes the error due to cosmic variance alone.

All 10 chains generate a total of 50,000 joint samples of a cleaned map and a theoretical CMB power spectrum. The burn-in phase in each chain is very brief. Visually, this phase does not appear to contain more than a few Gibbs iterations. However, we remove 50 initial Gibbs iterations from each chain as a conservative estimate of the burn-in period. After the burn-in rejection, we have a total of 49,500 samples from all chains.

4.1. Cleaned Maps

Using all samples after burn-in rejection, we estimate the marginalized probability density of the CMB temperature at each pixel given the observed data. The normalized probability densities obtained by division by the corresponding mode of the marginalized density function for some selected pixels over the sky are shown in Figure 1. These density functions are approximately symmetric, with some visible asymmetry near the tails. The positions of the mean temperatures are shown by the vertical blue lines for each pixel of this plot. We estimate the best-fit cleaned CMB map by taking the pixel temperatures corresponding to the location of the mode of the density for each pixel at N_side = 16. We show the best-fit map at the top panel of Figure 2. We compare the best-fit map with other cleaned CMB maps that are obtained using different methods by other science groups. We show the differences of the best-fit map from Planck Commander and SMICA-cleaned maps (Planck Collaboration et al. 2016a) in the middle left and middle right panels of the same figure, respectively, at 9° Gaussian beam resolution. The bottom left and bottom right panel show the difference of the best-fit cleaned map from needlet space ILC (NILC) and SEVEM-cleaned maps, respectively. We note that on large angular scales, the Commander method is designed to be optimal for reconstructing the CMB signal. We emphasize the striking similarity between the best-fit and NILC CMB map. The SEVEM-cleaned map contains significant positive residuals at the galactic plane. Interestingly, the best-fit map contains somewhat lower pixel temperatures at isolated locations along the galactic plane than the Commander-, NILC-, and SMICA-cleaned maps. In summary, the minor differences between our best-fit map and other Planck-cleaned maps (except for SEVEM) are mostly in the galactic region, where the foreground contamination dominates.⁸

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figure 1 we see that the mean and marginalized posterior maximum of the CMB temperature at different pixels agree closely with each other. We estimate the mean CMB map using all 49,500 samples for each pixel and show this in the top panel of Figure 3. The mean map matches the best-fit CMB map shown in the top panel of Figure 2 very well. We have plotted the difference between the best-fit and mean map in the bottom panel of Figure 3. Both maps agree with each other within an absolute difference of 1 μK. In order to quantify the reconstruction error in the cleaned CMB map obtained after each iteration of Gibbs sampling, we generate a standard deviation map using all 49,500 cleaned maps. We show this map in the bottom panel of Figure 3. From this panel we see that the reconstruction error is very small all over the sky. The maximum reconstruction error is visible along the galactic plane and toward the center of our galaxy, where the input frequency maps contain strong foreground contamination.

Zoom In Zoom Out Reset image size

4.2. Angular Power Spectrum

We estimate the marginalized probability density of the CMB theoretical angular power spectrum and show the results in Figure 4 for different multipoles. Like the density functions of the pixel temperatures, as discussed in Section 4.1, we normalize these densities to a value of unity at their peaks. The horizontal axis of each plot of this figure represents ${\ell }\left({\ell }+1\right){C}_{{\ell }}/(2\pi )$ in units of 1000 μK². The density functions show long asymmetric tails for low multipoles (e.g., ℓ = 2, 3, 5). For large multipoles (ℓ ≥ 20), the asymmetry of the densities is gradually reduced. The region within the two vertical lines in each plot shows the 1σ (68.27%) confidence interval for the CMB theoretical power spectrum for the corresponding multipole.

Zoom In Zoom Out Reset image size

In the top panel of Figure 5 we show the best-fit theoretical CMB angular power spectrum in brown, defined by the positions of the peaks of the marginalized angular power spectrum density functions (e.g., Figure 4). The asymmetric error bars at each ℓ show the 1σ confidence interval for the theoretical angular power spectrum. The black line shows the theoretical power spectrum consistent with Planck 2015 results (Planck Collaboration et al. 2016c). The best-fit theoretical angular power spectrum agrees well with the spectra estimated from Commander-, NILC-, and SMICA-cleaned maps, which are shown by purple, dark blue, and sky blue points, respectively. We see that due to the presence of residual foregrounds in the galactic plane, the angular power spectrum estimated from the SEVEM-cleaned map (pink points) has a larger power at all multipoles than our best-fit power spectrum. In the bottom panel of Figure 5 we show the difference of the best fit with other Planck power spectra. We observe that within the cosmic variance error, our power spectrum agrees well with Planck Commander, NILC, and SMICA power spectra. Our best-fit theoretical power spectrum agrees very well with the spectrum estimated from the NILC CMB map.

Zoom In Zoom Out Reset image size

Each of the Gibbs sampling chains for the estimation of the CMB signal and its theoretical angular power spectrum joint posterior consists of 4950 samplings after rejection of the burn-in phase. A diagnosis is necessary to be certain that these samples have converged to the actual targeted CMB posterior. This is a condition (the convergence) when satisfied an inference that is drawn about any parameter using the chains does not depend upon the point where the chain starts. Gelman & Rubin (1992) propose that lack of any such convergence is better diagnosed when we simulate a set of "parallel" chains than a single chain. Using all the 10 Gibbs chains, we therefore check for convergence by using the Gelman–Rubin statistic (Gelman & Rubin 1992). A detailed description of the statistic is given in Gelman & Rubin (1992) and Brooks & Gelman (1998), but for completeness, we define the statistic below.

Let us assume that we have generated M number of different chains, and let L be the number of steps in each chain after rejection of samples during the burn-in period.⁹ For a model parameter θ, let us assume that the sample posterior mean is given by ${\bar{\theta }}_{m}$ for mth chain using all L samples. Let the corresponding sample posterior variance be ${\bar{\sigma }}_{m}^{2}$. Then the between-chain (B/L) and within-chain variances (W) are given by

$$\begin{eqnarray}&&B=\displaystyle \frac{L}{M-1}\sum _{m=1}^{M}{\left({\bar{\theta }}_{m}-\bar{\theta }\right)}^{2}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&W=\displaystyle \frac{1}{M}\sum _{m=1}^{M}{\bar{\sigma }}_{m}^{2},\end{eqnarray} \tag{ 16 }$$

respectively, where $\bar{\theta }$ is the overall posterior mean of the samples estimated from all M chains and is given by $\bar{\theta }=\tfrac{1}{M}{\sum }_{m=1}^{M}{\theta }_{m}$. We define the pooled posterior variance following

$$\begin{eqnarray}&&\hat{V}=\displaystyle \frac{L-1}{L}W+\displaystyle \frac{M+1}{{ML}}B,\end{eqnarray} \tag{ 17 }$$

which can be used to compute the Gelman–Rubin statistic R as

$$\begin{eqnarray}&&R=\sqrt{\displaystyle \frac{\hat{V}}{W}}.\end{eqnarray} \tag{ 18 }$$

Following Gelman & Rubin (1992) and Brooks & Gelman (1998), a value of R close to unity implies that each of M Gibbs chains have converged to the target posterior density.

We have plotted the Gelman–Rubin statistic, R, for the theoretical angular power spectrum samples for the multipole range 2 ≤ ℓ ≤ 32 in the top panel of Figure 6. The value of R lies well within 0.9999 and 1.0003, implying convergence. In the bottom panel of Figure 6 we show the map of R − 1 all over the sky. R lies with in 0.99992 and 1.0002 for all the pixels, again implying convergence of the Gibbs chains.

Zoom In Zoom Out Reset image size

We generate a set of input maps at the 12 different WMAP and Planck frequencies at a Gaussian beam resolution 9° and N_side = 16 following the same procedure as described in Sudevan & Saha (2018). We do not reproduce the method in that article and refer to the above article for a description of the input frequency maps. We note that in the current work, we need to simulate only one random realization of the 12 input frequency maps. The random CMB realization used in the input frequency maps is generated using the CMB theoretical angular power spectrum consistent with Planck 2015 results. As is the case for our analysis of the Planck and WMAP observations, we remove both monopole and dipole from all the simulated input maps before sampling from the posterior density $P({\boldsymbol{S}},{C}_{{\ell }}| {\boldsymbol{D}})$. We simulate a total of 10 Gibbs chains. To draw the first cleaned map sample, we initialize the theoretical CMB power spectrum uniformly within ±3Δ_ℓ of the true theoretical spectrum, where Δ_ℓ represents the error induced by cosmic variance. As in the case of the analysis of the WMAP and Plank observed maps, for simulations, also we find that the burn-in period ends very rapidly. In particular, from the trace plots of the pixel temperature of the cleaned maps, we see that the burn-in phase is complete within a few samples. As a conservative approach, however, we reject the initial 50 samples from each Gibbs chain. After burn-in rejection, we have a set of 4950 joint samples of cleaned maps and theoretical angular power spectra from each Gibbs chain.

Using all cleaned map samples from all chains after burn-in rejection, we form a marginalized density of the CMB temperature at each N_side = 16 pixel. A CMB map formed from the pixel temperatures corresponding to the modes of these density functions defines the best-fit cleaned CMB map obtained from the simulation. We show the best-fit cleaned map in the top panel of Figure 7. The difference of the best-fit and input CMB realization is shown in the second panel of the same figure. Clearly, the best-fit CMB map matches the CMB realization used in the simulation very well. The maximum difference (15.15 μK) between the two maps is observed along the galactic plane. This shows that our method reliably removes the foreground. The third panel of Figure 7 shows the difference of the best-fit and mean CMB maps obtained from all Gibbs samples. These two maps agree very well with each other. The last panel of Figure 7 shows the standard deviation map computed from the Gibbs samples. The maximum error of 4.2 μK is observed at the galactic center. In summary, using the Monte Carlo simulations of our method, we see that the best-fit and mean CMB maps agree very well with the input CMB map, indicating that a reliable foreground minimization can be achieved with our method.

Zoom In Zoom Out Reset image size

We show the best-fit estimate of underlying CMB theoretical angular power spectrum from Monte Carlo simulations in Figure 8 in brown line. The asymmetric error-limits shown on this power spectrum indicate the 68.27% confidence intervals. The angular power spectrum for the input CMB map used in the simulation is shown in blue. The best-fit estimate agrees nicely with the input angular power spectrum The black line of this figure represents the underlying theoretical angular power spectrum that is used to generate the input CMB realization for the simulation.

Zoom In Zoom Out Reset image size

We test the Gibbs sequences obtained from the simulations for convergence as in Section 5 using the Gelman–Rubin statistic, R. The maximum and minimum values of R for the sampled CMB maps over all pixels are 1.00017 and 0.999924, respectively. The corresponding values for the sampled angular power spectra are 1.00027 and 0.999935, respectively. These values of R close to unity indicate convergence of the Gibbs sequences in the Monte Carlo simulations.

In a CMB experiment, the observed frequency maps are calibrated using various techniques (Bennett et al. 2003; Weiland et al. 2011; Adam et al. 2016; Ade et al. 2016; Aghanim et al. 2018; Akrami et al. 2018). Although these methods are advanced and they are already applied on the data, some residual calibration uncertainties exist in each frequency map. It is important to investigate the reconstruction error and bias effects that are translated into the final cleaned CMB signal and its angular power spectrum when such residual calibration errors are present in the input maps.

In order to understand the impact of these residual calibration errors on the cleaned CMB map and the angular power spectrum obtained by our method, we perform an additional two sets of independent Monte Carlo simulations. In each of these two sets to simulate the presence of residual calibration errors in all input frequency maps we multiply them by the actual calibration factors obtained by drawing the unit mean Gaussian random variable x_i (i = 1, .., n, n being the total number of input frequency maps), with the standard deviation equal to the chosen amount of calibration error. We chose a 0.5% error for each input frequency map for the first simulation set and a 1.0% error for the second simulation. This choice of error levels is well above the actual residual calibration errors expected from the WMAP and Planck satellite missions and hence allows us to understand the performance of our method when it is passed in simulations through a rather critical environment—in presence of unrealistically high levels of residual calibration errors.¹⁰ For the two sets of Monte Carlo simulations of the CMB reconstruction, we (incorrectly) assume that all input frequency maps are calibrated correctly, i.e., no residual calibration errors are present in order to propagate the effects of residual calibration errors to the final CMB products.

In Figure 9 we compare the cleaned CMB maps obtained from the two sets of simulations containing calibration errors with the corresponding results where calibration coefficients are exactly known. We note that the input frequency maps for these three sets of simulations contain the same CMB realization. Moreover, the noise realization for any given frequency for all the three sets is identical as well. In the top left panel of Figure 9 we show the difference between the best-fit cleaned CMB maps obtained from the third and first set of simulations. The bottom left subpanel shows the difference of best-fit cleaned CMB maps obtained using the third and second set of simulations. An increase (decrease) in maximum (minimum) pixel values of magnitude ∼2 μK in the bottom left panel is seen relative to the top left panel, indicating a larger reconstruction error in the best-fit cleaned map when the calibration uncertainty in the input frequency maps is increased from 0.5% to 1%. A similar trend is also visible for the difference of the mean cleaned CMB maps obtained using the two erroneous sets of frequency maps from the no-error case (top middle and bottom middle panels). In the top right panel we show the ratio of the standard deviation maps obtained from the 0.5% calibration error case and from the no-calibration-error case. Similarly, the bottom right panel shows the ratio of the standard deviation maps estimated from simulations with a 1% calibration error and without error. Most of the pixels of the top right panel have values close to unity, except for some localized regions near the galactic plane. Clearly, increasing the calibration error from 0% to 0.5% causes some pixels near the galactic plane to be erroneous by about 5% compared to the baseline values when input frequency maps are correctly calibrated. Increasing the calibration errors in the input frequency maps further to 1% causes more pixels in the cleaned CMB map to be erroneous at the same 5%, as shown in the bottom right panel. When we consider the maximum pixel error of 4.2 μK from the standard deviation map shown in the bottom panel of Figure 7, the error in the reconstructed CMB map obtained with our method increases only by less than 4.2 × 0.05 μK ∼ 0.2 μK even for a 1% calibration error in the input maps from the case when input frequency maps are exactly calibrated.

Zoom In Zoom Out Reset image size

We show the best-fit theoretical CMB angular power spectra in the top panel of Figure 10 for Monte Carlo simulations with 0.5% (blue line), 1% (green line), and no calibration error (red line). Even though for certain mulitpoles the estimated CMB power spectrum values from simulations with calibration errors differ slightly from the corresponding values with no calibration error, any such differences lie well inside the cosmic variance error. In order to compute these errors, we define relative error using the best-fit angular power spectra following

$$\begin{eqnarray}&&\mathrm{Err}=\displaystyle \frac{{C}_{{\ell }}^{0}-{C}_{{\ell }}^{i}}{{C}_{{\ell }}^{0}},\end{eqnarray} \tag{ 19 }$$

where ${C}_{{\ell }}^{0}$ and ${C}_{{\ell }}^{i}$ are the best-fit power spectrum obtained form simulations with no calibration error and with some calibration errors, respectively. In the bottom panel of Figure 10, we show the relative error in the angular power spectra for 0.5% and 1.0% calibration errors in the input frequency maps. We observe that with increase in calibration errors, the relative error for different mulitpoles varies somewhat randomly, without any indication of a significant systematic bias of the best-fit power spectrum. Based upon the simulation results, we safely conclude that for a realistic level of the calibration uncertainty in input frequency maps, the error budget in the cleaned CMB map and the angular power spectrum obtained by the global ILC method shows negligible increment. Furthermore, we do not see any indication of significant bias in the cleaned CMB maps and angular power spectrum due to presence of such calibration uncertainties. We defer a more detailed study of the issue of the calibration uncertainty in the input frequency maps to an accompanying publication, Sudevan & Saha (2020).

Zoom In Zoom Out Reset image size

In this article, we have presented a new method for estimating the CMB posterior density over the large angular scales of the sky, given the Planck and WMAP observations, using a global ILC method (Sudevan & Saha 2018) and Gibbs sampling (German & German 1984) as the basic tools. Our main results are joint estimates of the best-fit CMB signal and its theoretical angular power spectrum, along with the appropriate confidence intervals that can be directly used to estimate cosmological parameters. Therefore, our work effectively extends for the first time the ILC method for such purposes. We sample the CMB signal at each Gibbs iteration, conditioned on a set of CMB theoretical angular power spectra obtained in the previous Gibbs iteration. The CMB reconstruction step is independent of any explicit model of foreground components—which is a characteristic of the usual ILC method. However, considering the sampling of both CMB signal and its theoretical angular power spectrum, the new method extends the model-independent nature of the CMB reconstruction of the ILC method to the entire posterior density estimation at low resolution. Thus our method serves as a complementary route to the CMB posterior estimation, where detailed models of foregrounds are taken into account. We have implemented the CMB reconstruction method in the harmonic space, which significantly reduces computational time, unlike the pixel-space approach of Sudevan & Saha (2018).

There are some aspects of the method that need to be addressed in future investigations. In the current work we have assumed that the detector noise can be completely ignored, which is a valid assumption on the large angular scales for experiments such as WMAP and Planck. A general framework will be to formalize the method in the presence of detector noise. In the presence of detector noise, the blind foreground removal procedure will leave some foreground residuals on the cleaned maps. It would be interesting to see whether these residuals can be taken care of in a manner independent of the foreground model.

In the current method, where detector noise is assumed to be negligible, residual foregrounds will be present in the cleaned maps if the effective number of foreground components n_f present in the input frequency maps becomes larger than or equal to the number of the input frequency maps (n) available. For a large angular scale analysis like the one of this paper, n_f < n is a reasonable assumption outside the galactic plane. Along the plane, where the foreground spectral properties are expected to show a larger variation than outside plane, the effect of such residuals is mitigated by the smoothing of the input sky maps over the large angular scales. By performing detailed Monte Carlo simulations, we see that the method leaves a small residual along the galactic plane. By comparing our cleaned map and angular power spectrum results with those obtained by other science groups, we show that the level of such foreground residuals is low and of comparable magnitude to those present in CMB maps obtained by other methods. Using Monte Carlo simulations, we verify that the presence of residual calibration uncertainty in the input frequency maps does not cause any significant error or bias in the cleaned CMB map and the angular power spectrum obtained by using the method of this article.

We thank the anonymous referee of an earlier publication by the authors (Sudevan & Saha 2018) for suggestions to integrate Gibbs sampling into our method. Our work is based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We use the publicly available HEALPix (Górski et al. 2005) package (http://healpix.sourceforge.net) for some of the analysis of this work. We acknowledge the use of thePlanck 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.

Let us assume that S(p) denotes a random simulation of a pixellized CMB map (p denotes the pixel index) at some beam and pixel resolutions. When we use the spherical harmonic decomposition $S(p)={\sum }_{{\ell },m}{a}_{{\ell },m}{Y}_{{\ell }m}(p)$, the (p, q) element C_pq of the pixel-pixel CMB covariance matrix ${\boldsymbol{C}}$ can be written as

$$\begin{eqnarray}&&{C}_{{pq}}=\left\langle S(p)S(q)\right\rangle =\sum _{{\ell }m{\ell }^{\prime} m^{\prime} }\left\langle {a}_{{\ell }m}{a}_{{\ell }^{\prime} m^{\prime} }^{* }\right\rangle {Y}_{{\ell }m}(p){Y}_{{\ell }^{\prime} m^{\prime} }^{* }(q),\end{eqnarray} \tag{ A1 }$$

where $\left\langle ...\right\rangle$ represents the ensemble average and we assume that the beam- and pixel-smoothing effects are implicitly contained in the spherical harmonic coefficients ${a}_{{\ell }m}$. Using the statistical isotropy of the CMB, namely, $\left\langle {a}_{{\ell }m}{a}_{{\ell }^{\prime} m^{\prime} }^{* }\right\rangle ={C}_{{\ell }}^{{\prime} }{\delta }_{{\ell }{\ell }^{\prime} }{\delta }_{{mm}^{\prime} }$, where ${C}_{{\ell }}^{{\prime} }={C}_{{\ell }}{B}_{{\ell }}^{2}{P}_{{\ell }}^{2}$, (B_ℓ and P_ℓ are the beam and pixel window functions, respectively), we obtain

$$\begin{eqnarray}&&{C}_{{pq}}=\sum _{{\ell }m}{C}_{{\ell }}^{{\prime} }{Y}_{{\ell }m}(p){Y}_{{\ell }m}^{* }(q).\end{eqnarray} \tag{ A2 }$$

In matrix notation, we write Equation (A2) as

$$\begin{eqnarray}&&{\boldsymbol{C}}=\sum _{{\ell }m}{C}_{{\ell }}^{{\prime} }{{\boldsymbol{D}}}_{{\ell }m},\end{eqnarray} \tag{ A3 }$$

where the elements of matrix  D_{ℓ m} are given by ${D}_{{\ell }m}^{{pq}}={Y}_{{\ell }m}(p){Y}_{{\ell }^{\prime} m^{\prime} }^{* }(q)$. We can write ${{\boldsymbol{D}}}_{{\ell }m}={{\boldsymbol{Z}}}_{{\ell }m}{{\boldsymbol{Z}}}_{{\ell }m}^{c}$, where ${{\boldsymbol{Z}}}_{{\ell }m}$ is an N × 1 column vector with elements ${Z}_{{\ell }m}(p)={Y}_{{\ell }m}(p)$, and the superscript ${}^{c}$ represents the Hermitian conjugate. It is worth emphasizing that the right-hand side of the above equation represents a linear combination of ${{\boldsymbol{D}}}_{{\ell }m}$ matrices with the scalar amplitudes given by ${C}_{{\ell }}^{{\prime} }$. Clearly, therefore,

$$\begin{eqnarray}&&{{\boldsymbol{C}}}^{\dagger }=\sum _{{\ell }m}{{\boldsymbol{D}}}_{{\ell }m}^{\dagger }{C}_{{\ell }}^{{\prime} \dagger }.\end{eqnarray} \tag{ A4 }$$

Using the definition of the Moore–Penrose generalized inverse ${{\boldsymbol{x}}}^{\dagger }$ of a vector ${\boldsymbol{x}}$, ${{\boldsymbol{x}}}^{\dagger }={{\boldsymbol{x}}}^{c}/| | {\boldsymbol{x}}| {| }^{2}$, where $| | ...| {| }^{2}$ represents the squared norm of the vector, we obtain

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{{\ell }m}^{\dagger }=\displaystyle \frac{{{\boldsymbol{Z}}}_{{\ell }m}{{\boldsymbol{Z}}}_{{\ell }m}^{c}}{| | {{\boldsymbol{Z}}}_{{\ell }m}| {| }^{4}}.\end{eqnarray} \tag{ A5 }$$

When we use the orthogonality of spherical harmonics over discrete HEALPix pixels,

$$\begin{eqnarray}&&\left[\sum _{p=1}^{N}{Y}_{{\ell }m}(p){Y}_{{\ell }m}^{* }(p)\right]\displaystyle \frac{4\pi }{N}=1,\end{eqnarray} \tag{ A6 }$$

it is easy to find $| | {{\boldsymbol{Z}}}_{{\ell }m}| {| }^{4}={\left(N/(4\pi \right)}^{2}$. Using this result and Equation (A5) in Equation (A4), we obtain element wise

$$\begin{eqnarray}&&{C}_{{pq}}^{\dagger }=\sum _{{\ell }m}\displaystyle \frac{1}{{C}_{{\ell }}^{{\prime} }}\displaystyle \frac{{Y}_{{\ell }m}(p){Y}_{{\ell }m}^{* }(q)}{{\left(\tfrac{N}{4\pi }\right)}^{2}},\end{eqnarray} \tag{ A7 }$$

where we have used ${C}_{{\ell }}^{{\prime} \dagger }=1/{C}_{{\ell }}^{{\prime} }$. Expanding the input frequency maps ${{\boldsymbol{X}}}_{i}$ in spherical harmonic space, ${X}_{i}(p)={\sum }_{{{\ell }}_{1}{m}_{1}}{a}_{{{\ell }}_{1},{m}_{1}}^{i}{Y}_{{{\ell }}_{1}{m}_{1}}(p)$, using Equation (A7) and the orthogonality condition of spherical harmonics as mentioned in Equation (A6), after some algebra, we obtain

$$\begin{eqnarray}&&{\hat{A}}_{{ij}}={{\boldsymbol{X}}}_{i}^{T}{{\boldsymbol{C}}}^{\dagger }{{\boldsymbol{X}}}_{j}=\sum _{{\ell }}\left(2{\ell }+1\right)\displaystyle \frac{{\hat{\sigma }}_{{\ell }}^{{ij}}}{{C}_{{\ell }}^{{\prime} }},\end{eqnarray} \tag{ A8 }$$

where

$$\begin{eqnarray}&&{\hat{\sigma }}_{{\ell }}^{{ij}}=\sum _{m=-{\ell }}^{{\ell }}{a}_{{\ell }m}^{i}{a}_{{\ell }m}^{j* }/\left(2{\ell }+1\right).\end{eqnarray} \tag{ A9 }$$

We note that if x₁, x₂, ..., x_μ are identically distributed and independent Gaussian random variables with zero mean and unit variance, the new variable $x={\sum }_{k=1}^{\mu }{x}_{k}^{2}$ is distributed as a χ² random variable with μ degrees of freedom, with the probability density given by

$$\begin{eqnarray}&&P(x)=\displaystyle \frac{1}{{2}^{\nu }{\rm{\Gamma }}(\nu )}{x}^{\nu -1}\exp \left[-\displaystyle \frac{x}{2}\right],\end{eqnarray} \tag{ B1 }$$

where ν = μ/2. With this definition, the variable $x\equiv \left(2{\ell }+1\right)\tfrac{{\hat{C}}_{{\ell }}}{{C}_{{\ell }}}$ is distributed as Equation (B1), where ${\hat{C}}_{{\ell }}$ and C_ℓ denote the realization-specific and theoretical CMB angular power spectrum, respectively. To find the density function of ${\hat{C}}_{{\ell }}$, we first note using Equation (B1) that the density function Q(y) for the transformed variable y = βx (β = constant) follows $Q(y)=P(x){dx}/{dy}$, where on the right-hand side, x must be replaced by y using the inverse transformation $x=y/\beta$, so that we obtain a function of y as required. Using this concept and defining $y\equiv {\hat{C}}_{{\ell }}={C}_{{\ell }}x/(2{\ell }+1)$, so that $\beta ={C}_{{\ell }}/(2{\ell }+1)$, we find

$$\begin{eqnarray}\begin{array}{rcl}Q({\hat{C}}_{{\ell }}) & = & {\left[{2}^{\left(2{\ell }+1\right)/2}{\rm{\Gamma }}\left(\displaystyle \frac{2{\ell }+1}{2}\right)\right]}^{-1}{\left(\displaystyle \frac{2{\ell }+1}{{C}_{{\ell }}}\right)}^{\left(2{\ell }+1\right)/2}\\ & & \times \,{\hat{C}}_{{\ell }}^{\left(2{\ell }+1\right)/2-1}\exp \left[-\displaystyle \frac{{\hat{C}}_{{\ell }}\left(2{\ell }+1\right)}{2{C}_{{\ell }}}\right].\end{array}\end{eqnarray} \tag{ B2 }$$

When C_ℓ is assumed as a random variable, Equation (B2) represents the conditional probability density $Q\left({\hat{C}}_{{\ell }}| {C}_{{\ell }}\right)$. With Bayes theorem and a uniform prior on C_ℓ up to some irrelevant constant, the probability density of C_ℓ given some ${\hat{C}}_{{\ell }}$ can be obtained as

$$\begin{eqnarray}&&R\left({C}_{{\ell }}| {\hat{C}}_{{\ell }}\right)\propto {\left(\displaystyle \frac{1}{{C}_{{\ell }}}\right)}^{\left(2{\ell }+1\right)/2}\exp \left[-\displaystyle \frac{{\hat{C}}_{{\ell }}\left(2{\ell }+1\right)}{2{C}_{{\ell }}}\right].\end{eqnarray} \tag{ B3 }$$

Now defining a new variable $z={\hat{C}}_{{\ell }}\left(2{\ell }+1\right)/{C}_{{\ell }}$ and noting that the exponent of $1/{C}_{{\ell }}$ in Equation (B3) can be written as $(2{\ell }+1)/2=(2{\ell }-1)/2+1$, we can write Equation (B3) as

$$\begin{eqnarray}&&R\left(z| {\hat{C}}_{{\ell }}\right)\propto {z}^{-(2{\ell }-1)/2-1}\exp \left[-\displaystyle \frac{z}{2}\right],\end{eqnarray} \tag{ B4 }$$

where we have omitted some irrelevant constants. Comparing Equation (B4) with Equation (B1), we readily identify Equation (B4) as a χ² distribution of $2{\ell }-1$ degrees of freedom in the variable z.