# Astrophysical component separation

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

## CMBCosmic Microwave background and foreground separation

See the Component Separation paper Planck-2013-XII[1] Planck-2015-A09[2] for details.

### NILC

NILC is a linear method for combining the input frequency channels. It implements an ILC with weighting coefficients varying over the sky and over the multipole range up to $\ell=3200$ and it does so using 'needlets' which are spherical wavelets. A special procedure is used for processing the coarsest needlet scale which contains the large scale multipoles.

In practice, our NILC processing depends on several implementation choices, as follows:

Input channels
In this work, the NILC algorithm is applied to all Planck channels from 44 to 857 GHz omitting only the 30 GHz channel.
Pre-processing of point sources
Identical to the SMICA pre-processing.
The NILC map is actually produced in a three-step process. In a first step, the NILC weights are computed from covariance matrices evaluated using a Galactic mask removing about 2 % of the sky (and is apodized at 1◦). In a second step, those NILC weights are applied to needlet coefficients computed over the complete sky (except for point source masking/subtraction), yielding a NILC estimate over the full sky (except for the point source mask). In short, the weights are computed over a masked sky but are applied to a full sky (excluding point sources). In a final step, the pixels masked due to point source processing are replaced by the values of a constrained Gaussian realization (inpainting).
Spatial localization
The boundaries of the zones used for spatial localisation are obtained as iso-level curves of a low resolution map of Galactic emission.
Beam control and transfer function
As in the SMICA processing, the input maps are internally re-beamed to a 5′ resolution, so the resulting map is automatically synthesized with an effective Gaussian beam of five arcminutes, according to the unbiased nature of the ILC.
Using SMICA recalibration
In our current implementation, the NILC solution uses the values determined by SMICA for the spectrum.

### SEVEM

The aim of Sevem is to produce clean maps at one or several frequencies by using a procedure based on template fitting. The templates are internal, i.e., they are constructed from Planck data, avoiding the need for external data sets, which usually complicates the analyses and may introduce inconsistencies. The method has been successfully applied to Planck simulations (Leach et al., 2008[3]) and to WMAP polarisation data (Fernandez-Cobos et al., 2012[4]). In the cleaning process, no assumptions about the foregrounds or noise levels are needed, rendering the technique very robust.

The input maps used are all the Planck frequency channels. In particular, for intensity, we have cleaned the 100, 143 GHz and 217 GHz maps using four templates. Three or them are constructed as the difference of the following Planck channels (smoothed to a common resolution to remove the contribution): (30-44) GHz, (44-70) GHz, (545-353) GHz and a fourth template given by the 857 GHz channel (smoothed at the resolution of the 545 GHz channel). For polarization we clean maps at frequencies of 70, 100 and 143 GHz using three templates for each channel. In particular, we use (30-44) GHz smoothed to a common resolution, (353-217) at 10', and (217-143) GHz at 1 degree resolution to clean the 70 and 100 GHz maps. To clean the 143 GHz channel, the last template is replaced by (217-100) GHz at 1 degree resolution. Before constructing the templates, for both intensity and polarization, we perform inpainting in the positions of detected point sources to reduce their contamination of the final map.

A linear combination of the templates is then subtracted from the Planck sky map at the considered frequency, in order to produce the clean map. The coefficients of the linear combination are obtained by minimising the variance of the clean map outside a given mask. Although we exclude very contaminated regions during the minimization, the subtraction is performed for all pixels and, therefore, the cleaned maps cover the full-sky (although we expect that foreground residuals are present in the excluded areas). Inpainting of point sources is also carried out in the clean maps.

The final intensity map has then been constructed by combining the 143 and 217 GHz cleaned maps by weighting the maps in harmonic space taking into account the noise level, the resolution and a rough estimation of the foreground residuals of each map (obtained from realistic simulations). This final map has a resolution corresponding to a Gaussian beam of =5 arcminutes at $N_{side}$=2048. The final polarization map has been obtained by combining the 100 and 143 GHz clean maps at $N_{side}$=1024 and has a resolution of 10 arc minutes.

### SMICA

A linear method, SMICA reconstructs a map as a linear combination in the harmonic domain of $N_{chan}$ input frequency maps with weights that depend on multipole $\ell$. Given the $N_{chan} × 1$ vector $\mathbf{x}_{\ell m}$ of spherical harmonic coefficients for the input maps, it computes coefficients $s_{\ell m}$ for the map as

$\label{eq:smica:shat} \hat{s}_{\ell m} = \mathbf{w}^†_\ell \mathbf{x}_{\ell m}$

