Difference between revisions of "Astrophysical component separation"

From Planck PLA Wiki
Jump to: navigation, search
 
(15 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 
==CMB and foreground separation==
 
==CMB and foreground separation==
 +
 +
See the Component Separation paper {{PlanckPapers|planck2013-p06}} for details.
 +
  
 
===NILC===
 
===NILC===
 +
 +
NILC is a linear method for combining the input channels. It implements an ILC with weighting coefficients varying over the sky and over the multipole range up to <math>\ell=3200</math> 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.
  
 
===SEVEM===
 
===SEVEM===
Line 16: Line 22:
  
 
===SMICA===
 
===SMICA===
 +
 +
 +
A linear method, SMICA reconstructs a CMB map as a linear combination
 +
in the harmonic domain of <math>N_{chan}</math> input frequency maps
 +
with weights that depend on multipole <math>\ell</math>. Given the
 +
<math>N_{chan} × 1</math> vector <math>\mathbf{x}_{\ell m}</math> of
 +
spherical harmonic coefficients for the input maps, it computes
 +
coefficients <math>s_{\ell m}</math> for the CMB map as
 +
 +
: <math>\label{eq:smica:shat}
 +
\hat{s}_{\ell  m} = \mathbf{w}^†_\ell  \mathbf{x}_{\ell  m}</math>
 +
 +
where the <math>N_{chan} × 1</math> vector <math>\mathbf{w}_\ell
 +
</math> which contains the multipole-dependent weights is built to
 +
offer unit gain to the CMB with minimum variance. This is achieved
 +
with
 +
 +
: <math>\label{eq:smica:w}
 +
\mathbf{w}_\ell  = \frac{\mathbf{R}_\ell ^{-1} \mathbf{a}}{\mathbf{a}^† \mathbf{R}_\ell ^{-1} \mathbf{a}} </math>
 +
 +
where vector <math>\mathbf{a}</math> is the emission spectrum of the
 +
CMB evaluated at each channel (allowing for possible inter-channel
 +
recalibration factors) and <math> \mathbf{R}_\ell </math> is the
 +
<math>N_{chan} × N_{chan}</math> spectral covariance matrix of
 +
<math>\mathbf{x}_{\ell m}</math>. Taking <math>\mathbf{R}_\ell </math>
 +
in Eq. \ref{eq:smica:w} to be the sample spectral covariance matrix
 +
<math>\mathbf{\hat{R}}_\ell </math> of the observations:
 +
 +
: <math>\label{eq:smica:Rhat}
 +
\mathbf{\hat{R}}_\ell  = \frac{1}{2 \ell  + 1} \sum_m \mathbf{x}_{ \ell  m} \mathbf{x}_{\ell  m}^†</math>
 +
 +
would implement a simple harmonic-domain ILC. This is not what SMICA
 +
does. As discussed below, we instead use a model <math>\mathbf{R}_\ell
 +
(θ)</math> and determine the covariance matrix to be used in
 +
Eq. \ref{eq:smica:w} by fitting <math>\mathbf{R}_\ell (θ)</math> to
 +
<math>\mathbf{\hat{R}}_\ell </math>. This is done in the maximum
 +
likelihood sense for stationary Gaussian fields, yielding the best fit
 +
model parameters θ as
 +
 +
: <math>\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 (θ)).</math>
 +
 +
 +
SMICA models the data is a superposition of CMB, noise and
 +
foregrounds. The latter are not parametrically modelled; instead, we
 +
represent the total foreground emission by <math>d</math> templates
 +
with arbitrary frequency spectra, angular spectra and correlations:
 +
 +
: <math> \label{eq:smica:Rmodel}
 +
\mathbf{R}_\ell (θ) = \mathbf{aa}^† \, C_\ell  \, + \, \mathbf{A P}_\ell  \mathbf{A}^† \, + \, \mathbf{N}_\ell
 +
</math>
 +
 +
where <math>C_\ell </math> is the angular power spectrum of the CMB,
 +
<math>\mathbf{A}</math> is a <math>N_{chan} ×d</math> matrix,
 +
<math>\mathbf{P}_\ell </math> is a positive <math>d×d</math> matrix,
 +
and <math>\mathbf{N}_\ell </math> is a diagonal matrix representing
 +
the noise power spectrum. The parameter vector <math>θ</math> contains
 +
all or part of the quantities in Eq. (5).
 +
 +
 +
The above equations summarize the founding principles of SMICA; its
 +
actual operation depends on a choice for the spectral model
 +
<math>\mathbf{R}_\ell (θ)</math> 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  <math>\ell = 4000 </math>.
 +
; Fit
 +
