SILC: a new Planck internal linear combination CMB temperature map using directional wavelets

Abstract

We present new clean maps of the cosmic microwave background (CMB) temperature anisotropies (as measured by Planck) constructed with a novel internal linear combination (ILC) algorithm using directional, scale-discretized wavelets – scale-discretized, directional wavelet ILC or Scale-discretised, directional wavelet Internal Linear Combination (SILC). Directional wavelets, when convolved with signals on the sphere, can separate the anisotropic filamentary structures which are characteristic of both the CMB and foregrounds. Extending previous component separation methods, which use the frequency, spatial and harmonic signatures of foregrounds to separate them from the cosmological background signal, SILC can additionally use morphological information in the foregrounds and CMB to better localize the cleaning algorithm. We test the method on Planck data and simulations, demonstrating consistency with existing component separation algorithms, and discuss how to optimize the use of morphological information by varying the number of directional wavelets as a function of spatial scale. We find that combining the use of directional and axisymmetric wavelets depending on scale could yield higher quality CMB temperature maps. Our results set the stage for the application of SILC to polarization anisotropies through an extension to spin wavelets.

methods: data analysis, cosmic background radiation

1 INTRODUCTION

Accurate measurements of the cosmic microwave background (CMB) arguably form the bedrock of modern precision cosmology. In particular, the full-sky multifrequency CMB maps provided by three generations of satellite experiments – COBE (Mather et al. 1990; Boggess et al. 1992), Wilkinson Microwave Anisotropy Probe (WMAP; Bennett et al. 2003a) and Planck (Planck Collaboration I 2011) – represent milestones in our understanding of the cosmological model. However, to obtain a full-sky map of the CMB requires removing instrumental noise and signals due to astrophysical foregrounds (primarily in the Milky Way). Full-sky foreground-cleaned CMB maps are used for a wide variety of scientific purposes (see e.g. Planck Collaboration XVI 2015; Planck Collaboration XVII 2015).

There are numerous methods to perform foreground component separation. They broadly fall into two categories: blind methods which make minimal physical assumptions about the contributing signals and the so-called mixing matrix (which quantifies the strength of different components at different frequencies) and non-blind methods which are based on a physical modelling of the sky components. Examples of non-blind methods include the maximum entropy method (Hobson et al. 1998) and the parametric Bayesian CMB Gibbs sampler Commander (Eriksen et al. 2006, 2008). Correlated component analysis (Bedini et al. 2005) is a semiblind method that estimates the mixing matrix based on second-order statistics. Spectral Estimation via Expectation Maximization (SEVEM; Martínez-González et al. 2003; Leach et al. 2008; Fernández-Cobos et al. 2012) is a template-fitting technique. Examples of so-called blind source separation include the sparsity based method local-generalized morphological component analysis (Bobin et al. 2008, 2013) and the spectral matching independent component analysis (SMICA; Cardoso et al. 2008), although the latter work does discuss how the choice of component model affects the blindness of this method. Of particular interest to this work is another blind method, the internal linear combination (ILC), most recently implemented by the needlet ILC (NILC; Delabrouille et al. 2009). In its component separation analysis, the Planck Collaboration used Commander, NILC, SEVEM and SMICA (Planck Collaboration IX 2015). See, e.g. Delabrouille et al. (2009) and Bobin et al. (2013) for reviews of CMB component separation methods.

The ILC computes a weighted sum of CMB maps as measured at multiple frequencies. These weights are constrained to sum to unity at each point in the map, ensuring that the CMB signal is conserved, assuming that it is equal at each frequency. Under this constraint, the weights are calculated by minimizing the empirical variance of the ILC map, which in turn minimizes the variance of the error in CMB reconstruction (assuming the CMB and foregrounds and the CMB and noise are, respectively, uncorrelated). The variances we minimize are empirical in that they are calculated using the data themselves. In this way, the weights are calculated to remove foreground and noise, revealing the underlying primordial CMB anisotropies. The ILC method was originally used by the WMAP Collaboration (Bennett et al. 2003b) and then extended by Eriksen et al. (2004) through an analytical calculation of the weights. One limitation of the original ILC approach is the extent of localization of the weights. The initial versions calculated different weights in separate parts of the sky (e.g. Bennett et al. 2003b split the Galactic region into 11 parts). In order to further remove local contamination, the weights can be allowed to vary across the sky and also at different spatial scales. Tegmark, de Oliveira-Costa & Hamilton (2003) made an ILC map allowing the weights to vary at each multipole, as well as within different regions of the sky. A direct extension of this work is to make use of both spatial and frequency information simultaneously using wavelets. The weights are then defined across wavelet scales and within wavelet coefficient maps on the sky.

Wavelets are functions that are localized in both real and frequency space. To analyse full-sky CMB maps, wavelets defined on the surface of a sphere are required. A number of wavelet frameworks on the sphere have been developed recently (Antoine & Vandergheynst 1998, 1999; Barreiro et al. 2000; Wiaux, Jacques & Vandergheynst 2005; McEwen, Hobson & Lasenby 2006; Narcowich, Petrushev & Ward 2006; Sanz et al. 2006; Starck et al. 2006; Geller et al. 2008; McEwen & Scaife 2008; Marinucci et al. 2008; Wiaux et al. 2008; Baldi et al. 2009; Starck, Moudden & Bobin 2009; Geller & Marinucci 2010, 2011; McEwen, Wiaux & Eyers 2011; Leistedt et al. 2013, 2015; McEwen, Vandergheynst & Wiaux 2013; McEwen et al. 2014, 2015b; McEwen, Durastanti & Wiaux 2015a). In particular, needlets (Narcowich et al. 2006; Marinucci et al. 2008; Baldi et al. 2009) have been used in the latest generation of ILC methods. Needlets are a set of axisymmetric kernels defined on the surface of a sphere. Each member of the set has compact support in harmonic space over different multipole ranges. When each needlet is convolved with a signal defined on the sphere, the resulting signal (i.e. needlet coefficients) also has compact harmonic support. NILC (Delabrouille et al. 2009) computes its weights by considering needlet scales separately (harmonic localization) and then different parts of each needlet coefficient map separately (spatial localization). The needlets are constructed in such a way that the original signal can be recovered from its needlet coefficients with no loss of information (in practice, small losses can be introduced by approximate spherical harmonic transforms). NILC has been very successful at forming clean full-sky CMB maps, which contain very little residual foreground and noise contamination.

In this work, we introduce the scale-discretized, directional wavelet ILC or Scale-discretised, directional wavelet Internal Linear Combination (SILC), which extends the wavelet ILC framework by localizing the calculation of ILC weights in an additional domain. We use wavelets which are not only harmonically localized but also directional (Wiaux et al. 2008; McEwen et al. 2013; McEwen et al. 2015b). Unlike needlets, which are axisymmetric on the sphere, directional wavelets are non-axisymmetric, i.e. the kernels are ‘squeezed’. This means that for one wavelet scale, one axisymmetric kernel is replaced by a number of complementary directional kernels, each with a different orientation. When these directional wavelets are convolved with a signal on the sphere, (within each scale) different orientations of signal structure are separated. This directional localization allows the ILC weights to be additionally fine-tuned to better remove foreground and noise, in particular for signals with filamentary structure. Furthermore, directional wavelets exhibit exact reconstruction, allowing them to be embedded in an ILC such that no signal is lost.

SILC is being developed with the goal of analysing CMB polarization components through an extension to spin, directional wavelets (McEwen et al. 2014, 2015b; Leistedt et al. 2015), which are expected to be well suited to localizing the complex filamentary morphologies of polarized foregrounds. As a precursor step, in this work we test SILC on the scalar temperature field in order to demonstrate the quantitative consistency of its foreground cleaning performance compared with existing component separation methods, and to identify possible optimizations for the extension to spin fields.

Directional wavelets are explained briefly in Section 2. In Section 3, the SILC algorithm is explained in detail. Various sources of error in the method are considered in Section 4. In Section 5, we compare our method to previous component-separation methods. The application to Planck simulations (Section 6) is followed by application to Planck data (Section 7). We discuss the results and error estimation based on the data in Section 8 and conclude in Section 9.

2 DIRECTIONAL WAVELETS