where the $N_{chan} × 1$ vector $\mathbf{w}_\ell$ which contains the multipole-dependent weights is built to offer unit gain to the with minimum variance. This is achieved with

$\label{eq:smica:w} \mathbf{w}_\ell = \frac{\mathbf{R}_\ell ^{-1} \mathbf{a}}{\mathbf{a}^† \mathbf{R}_\ell ^{-1} \mathbf{a}}$

where vector $\mathbf{a}$ is the emission spectrum of the evaluated at each channel (allowing for possible inter-channel recalibration factors) and $\mathbf{R}_\ell$ is the $N_{chan} × N_{chan}$ spectral covariance matrix of $\mathbf{x}_{\ell m}$. Taking $\mathbf{R}_\ell$ in the second equation to be the sample spectral covariance matrix $\mathbf{\hat{R}}_\ell$ of the observations:

$\label{eq:smica:Rhat} \mathbf{\hat{R}}_\ell = \frac{1}{2 \ell + 1} \sum_m \mathbf{x}_{ \ell m} \mathbf{x}_{\ell m}^†$

would implement a simple harmonic-domain ILC. This is not what SMICA does. As discussed below, we instead use a model $\mathbf{R}_\ell (θ)$ and determine the covariance matrix to be used in the second equation by fitting $\mathbf{R}_\ell (θ)$ to $\mathbf{\hat{R}}_\ell$. This is done in the maximum likelihood sense for stationary Gaussian fields, yielding the best fit model parameters θ as

$\label{eq:smica:thetahat} \hat{θ} = \rm{arg \, min}_θ \sum_\ell (2\ell + 1) ( \mathbf{\hat{R}}_\ell \mathbf{R}_\ell (θ)^{-1} \, +\, log \, det \, \mathbf{R}_\ell (θ)).$

SMICA models the data is a superposition of , noise and foregrounds. The latter are not parametrically modelled; instead, we represent the total foreground emission by $d$ templates with arbitrary frequency spectra, angular spectra and correlations:

$\label{eq:smica:Rmodel} \mathbf{R}_\ell (θ) = \mathbf{aa}^† \, C_\ell \, + \, \mathbf{A P}_\ell \mathbf{A}^† \, + \, \mathbf{N}_\ell$

where $C_\ell$ is the angular power spectrum of the , $\mathbf{A}$ is a $N_{chan} ×d$ matrix, $\mathbf{P}_\ell$ is a positive $d×d$ matrix, and $\mathbf{N}_\ell$ is a diagonal matrix representing the noise power spectrum. The parameter vector $θ$ contains all or part of the quantities in the above equation.

The above equations summarize the founding principles of SMICA; its actual operation depends on a choice for the spectral model $\mathbf{R}_\ell (θ)$ and on several implementation-specific details.

The actual implementation of SMICA includes the following steps:

Inputs
All nine Planck frequency channels from 30 to 857 GHz, harmonically transformed up to $\ell = 4000$.
Fit
In practice, the SMICA fit, i.e. the minimization of the fourth equation, is conducted in three successive steps: We first estimate the spectral law by fitting all model parameters over a clean fraction of sky in the range $100 ≤ \ell ≤ 680$ and retaining the best fit value for vector $\mathbf{a}$. In the second step, we estimate the foreground emissivity by fixing a to its value from the previous step and fitting all the other parameters over a large fraction of sky in the range $4 ≤ \ell ≤ 150$ and retaining the best fit values for the matrix $\mathbf{A}$. In the last step, we fit all power spectrum parameters; that is, we fix $\mathbf{a}$ and $\mathbf{A}$ to their previously found values and fit for each $C_\ell$ and $\mathbf{P}_\ell$ at each $\ell$.
Beams
The discussion thus far assumes that all input maps have the same resolution and effective beam. Since the observed maps actually vary in resolution, we process the input maps in the following way. To the $i$-th input map with effective beam $b_i(\ell)$ and sampled on an grid with $N^i_{side}$, the sky multipole $s_{\ell m}$ actually contributes $s_{\ell m}a_i b_i(\ell) p_i(\ell)$, where $p_i(\ell)$ is the pixel window function for the grid at $N^i_{side}$. Seeking a final map at 5-arcmin resolution, the highest resolution of Planck, we work with input spherical harmonics re-beamed to 5 arcmins, $\mathbf{\tilde{x}}_{\ell m}$; that is, SMICA operates on vectors with entries $x ̃^i_{\ell m} = x^i_{\ell m} b_5(\ell) / b_i(\ell) / p_i(\ell)$, where $b_5(\ell)$ is a 5 arcmin Gaussian beam function. By construction, SMICA then produces an map with an effective Gaussian beam of 5 arcmin (without the pixel window function).
Pre-processing
We start by fitting point sources with SNR > 5 in the PCCS catalogue in each input map. If the fit is successful, the fitted point source is removed from the map; otherwise it is masked and the hole in-painted. This is done at all frequencies but 545 and 857 GHz, where all point sources with SNR > 7.5 are masked and in-painted.
In practice, SMICA uses a small Galactic mask leaving 97% of the sky. We deliver a full-sky map in which the masked pixels (Galactic and point-source) are replaced by a constrained Gaussian realization.
Binning
In our implementation, we use binned spectra.
High $\ell$
Since there is little point trying to model the spectral covariance at high multipoles, because the sample estimate is sufficient, SMICA implements a simple harmonic ILC at $\ell \gt 1500$; that is, it applies the filter (second equation) with $\mathbf{R}_\ell = \mathbf{\hat{R}}_\ell$.