: In practice, the SMICA fit,i.e.,the minimization of Eq. (4), is conducted in three successive steps: We first estimate the CMB spectral law by fitting all model parameters over a clean fraction of sky in the range <math> 100 ≤ \ell ≤ 680</math>  and retaining the best fit value for vector <math> \mathbf{a}</math>. 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 <math> 4 ≤ \ell  ≤ 150</math>  and retaining the best fit values for the matrix <math> \mathbf{A}</math>. In the last step, we fit all power spectrum parameters; that is, we fix <math>\mathbf{a}</math> and <math>\mathbf{A}</math> to their previously found values and fit for each <math> C_\ell </math>  and  <math>\mathbf{P}_\ell </math>  at each <math>\ell</math>.
 +
;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 <math>i</math>-th input map with effective beam <math>b_i(\ell)</math> and sampled on an HEALPix grid with <math>N^i_{side}</math>, the CMB sky multipole <math>s_{\ell m}</math> actually contributes <math>s_{\ell m}a_i b_i(\ell) p_i(\ell)</math>, where <math>p_i(\ell)</math> is the pixel window function for the grid at <math>N^i_{side}</math>. Seeking a final CMB map at 5-arcmin resolution, the highest resolution of Planck, we work with input spherical harmonics re-beamed to 5 arcmins, <math>\mathbf{\tilde{x}}_{\ell m} </math>; that is, SMICA operates on vectors with entries <math>x ̃^i_{\ell m} = x^i_{\ell m} b_5(\ell) / b_i(\ell) / p_i(\ell)</math>, where <math>b_5(\ell)</math> is a 5 arcmin Gaussian beam function. By construction, SMICA then produces an CMB 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. 
 +
; Masking and in-painting
 +
: In practice, SMICA uses a small Galactic mask leaving 97% of the sky. We deliver a full-sky CMB 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 <math>\ell</math>
 +
: 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 <math>\ell > 1500</math>; that is, it applies the filter (Eq. 2) with <math>\mathbf{R}_\ell = \mathbf{\hat{R}}_\ell</math>.
 +
 +
Viewed as a filter, SMICA can be summarized by the weights <math>\mathbf{w}_\ell</math> 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 <math>\mathbf{w}_\ell</math>. The SMICA weights used here are shown in figure below for input maps in units of K<math>_\rm{RJ}</math>. They show, in particular, the (expected) progressive attenuation of the lowest resolution channels with increasing multipole.
 +
 +
[[File:smica.jpg|thumb|center|600px|'''Weights <math>w_\ell</math> given by SMICA to the input maps, after they are re-beamed to 5 arcmin and expressed in K<math>_\rm{RJ}</math>, as a function of multipole.''']]
  
 
===Commander-Ruler===
 
===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 the their values from the low-resolution run.
 +
For the CMB-oriented analysis presented in this paper, we only use 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 HEALPix N<math>_\rm{side}</math> = 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<math>_\rm{side}</math> = 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.
 +
 +
<!--
 +
TBW.
 +
-->
 +
 +
== References ==
 +
 +
<References />
 +
  
==CO Maps==
 
  
[[Category:Data processing|0062]]
+
[[Category:HFI/LFI joint data processing|002]]

Latest revision as of 16:38, 23 July 2014

CMB and foreground separation[edit]

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


NILC[edit]

NILC is a linear method for combining the input channels. It implements an ILC with weighting coefficients varying over the sky and over the multipole range up to [math]\ell=3200[/math] 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.

SEVEM[edit]

The aim of Sevem is to produce clean CMB 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 and to WMAP polarisation data Fernandez-Cobos et al., 2012. 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, we have cleaned the 100, 143 GHz and 217 GHz maps using four templates constructed as the difference of the following Planck channels (smoothed to a common resolution): (30-44) GHz, (44-70) GHz, (545-353) GHz and (857-545)GHz.

The templates are constructed by subtracting two neighbouring Planck frequency channel maps, after first smoothing them to a common resolution to ensure that the CMB signal is properly removed. A linear combination of the templates is then subtracted from the Planck sky map at the frequency to be cleaned, in order to produce the clean CMB. 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).

An additional level of flexibility can also be considered: the linear coefficients can be the same for all the sky, or several regions with different sets of coefficients can be considered. The regions are then combined in a smooth way, by weighting the pixels at the boundaries, to avoid discontinuities in the clean maps. In order to take into account the different spectral behaviour of the foregrounds at low and high galactic latitudes, we have chosen to use two regions: the region with the 3 per cent brightest Galactic emission, and the region with the remaining 97 per cent of the sky.

The final CMB 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 FWHM =5 arcminutes.

SMICA[edit]

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