Directional, scale-discretized wavelets on the sphere that support exact reconstruction have been developed in Wiaux et al. (2008), McEwen et al. (2013) and McEwen et al. (2015b), while their localization properties have been studied in McEwen et al. (2015a). Fig. 1 shows an example of the spatial localization of directional wavelets. Larger wavelet scales have larger kernels, and when these are convolved with signals defined on the sphere (such the CMB and astrophysical foregrounds), signal structure with the same scale and orientation as the wavelet is isolated. The kernels in Fig. 1 are shown for a single direction and (for complete reconstruction) would be complemented by two more sets of kernels of the same sizes but rotated to different orientations. Fig. 2 shows an example of the harmonic localization of directional wavelets (for the wavelets used in this work). The harmonic supports of the wavelets overlap, with each wavelet covering a finite set of multipoles. Fig. 3 shows an example of directional wavelet decomposition as applied to the CMB. Although the CMB anisotropies are statistically Gaussian, the CMB spots on the sky demonstrate anisotropy as a function of scale (Bond & Efstathiou 1987). When the CMB is convolved with directional wavelets, structure of different orientations is separated. This further supports the use of directional wavelets in CMB analysis: both filamentary foreground structure and the CMB itself are better localized. This particularly applies in the case of polarization, as will be discussed in Section 9 when we consider extensions to our method. For a mathematical description of directional wavelets, see Section 3.4.

Figure 1.

The spatial localization on the sphere of directional, scale-discretized wavelets. Each sub-plot shows a representation of a directional wavelet kernel at different scales, where red, raised parts show positive wavelet response and blue, depressed parts show negative wavelet response. From left to right, top to bottom: wavelet scale index j decreases. The number of directions per wavelet scale N = 3. Therefore, for complete reconstruction at each scale, the above wavelets would be complemented by two more wavelets of the same size but of a different orientation on the sphere. This figure is adapted from McEwen et al. (2013).

Figure 2.

The harmonic response of the directional wavelets used in this work, where j specifies the wavelet scale. Increasing j corresponds to a smaller wavelet kernel and so a multipole range on smaller scales (i.e. larger multipoles ℓ). The largest wavelet scale (Scal.) is the scaling function (Section 3.4). The two smallest wavelets are harmonically truncated at ℓ = 3600 but are smoothly tapered to zero from ℓ = 3400 to ℓ = 3600 (the two dotted lines) by the beam tapering discussed in Section 3.3. The band-limits of the above wavelets are given in Table 1.

Figure 3.

The CMB (top map) decomposed into directional wavelet coefficient maps (bottom section). The wavelet kernels are shown (middle section), where red indicates positive response and blue indicates negative. In the full analysis, we also include smaller wavelets than we show above.

In the spherical harmonic transforms used in the computation of directional wavelet coefficient maps, we adopt the sampling scheme on the sphere of McEwen & Wiaux (2011, hereafter MW sampling), rather than, e.g. healpix sampling (Górski et al. 2005), although in principle healpix could be used if desired. The corresponding sampling theorem of McEwen & Wiaux (2011) shows that the MW sampling scheme requires fewer samples for a band-limited signal than any other sampling theorem. Additionally, the use of a separation of variables and fast Fourier transforms (FFTs) yields a numerically efficient algorithm. In particular, our spherical harmonic transforms are theoretically exact, unlike healpix. This allows one to manipulate signals with the minimal number of samples and to perform the numerous spherical harmonic transforms involved in the ILC algorithm without any loss of information (other than that due to the finite representation of floating point numbers). Our final map products, however, are provided in healpix format. Finally, MW sampling of spin signals on the sphere requires no additional computational complexity and this will be vital in the extension of our method to polarization E and B modes (Section 9). Further details on MW sampling are given in Section 3.4.

3 METHOD

We start by outlining the SILC algorithm. The steps are explained in more detail in the subsequent subsections (Sections 3.1–3.6). We discuss our numerical implementation in Section 3.7.

1. The raw input data are multifrequency full-sky maps of CMB temperature fluctuations. These maps use the healpix format. (See Section 3.1.) The model we employ for the raw data is explained in Section 3.2.

2. The maps are ‘pre-processed’ by inpainting in a small point source mask (see Section 3.6).

3. The input maps are converted to thermodynamic (CMB) temperature (if necessary). For Planck temperature data, the 545 and 857 GHz maps are converted from spectral flux density per unit solid angle (MJy Sr⁻¹) to CMB temperature (K_CMB) by the unit conversions given in the Planck 2015 Release Explanatory Supplement.¹

4. The maps are each convolved to have the same effective beam (see Section 3.3).

5. Each input map is converted into a set of wavelet coefficient maps. This separates both the scale and orientation of structure within each map. These wavelet coefficient maps use MW sampling (McEwen & Wiaux 2011). (See Section 3.4.)

6. The ILC method is then applied separately to each wavelet scale and orientation. For each scale and orientation, the multifrequency wavelet coefficient maps are weighted and added to form a single wavelet coefficient map that contains mainly CMB signal, as well as some residual foreground and noise. These weights are allowed to vary from wavelet coefficient to wavelet coefficient. The calculation of these weights is explained in Section 3.5.

7. The final ILC wavelet coefficient maps are synthesized to form the final product: a full-sky map of CMB temperature fluctuations (with some residual foreground and noise). The final map uses the healpix format. (See Section 3.4.)

8. The final map is inpainted in a small point source mask (see Section 3.6).

3.1 Input data

Our main CMB temperature map products use full-mission 2015 release Planck temperature maps as their input.² All nine frequency channels are used. At 70 GHz, we use the higher resolution version at N_side = 2048. We also use the full-mission full focal plane 8 (FFP8) simulations (Planck Collaboration XII 2015) without bandpass mismatch.³

3.2 Data model

Each full-sky temperature map can be modelled (e.g. Basak & Delabrouille 2012) as

\begin{equation} T^{\mathrm{OBS},c}(\hat{\boldsymbol {n}}) = \int _{\hat{\boldsymbol {n}}^{\prime }} \mathrm{d}\hat{\boldsymbol {n}}^{\prime } B^c(\hat{\boldsymbol {n}},\hat{\boldsymbol {n}}^{\prime }) T^{\mathrm{SIG},c}(\hat{\boldsymbol {n}}^{\prime }) + T^{\mathrm{N},c}(\hat{\boldsymbol {n}}), \end{equation}

(1)

where the signal component can further be decomposed as

\begin{equation} T^{\mathrm{SIG},c}(\hat{\boldsymbol {n}}) = a^c T^{\mathrm{CMB}}(\hat{\boldsymbol {n}}) + T^{\mathrm{FG},c}(\hat{\boldsymbol {n}}). \end{equation}

(2)

|$T^{\mathrm{CMB}}(\hat{\boldsymbol {n}})$| is the CMB component at a point on the sky |$\hat{\boldsymbol {n}}$|⁠. |$T^{\mathrm{FG},c}(\hat{\boldsymbol {n}})$| and |$T^{\mathrm{N},c}(\hat{\boldsymbol {n}})$| are, respectively, the foreground and detector noise components for frequency channel c. a^c is the calibration coefficient for the CMB for each channel. The overall signal component is smoothed by a beam function |$B^c(\hat{\boldsymbol {n}},\hat{\boldsymbol {n}}^{\prime })$| due to the finite resolution of the observations. However, the noise component is not smoothed by the beam. Here, we assume the beam to be circularly symmetric. Therefore, the beam can be represented as a sum over Legendre polynomials,

\begin{equation} B^c(\hat{\boldsymbol {n}},\hat{\boldsymbol {n}}^{\prime }) = \sum _{\ell =0}^\infty \frac{2\ell +1}{4\pi } B_\ell ^c P_\ell (\hat{\boldsymbol {n}}.\hat{\boldsymbol {n}}^{\prime }). \end{equation}

(3)

We can recast equation (1) in the spherical harmonic representation as

\begin{equation} a_{\ell m}^{\mathrm{OBS},c} = a^c B_\ell ^c a_{\ell m}^\mathrm{CMB} + B_\ell ^c a_{\ell m}^{\mathrm{FG},c} + a_{\ell m}^{\mathrm{N},c}, \end{equation}

(4)