Viewed as a filter, SMICA can be summarized by the weights $\mathbf{w}_\ell$ applied to each input map as a function of multipole. In this sense, SMICA is strictly equivalent to co-adding the input maps after convolution by specific axi-symmetric kernels directly related to the corresponding entry of $\mathbf{w}_\ell$. The SMICA weights used here are shown in figure below for input maps in units of K$_\rm{RJ}$. They show, in particular, the (expected) progressive attenuation of the lowest resolution channels with increasing multipole.

Weights $w_\ell$ given by SMICA to the input maps, after they are re-beamed to 5 arcmin and expressed in K$_\rm{RJ}$, as a function of multipole.

### Commander-Ruler

The Commander-Ruler (C-R) approach implements Bayesian component separation in pixel space, fitting a parametric model to the data by sampling the posterior distribution for the model parameters. For computational reasons, the fit is performed in a two-step procedure: First, both foreground amplitudes and spectral parameters are found at low-resolution using MCMC/Gibbs sampling algorithms (Jewell et al. 2004; Wandelt et al. 2004; Eriksen et al. 2004, 2007, 2008). Second, the amplitudes are recalculated at high resolution by solving the generalized least squares system (GLSS) per pixel with the spectral parameters fixed to their values from the low-resolution run. For the -oriented analysis presented in this paper, we use only the seven lowest Planck frequencies, i.e., from 30 to 353 GHz. We first downgrade each frequency map from its native angular resolution to a common resolution of 40 arcminutes and re-pixelize at N$_\rm{side}$ = 256. Second, we set the monopoles and dipoles for each frequency band using a method that locally conserves spectral indices (Wehus et al. 2013, in preparation). We approximate the effective instrumental noise as white with an RMS per pixel given by the Planck scanning pattern and an amplitude calibrated by smoothing simulations of the instrumental noise including correlations to the same resolution. For the high-resolution analysis, the important pre-processing step is the upgrading of the effective low-resolution mixing matrices to full Planck resolution: this is done by repixelizing from N$_\rm{side}$ = 256 to 2048 in harmonic space, ensuring that potential pixelization effects from the low-resolution map do not introduce sharp boundaries in the high-resolution map.

## References

1. Planck 2013 results. XI. Component separation, Planck Collaboration, 2014, A&A, 571, A11
2. Planck 2015 results. XI. Diffuse component separation: CMB maps, Planck Collaboration, 2016, A&A, 594, A9.
3. Component separation methods for the PLANCK mission, S. M. Leach, J.-F. Cardoso, C. Baccigalupi, R. B. Barreiro, M. Betoule, J. Bobin, A. Bonaldi, J. Delabrouille, G. de Zotti, C. Dickinson, H. K. Eriksen, J. González-Nuevo, F. K. Hansen, D. Herranz, M. Le Jeune, M. López-Caniego, E. Martínez-González, M. Massardi, J.-B. Melin, M.-A. Miville-Deschênes, G. Patanchon, S. Prunet, S. Ricciardi, E. Salerno, J. L. Sanz, J.-L. Starck, F. Stivoli, V. Stolyarov, R. Stompor, P. Vielva, A&A, 491, 597-615, (2008).
4. Multiresolution internal template cleaning: an application to the Wilkinson Microwave Anisotropy Probe 7-yr polarization data, R. Fernández-Cobos, P. Vielva, R. B. Barreiro, E. Martínez-González, MNRAS, 420, 2162-2169, (2012).