[math]\label{eq:smica:shat} \hat{s}_{\ell m} = \mathbf{w}^†_\ell \mathbf{x}_{\ell m}[/math]

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

[math]\label{eq:smica:w} \mathbf{w}_\ell = \frac{\mathbf{R}_\ell ^{-1} \mathbf{a}}{\mathbf{a}^† \mathbf{R}_\ell ^{-1} \mathbf{a}} [/math]

where vector [math]\mathbf{a}[/math] is the emission spectrum of the CMB evaluated at each channel (allowing for possible inter-channel recalibration factors) and [math] \mathbf{R}_\ell [/math] is the [math]N_{chan} × N_{chan}[/math] spectral covariance matrix of [math]\mathbf{x}_{\ell m}[/math]. Taking [math]\mathbf{R}_\ell [/math] in Eq. \ref{eq:smica:w} to be the sample spectral covariance matrix [math]\mathbf{\hat{R}}_\ell [/math] of the observations:

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

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

[math]\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 (θ)).[/math]


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

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

where [math]C_\ell [/math] is the angular power spectrum of the CMB, [math]\mathbf{A}[/math] is a [math]N_{chan} ×d[/math] matrix, [math]\mathbf{P}_\ell [/math] is a positive [math]d×d[/math] matrix, and [math]\mathbf{N}_\ell [/math] is a diagonal matrix representing the noise power spectrum. The parameter vector [math]θ[/math] contains all or part of the quantities in Eq. (5).


The above equations summarize the founding principles of SMICA; its actual operation depends on a choice for the spectral model [math]\mathbf{R}_\ell (θ)[/math] 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 [math]\ell = 4000 [/math].
Fit
In practice, the SMICA fit,i.e.,the minimization of Eq. (4), is conducted in three successive steps: We first estimate the CMB spectral law by fitting all model parameters over a clean fraction of sky in the range [math] 100 ≤ \ell ≤ 680[/math] and retaining the best fit value for vector [math] \mathbf{a}[/math]. 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 [math] 4 ≤ \ell ≤ 150[/math] and retaining the best fit values for the matrix [math] \mathbf{A}[/math]. In the last step, we fit all power spectrum parameters; that is, we fix [math]\mathbf{a}[/math] and [math]\mathbf{A}[/math] to their previously found values and fit for each [math] C_\ell [/math] and [math]\mathbf{P}_\ell [/math] at each [math]\ell[/math].
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 [math]i[/math]-th input map with effective beam [math]b_i(\ell)[/math] and sampled on an HEALPix grid with [math]N^i_{side}[/math], the CMB sky multipole [math]s_{\ell m}[/math] actually contributes [math]s_{\ell m}a_i b_i(\ell) p_i(\ell)[/math], where [math]p_i(\ell)[/math] is the pixel window function for the grid at [math]N^i_{side}[/math]. Seeking a final CMB map at 5-arcmin resolution, the highest resolution of Planck, we work with input spherical harmonics re-beamed to 5 arcmins, [math]\mathbf{\tilde{x}}_{\ell m} [/math]; that is, SMICA operates on vectors with entries [math]x ̃^i_{\ell m} = x^i_{\ell m} b_5(\ell) / b_i(\ell) / p_i(\ell)[/math], where [math]b_5(\ell)[/math] is a 5 arcmin Gaussian beam function. By construction, SMICA then produces an CMB 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.
Masking and in-painting
In practice, SMICA uses a small Galactic mask leaving 97% of the sky. We deliver a full-sky CMB 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 [math]\ell[/math]
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 [math]\ell \gt 1500[/math]; that is, it applies the filter (Eq. 2) with [math]\mathbf{R}_\ell = \mathbf{\hat{R}}_\ell[/math].

Viewed as a filter, SMICA can be summarized by the weights [math]\mathbf{w}_\ell[/math] 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 [math]\mathbf{w}_\ell[/math]. The SMICA weights used here are shown in figure below for input maps in units of K[math]_\rm{RJ}[/math]. They show, in particular, the (expected) progressive attenuation of the lowest resolution channels with increasing multipole.

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

Commander-Ruler[edit]

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 the their values from the low-resolution run. For the CMB-oriented analysis presented in this paper, we only use 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 HEALPix N[math]_\rm{side}[/math] = 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[math]_\rm{side}[/math] = 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[edit]

  1. Planck 2013 results: Component separation, Planck Collaboration XII, A&A, in press, (2014).

Cosmic Microwave background

Full-Width-at-Half-Maximum

(Hierarchical Equal Area isoLatitude Pixelation of a sphere, <ref name="Template:Gorski2005">HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere, K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, M. Bartelmann, ApJ, 622, 759-771, (2005).