where a_ℓm are the coefficients of spherical harmonics |$Y_{\ell m}(\hat{\boldsymbol {n}})$|⁠.

3.3 Beam convolution

Equation (4) shows that each frequency channel c has a different beam transfer function |$B_\ell ^c$|⁠. To replace each beam with a channel-independent resolution, we perform a deconvolution/convolution procedure to give spherical harmonic coefficients

\begin{equation} a_{\ell m}^c = \frac{B_\ell ^\mathrm{EFF}}{B_\ell ^c} a_{\ell m}^{\mathrm{OBS},c} \,, \end{equation}

(5)

where |$B_\ell ^\mathrm{EFF}$| is the final (effective) beam transfer function of our map products. For Planck data, we use a Gaussian beam with a full width at half-maximum (FWHM) of 5 arcmin as our input beam. We taper this beam to zero from ℓ = 3400 to ℓ = 3600 using a Fermi function. This suppresses any small-scale power aliasing due to having harmonically truncated wavelets in this multipole range. Convolving with beam transfer functions ignores the non-axisymmetric component of the beams; these will remain in the input maps but are assumed to be small.

This deconvolution/convolution procedure does not correctly handle the noise component of our input maps. Equation (5) can be expanded (using equation 4) as

\begin{equation} a_{\ell m}^c = B_\ell ^\mathrm{EFF} (a^c a_{\ell m}^\mathrm{CMB} + a_{\ell m}^{\mathrm{FG},c}) + \frac{B_\ell ^\mathrm{EFF}}{B_\ell ^c} a_{\ell m}^{\mathrm{N},c}, \end{equation}

(6)

where a^c are the CMB calibration coefficients (not to be confused with the inverse spherical harmonic transform of harmonic coefficients |$a_{lm}^c$|⁠). The final resolution of an ILC map is usually chosen to match the best resolution of the input maps. Therefore, for all but the highest resolution channel and for all ℓ, |$B_\ell ^\mathrm{EFF} > B_\ell ^c$|⁠. This has the effect of increasing the noise contribution of the input maps, particularly at high ℓ and for low-resolution maps, where |$B_\ell ^\mathrm{EFF} \gg B_\ell ^c$|⁠. We use the Planck beam transfer functions as provided in the Reduced Instrument Model (RIMO).⁴ For the Low Frequency Instrument (LFI) beams, we use Gaussian approximations with FWHM 32.33, 27.01 and 13.25 arcmin for 30, 44 and 70 GHz, respectively. Following Planck Collaboration XII (2014), the deconvolved beams are thresholded such that the |$B_l^c$| is set to the value given in the RIMO or 0.001, whichever is larger. This prevents the last term in equation (6) from becoming so large that numerical errors occur. Although we lose accuracy in the deconvolution process, the contribution of the channels in the multipole ranges affected is highly attenuated in the ILC weights in any case.

3.4 Wavelet analysis and synthesis

The wavelet ILC method requires the decomposition of each band-limited temperature map |$T^c(\hat{\boldsymbol {n}})$| into a set of wavelet coefficient maps |$W^{\Psi ^j}$|⁠: in our case, directional, scale-discretized wavelets (Wiaux et al. 2008; McEwen et al. 2013; McEwen et al. 2015b). A general introduction was provided in Section 2 – here, we provide some technical details of the implementation. We drop the c superscript on T for the rest of this subsection since each map is analysed using the same wavelets. The wavelet coefficients are defined as the directional convolution of T with wavelets defined on the sphere |$\Psi ^j \in \mathrm{L}^2(\mathbb {S}^2)$| [specifically those shown in Fig. 2], where index j denotes the wavelet scale. Importantly, directional wavelets yield coefficients |$W^{\Psi ^j} (\hat{\boldsymbol {\rho }})$| that live on the space of three-dimensional rotations, i.e. the rotation group SO(3):

\begin{equation} W^{\Psi ^j} (\hat{\boldsymbol {\rho }}) \equiv \langle T, \mathcal {R}_{\hat{\boldsymbol {\rho }}} \Psi ^j \rangle = \int _{\mathbb {S}^2} \mathrm{d}\hat{\boldsymbol {n}} \ T(\hat{\boldsymbol {n}}) (\mathcal {R}_{\hat{\boldsymbol {\rho }}} \Psi ^j)^{{\ast }} (\hat{\boldsymbol {n}}), \end{equation}

(7)

where |$\mathrm{d}\hat{\boldsymbol {n}}$| is the usual invariant measure on the sphere, |$(\mathcal {R}_{\hat{\boldsymbol {\rho }}} \Psi ^j)^{{\ast }}$| denotes complex conjugation and the rotation operator is defined by

\begin{equation} (\mathcal {R}_{\hat{\boldsymbol {\rho }}} \Psi ^j)(\hat{\boldsymbol {n}}) \equiv \Psi ^j (\bf{R}_{\hat{\boldsymbol {\rho }}}^{-1} \hat{\boldsymbol {n}}), \end{equation}

(8)

where |$\bf{R}_{\hat{\boldsymbol {\rho }}}$| is the three-dimensional rotation matrix corresponding to |$\mathcal {R}_{\hat{\boldsymbol {\rho }}}$|⁠. In these equations, |$\hat{\boldsymbol {\rho }} = (\theta ,\phi ,\chi ) \in \mathrm{SO}(3)$| denotes the Euler angles (in the zyz convention) with colatitude θ ∈ [0, π], longitude ϕ ∈ [0, 2π) and direction χ ∈ [0, 2π).⁵ In other words, the wavelet coefficients probe directional structure in T with χ corresponding to the orientation about each point (θ, ϕ) on the sphere.

Following the directional construction of scale-discretized wavelets (Wiaux et al. 2008; McEwen et al. 2013; McEwen et al. 2015b), wavelets are defined by their spherical harmonic coefficients in factorized form:

\begin{equation} \Psi _{\ell n}^j \equiv \kappa ^j(\ell ) s_{\ell n}, \end{equation}

(9)

where κ^j(ℓ) sets the harmonic localization (Fig. 2) and s_ℓn sets the directional localization.

In the original definition of scale-discretized wavelets, the size of all harmonic kernels (setting the harmonic localization of the wavelets) is parametrized by a unique wavelet dilation parameter |$\lambda \in \mathbb {R}_{{\ast }}^{+}$|⁠, λ > 1. Similarly, the number of directions is set by a unique azimuthal band-limit N. These two parameters, respectively, characterize κ^j(ℓ) and s_ℓn for all j. In this work, we vary λ as a function of multipole in order to allow more flexible harmonic localization. We achieve this by defining different values of λ in different multipole regions and then stitching together harmonically truncated wavelets at the region boundaries. We use the values λ = 2, 1.3, 1.2 with transitions at the multipoles ℓ = 512, 2015. If at a transition multipole the harmonic peak of the larger wavelet does not equal the peak of the smaller wavelet, then a small amount of unit response is used so that the two wavelets can be continuously combined. Wavelets constructed in this manner satisfy the standard admissibility criterion required for exact reconstruction. The harmonic tiling of the resulting wavelets is shown in Fig. 2. The technical details of the construction of each kernel is described in McEwen et al. (2015a). Finally, we use a single parameter N for all scales, i.e. each wavelet is divided into the same number of directions. However, a possible extension of this work is to vary N as a function of scale j, e.g. by using curvelet kernels (Chan et al. 2015) or other directional optimizations.

In the case of a single parameter λ, the limits of the wavelet harmonic window for scale j are simply |$(\ell _\mathrm{min}^j,\ell _\mathrm{max}^j) = (\lambda ^{j-1},\lambda ^{j+1})$|⁠, with their peak response at λ^j. In our hybrid scheme, this property remains but j and λ must be adjusted in each harmonic region. The full details of our tiling are given in Table 1. When the limits of the harmonic windows of the maximum wavelet scales extend beyond the overall band-limit ℓ_max, the windows are truncated at ℓ_max. Finally, note that a scaling function W^Φ is needed to capture the very low frequency content of the signal. It is axisymmetric and the corresponding scaling coefficients therefore live on the sphere. Here, we do not give the full details of the construction of the scaling function or the factors κ^j(ℓ) and s_ℓn, since these can be straightforwardly reproduced by following previous approaches (Wiaux et al. 2008; McEwen et al. 2013; McEwen et al. 2015b) and using Table 1.

Table 1.

The harmonic band-limits |$[\ell ^j_{\mathrm{min}},\ell ^j_{\mathrm{max}}]$| of the directional wavelets used in this work. |$\ell ^j_{\mathrm{peak}}$| is the multipole at which each wavelet has its maximum response. The final column shows the number of equiangular samples per wavelet coefficient map |$N^j_\mathrm{samp}$|⁠.

Wavelet scale j	|$\ell ^j_{\mathrm{min}}$|	|$\ell ^j_{\mathrm{peak}}$|	|$\ell ^j_{\mathrm{max}}$|	|$N^j_\mathrm{samp}$|
Scal.	0	64	64	8385
0	32	64	128	33 153
1	64	128	256	131 841
2	128	256	512	525 825
3	256	512	706	998 991
4	542	705	918	1688 203
5	705	917	1193	2850 078
6	917	1192	1551	4815 856
7	1192	1550	2015	8126 496
8	1550	2015	2540	12 910 821
9	2116	2539	3048	18 589 753
10	2539	3047	3600	25 930 801
11	3047	3600	3600	25 930 801

To apply the ILC algorithm, the above continuous wavelet coefficients must be discretized. Since they live on the rotation group SO(3), we represent them using the sampling scheme of McEwen et al. (2015c), which is a generalization of the MW sampling scheme (McEwen & Wiaux 2011). Because our wavelets have well-defined band-limits, this approach allows a multiresolution scheme where each scale is pixellated with a minimal number of samples. In practice, the j-th wavelet scale has a band-limit |$\ell _\mathrm{max}^j$| and is only evaluated at locations |$(\theta _t^j, \phi _p^j, \chi _n)$| with |$t \in \lbrace 0,1,\ldots ,\ell _\mathrm{max}^j\rbrace$|⁠, |$p \in \lbrace 0,1,\ldots ,2\ell _\mathrm{max}^j\rbrace$| and n ∈ {0, 1, …, N − 1}. Although wavelet coefficients are evaluated at discrete samples only, for a band-limited signal they capture the total information content of the underlying continuous wavelet coefficient representation, probed up to harmonic band-limit |$\ell _\mathrm{max}^j$| and azimuthal band-limit N. Thanks to the sampling theory on the rotation group SO(3) of McEwen et al. (2015c). In the full ensemble of realizations, the ILC (see Section 3.5 for details) has no sensitivity to the choice of coordinate convention for directions χ. In a single realization, there will be some marginal sensitivity to this choice manifesting in the localization of the empirical covariances we use. However, this effect is sub-dominant to the choice of N, on which we concentrate our analysis.

After the ILC method (see Section 3.5) has been applied to the sets of wavelet coefficient maps, there is one final map |$W^{\Psi ^j,\mathrm{ILC}}(\hat{\boldsymbol {\rho }})$|⁠, for each wavelet scale j, living on SO(3) and including the multiple orientations χ₀, …, χ_N−1. The additional axisymmetric scaling coefficients |$W^{\Phi ,\mathrm{ILC}}(\hat{\boldsymbol {n}})$| live on the sphere. The final temperature map |$T^\mathrm{ILC}(\hat{\boldsymbol {n}})$| is synthesized by

\begin{eqnarray} T^\mathrm{ILC}(\hat{\boldsymbol {n}})& = &\int _{\mathbb {S}^2} \mathrm{d}\hat{\boldsymbol {n}}^{\prime } W^{\Phi ,\mathrm{ILC}}(\hat{\boldsymbol {n}}^{\prime }) (\mathcal {R}_{\hat{\boldsymbol {n}}^{\prime }} \Phi ) (\hat{\boldsymbol {n}}) \nonumber\\ &&+ \sum _{j=j_\mathrm{min}}^{j_\mathrm{max}} \int _{\mathrm{SO}(3)} \mathrm{d} \hat{\boldsymbol {\rho }} \ W^{\Psi ^j,\mathrm{ILC}}(\hat{\boldsymbol {\rho }}) (\mathcal {R}_{\hat{\boldsymbol {\rho }}} \Psi ^j) (\hat{\boldsymbol {n}}), \end{eqnarray}

(10)

where |$\mathrm{d}\hat{\boldsymbol {\rho }}$| is the usual invariant measure on the rotation group. This final ILC temperature map is pixellated using the healpix format from its spherical harmonic coefficients |$T^\mathrm{ILC}_{\ell m}$|⁠.

The wavelet analysis and synthesis are performed using the latest version of the s2let⁶ code (Leistedt et al. 2013; McEwen et al. 2015c), which in turn relies on the ssht⁷ (McEwen & Wiaux 2011) and so3⁸ (McEwen et al. 2015c) codes to compute spin spherical harmonics and Wigner transforms exactly and efficiently using the MW sampling scheme. Thanks to the sampling theorem, the wavelet coefficients can be transformed using Wigner transforms without any loss of information (McEwen et al. 2015c).

3.5 ILC method

Following the wavelet analysis of the input maps (see Section 3.4), there is a wavelet coefficient map |$W_{jnk}^c$| for each channel c, scale j and orientation n with a pixel index k. Using this more compact notation, we conflate the scaling coefficient map with the wavelet coefficient maps as the ILC method applies in exactly the same fashion. The ILC estimate of the CMB signal at each wavelet scale and orientation is defined as a weighted sum of the wavelet coefficient maps at that scale and orientation

\begin{equation} W_{jnk}^\mathrm{ILC} \equiv \sum _{c=1}^{N_{\rm c}} \omega _{jnk}^c W_{jnk}^c \,, \end{equation}

(11)

where |$\omega _{jnk}^c$| are the weights (which are allowed to vary across the scale and orientation of the signal as well as pixel space) and N_c is the number of input channels.

We impose a constraint on the weights (to ensure that the CMB signal is preserved) such that

\begin{equation} \sum _{c=1}^{N_{\rm c}} a^c \omega _{jnk}^c = 1. \end{equation}

(12)

Assuming that the CMB and foregrounds and the CMB and noise are, respectively, uncorrelated, the variance of the error in the result is minimized when the variance of the ILC map itself is minimized. The resulting weights are given by

\begin{equation} \omega _{jnk}^c = \frac{\sum _{c^{\prime }=1}^{N_{\rm c}} (R_{jnk}^{-1})^{cc^{\prime }} a^{c^{\prime }}}{\sum _{c=1}^{N_{\rm c}} \sum _{c^{\prime }=1}^{N_{\rm c}} a^c (R_{jnk}^{-1})^{cc^{\prime }} a^{c^{\prime }}}, \end{equation}

(13)

where the true covariance matrices at scale j, orientation n and pixel k, |$(R_{jnk})^{cc^{\prime }} = \langle W_{jnk}^c W_{jnk}^{c^{\prime }} \rangle$| (where the angled brackets indicate an ensemble average). For a derivation of equation (13), see Tegmark et al. (2003) and Eriksen et al. (2004).

In this work, we estimate covariance matrices empirically by the following procedure (as used in Basak & Delabrouille 2012; Planck Collaboration IX 2015). We start by calculating at each pixel k:

\begin{equation} (R_{jnk}^\mathrm{approx})^{cc^{\prime }} = W_{jnk}^c W_{jnk}^{c^{\prime }}. \end{equation}

(14)

We then smooth each element of the above matrix by a Gaussian beam w_j(k, k′) in pixel space to form the empirical estimates of covariance matrices

\begin{equation} (\hat{R}_{jnk})^{cc^{\prime }} = \sum _{k^{\prime }=1}^{N_\mathrm{samp}^j} w_j(k,k^{\prime }) (R_{jnk^{\prime }}^\mathrm{approx})^{cc^{\prime }}, \end{equation}

(15)

where |$N_\mathrm{samp}^j$| is the total number of pixels in a given map at scale j. For computational efficiency, we perform this smoothing in harmonic space:

\begin{equation} (\hat{R}_{jnk})^{cc^{\prime }} = \sum _{\ell =0}^{2\ell _\mathrm{max}^j} \sum _{m=-\ell }^{\ell } w_j^\ell (r_{jn}^{\ell m})^{cc^{\prime }} Y_k^{\ell m}, \end{equation}

(16)

where |$(r_{jn}^{\ell m})^{cc^{\prime }}$| are the harmonic coefficients of the maps formed by the elements of matrices |$(R_{jnk}^\mathrm{approx})^{cc^{\prime }}$|⁠, |$w_j^\ell$| is a Gaussian beam transfer function and |$Y_k^{\ell m}$| are the spherical harmonics evaluated at pixel k.

The size of the Gaussian kernel used to smooth the covariance matrices is chosen to be proportional to the size of the wavelet used to form a particular set of wavelet coefficient maps.⁹ In general, the estimation of covariance matrices in ILC methods could be further optimized. It may be preferable to dynamically adapt the smoothing kernel used based on local data. Delabrouille et al. (2009) suggested using a larger kernel at high Galactic latitudes, where Galactic emission does not vary so much and a smaller kernel towards the Galactic equator, where emission is more complex. It could involve masking equatorial regions when estimating the covariance at higher latitudes [somewhat akin to Planck Collaboration XII (2014)]. It could involve convolving the maps of elements of covariance matrices with the same directional wavelet in order to pick out how the local covariance follows the directional structure of the underlying signal. As mentioned above, in this work, we use a similar method as in previous work for ease of comparison.

It is also worth discussing the upper limit on the summation over ℓ in equation (16). We first note the general rule that for the product of two spherical harmonics (Driscoll & Healy 1994)

\begin{equation} Y_{\ell _1,m_1}(\hat{\boldsymbol {n}}) Y_{\ell _2,m_2}(\hat{\boldsymbol {n}}) = \sum _{L = |\ell _1 - \ell _2|}^{\ell _1 + \ell _2} a_{L,m_1+m_2} Y_{L,m_1+m_2}(\hat{\boldsymbol {n}}), \end{equation}

(17)

where |$Y_{L,m_1+m_2}(\hat{\boldsymbol {n}})$| is defined to be zero, if |m₁ + m₂| > L. It follows that for the product of two band-limited maps, |$M(\hat{\boldsymbol {n}}) = \sum _{\ell =\ell _1}^{\ell _2} \sum _{m=-\ell }^\ell m_{\ell m} Y_{\ell m}(\hat{\boldsymbol {n}})$| and |$N(\hat{\boldsymbol {n}}) = \sum _{\ell =\ell _3}^{\ell _4} \sum _{m=-\ell }^\ell n_{\ell m} Y_{\ell m}(\hat{\boldsymbol {n}})$| (where, without loss of generality, ℓ₁ ≤ ℓ₄):

\begin{equation} M(\hat{\boldsymbol {n}}) N(\hat{\boldsymbol {n}}) = \sum _{L=0}^{\ell _2 + \ell _4} \sum _{M=-L}^L p_{LM} Y_{LM}(\hat{\boldsymbol {n}}) \end{equation}

(18)

for ℓ₃ < ℓ₂, i.e. the limits on ℓ in the two maps overlap (the p_LM are the new harmonic coefficients). (If the limits do not overlap, the lower limit on L in equation (18) becomes ℓ₃ − ℓ₂.) The limits on ℓ in the wavelet coefficient maps |$W_{jnk}^c$| are |$(\ell _\mathrm{min}^j,\ell _\mathrm{max}^j)$| [see Section 3.4]. Therefore, as per equations (18) and (15), the limits on ℓ in the covariance matrix element maps |$(\hat{R}_{jnk})^{cc^{\prime }}$| are |$(0,2\ell _\mathrm{max}^j)$|⁠; hence the limits on ℓ in equation (16).

Having established the main equations governing the ILC method, we now present the main steps in the ILC algorithm that we use.

• Form the |$(R_{jnk}^\mathrm{approx})^{cc^{\prime }}$| by equation (14).

• Smooth the |$(R_{jnk}^\mathrm{approx})^{cc^{\prime }}$| in harmonic space by equation (16).

• Take the inverse of each covariance matrix at each pixel to form |$(\hat{R}_{jnk}^{-1})^{cc^{\prime }}$|⁠.

• Calculate the ILC weights |$\omega _{jnk}^c$| by equation (13), where we assume that a^c = 1 for all c and we substitute the empirical estimates for the inverse covariance matrices.

• Finally, calculate the ILC estimate wavelet coefficient maps |$W_{jnk}^\mathrm{ILC}$| by applying equation (11).

3.6 Point source masking

The input frequency maps are diffusively inpainted in a small point source mask following the method employed by Planck Collaboration XVII (2015). This recognizes that the ILC fails when the CMB is obscured by bright extragalactic sources or complex emission near the Galactic equator. The inpainting removes these sources and attempts to replace them with an extrapolation of the surrounding signal. The mask supplied is taken from the NILC section of Planck Collaboration I (2014) and is constructed from the Planck Catalogue of Compact Sources (Planck Collaboration XXVIII 2014; Planck Collaboration XXVI 2015).¹⁰ It masks about 2.2 per cent of the whole sky, predominantly along the Galactic equator towards the Galactic Centre.

Because of this inpainting, the final ILC map is inpainted within the point source mask. For the purposes of this inpainting, we have split the mask into two, based on the size of its constituent individual contiguous holes.¹¹ For holes consisting of less than or equal to 800 pixels, we inpaint with a constrained Gaussian realization following the method of Benoit-Lévy et al. (2013), itself an approximate implementation of the Hoffman–Ribak algorithm (Hoffman & Ribak 1991). For holes consisting of more than 800 pixels (the largest 131 out of 10 031), we inpaint with a standard diffusive algorithm (in particular, following the method employed by Planck Collaboration XVII 2015). The result is that the ILC map is 1.3 per cent Gaussian inpainted and 0.9 per cent diffusively inpainted. This follows Benoit-Lévy et al. (2013), who do not recommend using their Gaussian inpainting for large holes near the Galactic equator.

3.7 Numerical implementation

SILC is implemented in python and is parallelized. At full Planck resolution (N_side = 2048, ℓ_max = 3600), when run on a 60-core symmetric multiprocessor (SMP) with 1.5 TB RAM and a 24-core cluster node with 256 GB RAM,¹² the pipeline takes approximately 12 h per direction. As shown by equation (16), we perform spherical harmonic transforms to 2ℓ_max = 7200. For a given number of directions N, the full pipeline takes approximately N times as long as the axisymmetric limit of our method (when N = 1). In our infrastructure, the code was usually memory-limited (due to the very high resolution of the covariance matrix maps |$(R_{jnk})^{cc^{\prime }}$| at harmonic band-limit |$2 \ell _\mathrm{max}^j$|⁠); the amount of parallelization sometimes had to be reduced to prevent memory overloads on a single node. As mentioned in Section 3.4, the wavelet transforms employ the latest version of s2let (Leistedt et al. 2013; McEwen et al. 2015c), written in c with python wrappers, itself employing ssht (McEwen & Wiaux 2011) and so3 (McEwen et al. 2015c). Despite the use of MW sampling and FFTs, spherical harmonic transforms are the most time-consuming part of the pipeline, again due to the very high resolution of the |$(R_{jnk})^{cc^{\prime }}$| [for the smallest wavelets, these covariance maps are band-limited at ℓ = 7200]. There is scope to further optimize the implementation. Our wavelet analysis and synthesis functions do not, respectively, output and take as input wavelet coefficient maps at double-resolution (i.e. a map band-limited at |$\ell _j^\mathrm{max}$| sampled at |$2 \ell _j^\mathrm{max}$|⁠), requiring additional spherical harmonic transforms to double the resolution. Also, our spherical harmonic transform function does not calculate harmonic coefficients to a multipole less than the band-limit of the input map (i.e. to only calculate a_ℓm for ℓ < L where |$L < \ell _\mathrm{max}^j$|⁠), resulting in excess computation at certain steps in the algorithm. These optimizations are left as further work.

4 SOURCES OF ERROR IN THE ILC

By the linearity of the wavelet transform in equation (7), the data model in equations (1) and (2) can be recast in wavelet space as

\begin{equation} W_{jnk}^c = a^c W_{jnk}^{\mathrm{CMB}} + W_{jnk}^{\mathrm{FG},c} + W_{jnk}^{\mathrm{N},c} \,, \end{equation}

(19)

where |$W_{jnk}^{\mathrm{CMB}}$|⁠, |$W_{jnk}^{\mathrm{FG},c}$| and |$W_{jnk}^{\mathrm{N},c}$| are, respectively, the CMB, foreground and instrumental noise contributions to each wavelet coefficient map. The beams within each component have been absorbed into the component wavelet coefficient maps. Substituting equation (19) into equation (11) gives

\begin{eqnarray} W_{jnk}^{\mathrm{ILC}} &=& \sum _{c=1}^{N_{\rm c}} a^c \omega _{jnk}^c W_{jnk}^{\mathrm{CMB}} + \sum _{c=1}^{N_{\rm c}} \omega _{jnk}^c (W_{jnk}^{\mathrm{FG},c} + W_{jnk}^{\mathrm{N},c}) \nonumber\\ &=& W_{jnk}^{\mathrm{CMB}} + \frac{\sum _{c,c^{\prime }=1}^{N_{\rm c}} (W_{jnk}^{\mathrm{FG},c} + W_{jnk}^{\mathrm{N},c}) (R_{jnk}^{-1})^{cc^{\prime }} a^{c^{\prime }}}{\sum _{c,c^{\prime }=1}^{N_{\rm c}} a^c (R_{jnk}^{-1})^{cc^{\prime }} a^{c^{\prime }}} \,, \end{eqnarray}

(20)

where the second equality follows by applying the constraint given in equation (12) and expanding the weights as given in equation (13). Even when the calibration a^c and the covariance matrices |$(R_{jnk})^{cc^{\prime }}$| are correct, there is always residual signal in the final ILC wavelet coefficient maps, given by the second term on the right-hand side of equation (20). Due to the linearity of the inverse wavelet transforms, this residual signal will propagate linearly into the final ILC temperature map as calculated by equation (10). As explained in Section 3.5, this error term is reduced by minimizing the empirical variance of the ILC map assuming that the CMB and foregrounds and the CMB and noise are, respectively, uncorrelated.

There are additional sources of error in the ILC method. The first is due to inaccuracy in the calculation of covariance matrices |$(R_{jnk})^{cc^{\prime }}$|⁠, i.e. deviations in the empirical estimate |$(\hat{R}_{jnk})^{cc^{\prime }}$| from the true covariance |$(R_{jnk})^{cc^{\prime }}$|⁠. Delabrouille et al. (2009) estimated the first order expansion of the reconstruction error in the ILC map estimate due to this covariance error. They showed that the covariance of the ILC error with the CMB is inversely proportional to the number of ‘effective modes’ used in the ILC calculation. This covariance bias is negative. In our directional wavelet decomposition, our ‘effective modes’ are spherical harmonic coefficients weighted to take account of the fact that the harmonic responses of wavelets overlap in both scale and direction. As N, the number of orientations probed, increases and so does the number of wavelets, each wavelet coefficient map contains fewer ‘effective modes’. We therefore conclude that our directional wavelet ILC may be susceptible to this negative ILC bias by increasing N. Delabrouille et al. (2009) also showed that due to chance correlations between the CMB and foregrounds, the variance minimization leads to the unintentional cancellation of N^c − 1 CMB modes. For Planck, N^c = 9, whereas for WMAP, N^c = 5. We therefore expect the magnitude of this negative bias to double simply by using more input frequency channels. Also, since this covariance bias is due to the cancellation of CMB modes, Delabrouille et al. (2009) showed that the absolute value is proportional to the CMB power. Therefore, the absolute value of the bias is greatest on large scales where CMB power is concentrated. In general, these biases are best estimated through suites of Monte Carlo simulations.

Another source of error is due to inaccuracy in the calibration a^c of the CMB. Dick, Remazeilles & Delabrouille (2010) calculated the consequence of a first-order error in a^c on a multiplicative correction to the CMB term in equation (20). They showed that even a small error in calibration can lead to a significant negative multiplicative bias in the CMB term, when the signal-to-noise ratio is large. (Here, the noise in this ratio also includes foreground signal.) They consider the implications for using an ILC on Planck data, where the signal-to-noise ratio is larger than for WMAP data. They estimate that a 1 per cent error in a^c can cancel about a third of the CMB signal, while even a 0.1 per cent error in a^c can remove about 1 per cent of the CMB. Since our main map products use Planck data as input, they will be susceptible to this additional negative calibration bias. In this work, we assume that the CMB is calibrated to have unit response for all frequency channels, i.e. a^c = 1 for all c.

As mentioned in Sections 3.2 and 3.3, we assume all beams to be circularly symmetric. Therefore, non-axisymmetric beam components will propagate into the ILC calculation but are assumed sufficiently small to be ignored.

5 COMPARISON TO PREVIOUS WORK

We now consider how SILC compares with existing component separation methods, particularly those adopted for the Planck analysis. We applied the axisymmetric limit (when N = 1) of SILC to full-mission Planck data and compared the results to existing Planck analyses using the NILC and SMICA methods: the former because it is the closest in spirit to SILC, and the latter because it is the baseline method adopted by the Planck Collaboration for high-resolution analyses. Fig. 4 shows the CMB reconstructed by the axisymmetric limit of SILC, while Fig. 5 shows the differences between this map and the NILC and SMICA (full-mission 2015 release) CMB maps and the difference between NILC and SMICA. The differences between the three maps are small in magnitude and mostly concentrated at the edges of the Galactic mask towards the Galactic Centre, where foreground emission is most intense and complex. Quantitatively, we can compare the mean values and standard deviations of the masked difference maps. The mean values of Figs 5(a), (b) and (c) are, respectively, 0.44, −0.63 and −1.07 μK, while the standard deviations are, respectively, 4.24, 3.38 and 3.43 μK². These values are small and similar, suggesting a strong consistency between the three methods. These difference maps have been formed from maps which have been smoothed and downgraded in resolution and so visually highlight differences at the lowest multipoles.

Figure 4.

Planck data. The CMB temperature anisotropies reconstructed using SILC in the axisymmetric limit (N = 1, FWHM = 5 arcmin, N_side = 2048). The grey pixels are the point source mask.

Figure 5.

Planck data. Differences between the axisymmetric limit (N = 1) of SILC, NILC and SMICA. The maps have been smoothed to FWHM = 80 arcmin and downgraded to N_side = 128. The grey pixels are the UT78 confidence mask from Planck Collaboration IX (2015), which masks the regions of the NILC and SMICA maps not recommended for cosmological analysis. The differences are [from top to bottom] (a) SILC (N = 1) − NILC, (b) SILC (N = 1) − SMICA and (c) NILC − SMICA.

Fig. 6 compares point source masked TT angular power spectra (D_ℓ = ℓ(ℓ + 1)C_ℓ/2π) at the full multipole range of the three maps (up to ℓ = 3400)¹³ with a CMB spectrum derived from the Planck 2015 TT and low TEB likelihood.¹⁴ The SILC spectrum is remarkably similar to that of NILC. This is unsurprising since the axisymmetric limit of SILC (when N = 1) is very similar to the NILC method. None the less, there are a number of pipeline differences. In particular, we use a different set of wavelets than the needlets employed in NILC (as discussed in Section 3.4), even in the axisymmetric limit, with different harmonic responses. Fig. 2 shows the harmonic response of the wavelets used in this work and Table 1 lists their harmonic band-limits |$\ell ^j_{\mathrm{min}}$| and |$\ell ^j_{\mathrm{max}}$|⁠. The SMICA spectrum has lower residuals at higher multipoles than both the axisymmetric limit of SILC and NILC.

Figure 6.

Planck data. TT angular power spectra comparing the axisymmetric limit (N = 1) of SILC to NILC and SMICA. The top panel (a) shows point source masked spectra. The middle panel (b) shows residuals after subtracting the best-fitting Λ cold dark matter (ΛCDM) model from the Planck 2015 likelihood. The bottom panel (c) shows the same residuals at low multipoles only (ℓ < 1500).

Fig. 7 compares full-sky angular power spectra of the three maps, including the inpainted point source regions. The spectra are similar to those in Fig. 6. The main difference is the lower noise tail in the SILC map at high multipoles above ℓ = 1500 (where all component separation CMB maps are dominated by residual instrumental noise). This is because, unlike NILC and SMICA, we do not Gaussian inpaint the very largest point source holes, but rather use diffusive inpainting (as discussed in Section 3.6). The Gaussian inpainting of large irregular holes is poorly constrained and adds residual noise relative to diffusive inpainting.

Figure 7.

Planck data. TT angular power spectra comparing the axisymmetric limit (N = 1) of SILC to NILC and SMICA. The top panel (a) shows full-sky spectra of inpainted maps. The middle panel (b) shows residuals after subtracting the best-fitting ΛCDM model from the Planck 2015 likelihood. The bottom panel (c) shows the same residuals at low multipoles only (ℓ < 1500).

We have shown that the axisymmetric limit of SILC gives comparable performance to NILC and SMICA. In Sections 6 and 7, we ‘turn on’ the directionality of the wavelets and consider the impact on CMB reconstruction from simulated and real data, respectively.

6 Application to PLANCK simulations

We now apply SILC to the fiducial full-mission Planck FFP8 simulated sky maps, focusing on the impact on CMB reconstruction by increasing directionality as a function of scale. Fig. 8 shows the difference between our reconstructed CMB (using N = 1) and the input simulated CMB. There are small-magnitude differences particularly at the edge of the Galactic mask where the strength and complexity of foreground emission is greatest. As in Fig. 5, this difference map is at low resolution and so highlights residuals at the lowest multipoles. Fig. 9 compares point source masked TT angular power spectra (up to ℓ = 3400) of CMB maps reconstructed using values of N from 1 to 5. It can be seen that the introduction of directionality has the greatest effect at multipoles around ℓ = 800; the residuals are beginning to converge for |$\ell \gtrapprox 2000$|⁠. Fig. 10 shows the differences between simulated CMB maps reconstructed using N = [2, 3, 4, 5] minus the input CMB. The four maps and the axisymmetric difference map in Fig. 8 are almost identical with small magnitude residuals. This is because these low-resolution difference maps again highlight residuals on the very largest scales. However, as discussed in Section 3.4, the wavelets we use are constructed to have an axisymmetric scaling function at the very lowest multipoles. The scaling function we use (as detailed in Table 1) means that no directionality is applied for ℓ < 32.

Figure 8.

Planck simulations. Difference between output ILC and input CMB temperature maps from FFP8 simulations. The maps have been smoothed to FWHM = 80 arcmin and downgraded to N_side = 128. The grey pixels are the UTA76 confidence mask from Planck Collaboration IX (2015), which masks the Galactic region in FFP8 simulations where foreground emission is strongest.

Figure 9.

Planck simulations. TT angular power spectra comparing output ILC using different values of N and input CMB from FFP8 simulations. The top panel (a) shows point source masked spectra. The middle panel (b) shows residuals after subtracting the input CMB spectrum. The bottom panel (c) shows the same residuals at low multipoles only (ℓ < 1500).

Figure 10.

Planck simulations. Differences between output ILC reconstructed using different values of N and input CMB temperature maps from FFP8 simulations. The maps have been smoothed to FWHM = 80 arcmin and downgraded to N_side = 128. The grey pixels are the UTA76 confidence mask. The differences are [ from left to right, top to bottom] (a) N = 2, (b) N = 3, (c) N = 4, (d) N = 5 minus the input CMB.

Fig. 11 shows equivalent difference maps as in Fig. 10 but for the simulated CMB reconstructed using directional wavelets at all scales,¹⁵ including for ℓ < 32. It can be seen that the reconstruction errors are significantly larger in magnitude and cover almost the entire sky. The errors are also dominated by the largest scales, in particular a large error in the quadrupole increasing with magnitude as the amount of directionality N increases. We attribute this effect most probably to the ILC ‘biases’ discussed in Section 4, in particular the cancellation of CMB modes due to chance correlations with foregrounds in the ILC variance minimization. Delabrouille et al. (2009) showed that the absolute value of this effect is largest on large scales where CMB power is concentrated, since the cancelled CMB modes on large scales have the greatest magnitude. Further, as discussed in Section 4 and shown in Fig. 11, these errors are expected to increase in magnitude as a function of N. This is because as N increases, each directional wavelet coefficient map (the space in which our ILC operates) contains fewer ‘effective modes’ of the input data and so the error in our empirical covariance estimation is expected to increase. This error propagates to the final maps.

Figure 11.

Planck simulations. Same as Fig. 10 (which uses the recommended wavelets) but here using directional wavelets on large scales (ℓ < 32), which is not recommended as it leads to increased CMB reconstruction errors as seen above.

These map reconstruction errors due to the implementation of directionality on the very largest scales are accompanied by increasingly negative power spectrum residuals as N increases, in particular in the first multipole bin from ℓ = 2 to ℓ = 11. This is also indicative of the negative ILC bias due to empirical CMB cancellation, as discussed in Section 4 and Delabrouille et al. (2009). The results in Fig. 11 thus motivate the use of an axisymmetric scaling function, which ensures that no directionality is used for ℓ < 32 and so reduces the errors in CMB reconstruction. In principle, these biases can be estimated and corrected through large suites of simulations, which is beyond the scope of this work.

7 Application to PLANCK data

We now study the application of SILC with increasing directionality to the full-mission Planck sky maps. The left column of Fig. 12 shows the full-resolution reconstructed CMB maps as calculated with different values of N from 2 to 5, which visually appear very similar. The right column of Fig. 12 shows the differences between the CMB reconstructed using N = [2, 3, 4, 5] minus the axisymmetric limit (when N = 1), highlighting the differences at the lowest multipoles. The differences are of largest magnitude towards the Galactic plane where foreground emission is concentrated. This shows how the different wavelet kernels are localizing the ILC weights differently in response to the directional structure of the foregrounds and CMB. The differences are small, reflecting the implementation of an axisymmetric scaling function, meaning that no directionality is applied at ℓ < 32. Fig. 13 compares point source masked TT angular power spectra of the CMB reconstructed using values of N from 1 to 5. The power spectrum residuals from a Planck best-fitting ΛCDM model remain small for most scales until the reconstructed spectra reach a characteristic noise spectrum for |$\ell \gtrapprox 1500$| where the different values of N converge. At these high multipoles, the ILC solution is dominated by residual instrumental noise. We see the biggest impact from directionality at intermediate multipoles (from ℓ = 400 to ℓ = 1500). For comparison, we plot the SMICA power spectrum. In further support to the discussion in Section 5, SILC matches the performance of SMICA. We note that, as with the simulations in Section 6, directionality changes the reconstructed CMB power spectrum most significantly at intermediate multipoles around ℓ = 800.

Figure 12.

Planck data. Left: CMB temperature anisotropies reconstructed using SILC with different values of N (FWHM = 5 arcmin, N_side = 2048). Right: differences between CMB temperature maps reconstructed using different values of N minus the axisymmetric limit N = 1. The maps have been smoothed to FWHM = 80 arcmin and downgraded to N_side = 128. In both columns: the grey pixels are the point source mask (downgraded in resolution as appropriate). From top to bottom: (a) N = 2, (b) N = 3, (c) N = 4, (d) N = 5.

Figure 13.

Planck data. TT angular power spectra comparing different values of N from 1 to 5 and SMICA. The top panel (a) shows point source masked spectra. The middle panel (b) shows residuals after subtracting the best-fitting ΛCDM model from the Planck 2015 likelihood. The bottom panel (c) shows the same residuals at low multipoles only (ℓ < 1500).

8 DISCUSSION

The comparisons in Section 5 demonstrate that SILC matches the performance of two previous methods, NILC and SMICA, in both maps and power spectra, with particular similarity between the axisymmetric limit of SILC and NILC, as expected. Both map residuals and power spectra in Sections 6 and 7 show that switching on directionality changes CMB reconstruction most significantly at intermediate multipoles ℓ = 400–1500. There appears to be little benefit in localizing the ILC with directional wavelets at the very smallest scales, where the ILC result is noise-limited. We also adopt an axisymmetric scaling function on the very largest scales, meaning that there is no directionality at ℓ < 32. In Fig. 11, we show the large CMB reconstruction errors arising from using directionality on the largest scales. This motivates the use of an axisymmetric scaling function, which significantly reduces the errors as seen in Fig. 10. In Section 6, we sketched out an argument that attributes these errors to empirical CMB cancellation (Section 4). However, the precise source and exact magnitude of any ILC errors are best estimated through suites of simulations.

We have presented this analysis by producing CMB maps (in Sections 6 and 7) each with a different single value of N at all wavelet scales. Our method can be simply extended to allow different values of N at each wavelet scale. In the same way that each wavelet scale has different harmonic band-limits, they can also have different azimuthal band-limits, optimized as identified above to reduce foreground and noise residuals.

The negative ILC power spectrum biases discussed in Section 4 must be quantified in parallel to this directionality optimization if using the resulting map for power spectrum analyses. It is possible to estimate variance biases in the data through suites of realistic simulations. However, we can also calculate this using the data itself and a fiducial CMB spectrum. In wavelet space, the variance estimator at each wavelet coefficient is |$\langle W^\mathrm{ILC} W^{\mathrm{ILC}\dagger }\rangle = \boldsymbol {\omega }^\dagger \langle \boldsymbol {W} \boldsymbol {W}^\dagger \rangle \boldsymbol {\omega } = \boldsymbol {\omega }^\dagger \bf{R} \boldsymbol {\omega } = (\boldsymbol {a}^\dagger \bf{R}^{-1} \boldsymbol {a})^{-1} = (\sum _{c,c^{\prime }=1}^{N_{\rm c}} (R^{-1})^{cc^{\prime }})^{-1}$|⁠, where each equality follows by applying in turn equations (11), (14) and (13) [from Section 3.5] and then expanding the vector notation [the vectors span the space defined by the number of input channels; explicitly, we assume unit CMB calibration |$\boldsymbol {a} = (1,1,1,1,1,1,1,1,1)$|]. In order to calculate the variance bias in the ILC, we can subtract the expected CMB variance |$\sum _{\ell m} C_\ell {| \Psi _{\ell m}^j |}^2$|⁠. If analysis was done in wavelet space, the above would define the variance bias. If only considering the diagonal terms in wavelet space, it is possible to straightforwardly transform this estimate to real space through an inverse wavelet transform as in equation (10), substituting |$\Phi (\hat{\boldsymbol {n}})$| and |$\Psi ^j (\hat{\boldsymbol {\rho }}),$| respectively, for |${| \Phi (\hat{\boldsymbol {n}}) |}^2$| and |${| \Psi ^j (\hat{\boldsymbol {\rho }}) |}^2$|⁠. However, for a full treatment of the variance bias including off-diagonal terms, full wavelet space covariances need to be calculated. Although many off-diagonal terms would decay, this would still be computationally demanding and will be the focus of future research. However, we reiterate that if analysis is carried out in wavelet space, then variance biases can be straightforwardly calculated from the information already contained in the results.

9 CONCLUSIONS

We have presented SILC, a new form of ILC that uses directional, scale-discretized wavelets to localize the ILC weighting according to the frequency, spatial, harmonic and, for the first time, morphological information in the CMB and its foregrounds. This is motivated by the anisotropic or filamentary morphology of both the CMB and astrophysical foregrounds in the microwave sky. We have tested SILC on 2015-release Planck data and simulations, demonstrating comparable performance to two existing component separation algorithms, NILC and SMICA, and investigated how to optimize the use of morphological information through directionality. We have explored increasing the amount of directionality in the algorithm, showing that on the largest and the smallest scales, the axisymmetric limit of the ILC works well, while at intermediate multipoles (from ℓ = 400–1500), increasing N (the number of directions per scale) leads to lower residuals. At high multipoles (⁠|$\ell \gtrapprox 1500$|⁠), the input data are already noise-limited, as is the ILC reconstruction, and directionality does not reduce the reconstruction error, as instrumental noise has no directional structure. We adopt an axisymmetric scaling function to analyse the largest scales, meaning that we use no directionality for ℓ < 32. This is motivated by the observation that increasing directionality on large scales gives increased reconstruction errors over the axisymmetric limit. We argue that these errors are due to empirical CMB cancellation in the ILC calculation, though the exact source must be estimated through large suites of realistic simulations. Allowing N to vary with wavelet scale is analogous to the choice of different harmonic band-limits at different scales.

We conclude that the introduction of directional wavelets allows greater flexibility in the ILC to make use of morphological information at targeted scales. Our multiprocessing implementation takes advantage of the wavelet scales to allow large-scale results to be analysed while small scales are still being processed. Moreover, our wavelet transforms are quick and exact, using MW sampling and FFTs (see Section 2). As discussed in Section 4, the ILC is prone to several sources of error and variance bias, including empirical CMB cancellation. This bias can be estimated through suites of Monte Carlo simulations, but we have also outlined (in Section 8) the ability to estimate biases directly from the data, most straightforwardly in wavelet space. We make our map products available at http://www.silc-cmb.org.¹⁶

This work on scalar signals (i.e. the temperature I component of the CMB) can be extended to spin signals (i.e. the polarization Q and U components of the CMB, or, equivalently, the E and B modes), by using spin wavelets (McEwen et al. 2014, 2015b; Leistedt et al. 2015). These are an extension of directional, scale-discretized wavelets to represent spin signals, such as CMB polarization, a spin ±2 signal. We expect that the directionality will be particularly suited to the anisotropic, filamentary nature of polarized foregrounds when observed on the sky, and in future work will present the application of SILC to CMB polarization data.

KKR thanks Franz Elsner for valuable discussions. KKR was supported by the Science and Technology Facilities Council. HVP and BL were partially supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 306478-CosmicDawn. JDM was partially supported by the Engineering and Physical Sciences Research Council (grant number EP/M011852/1). AP was supported by the Royal Society. 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.

1

Planck 2015 Release Explanatory Supplement: UC CC tables (http://wiki.cosmos.esa.int/planckpla2015/index.php/UC_CC_Tables). For the 545 GHz map, the unit conversion is (58.0356 ± 0.0278) MJy Sr⁻¹ K_CMB^{− 1} and for the 857 GHz map, the unit conversion is (2.2681 ± 0.0270) MJy Sr⁻¹ K_CMB^{− 1}.

2

http://pla.esac.esa.int/pla

3

FFP8 simulations are also available with bandpass mismatch, accounting for differences in the bandpasses of detectors nominally at the same frequency, leading to spurious signals in the frequency maps.

4

Planck 2015 Release Explanatory Supplement: the 2015 instrument model (http://wiki.cosmos.esa.int/planckpla2015/index.php/The_RIMO).

5

We adopt the zyz Euler convention corresponding to the rotation of a physical body in a fixed coordinate system about the z, y and z axes by χ, θ and ϕ, respectively.

6

http://www.s2let.org

7

http://www.spinsht.org

8

http://www.sothree.org

9

|$\mathrm{FWHM}^j = 50 \sqrt{\frac{1200}{N_\mathrm{samp}^j}}$|⁠. This value is the same as used in the NILC implementation on Planck data.

10

The details of its construction are given in Planck Collaboration I (2014). It can be downloaded from http://pla.esac.esa.int/pla and is labelled I_MASK in the NILC data products.

11

Note that many holes can be large and irregularly shaped due to the overlapping of smaller circular holes.

12

The exact specification for our infrastructure is an Intel Xeon E7-4890 2.8 GHz SMP with 4 × 15-core CPUs with 25.6 GB RAM per core, and an Intel Xeon E5-2697 2.7 GHz node with 2 × 12-core CPUs with 10.7 GB RAM per core.

13

In order to estimate full-sky spectra from a masked map, we correct the C_ℓ by dividing by f_sky = 0.978, a good approximation for a small mask. We elect to use point source masked spectra in order to concentrate our analysis on foreground and noise removal, rather than how maps are inpainted; all three maps are inpainted (at least) within the mask used.

14

The parameters come from the base_plikHM_TT_lowTEB likelihood. The values are available in the Planck 2015 Release Explanatory Supplement: 2015 Cosmological parameters and MC chains (http://wiki.cosmos.esa.int/planckpla2015/images/f/f7/Baseline_params_table_2015_limit68.pdf).

15

In particular, the scaling function and j = 0 wavelet are replaced by two directional wavelets with harmonic band-limits [1, 60] and [1, 128].

16

The DOI for our data release is 10.5281/zenodo.44373.