SPECTRAL IMAGING OF GALAXY CLUSTERS WITH PLANCK

Across the radio electromagnetic spectrum, the SZ signal is a mixture of CMB radiation, and several extended emissions of Galactic (free–free, synchrotron and thermal dust continua, CO lines) and extragalactic origin (cosmic infrared background). When looking at all-sky gigahertzian frequency maps ${\boldsymbol{f}}$that are contaminated by their instrumental noise ${\boldsymbol{N}}$, separating these components is an inverse problem that is traditionally stated as:

$$\begin{eqnarray}&&{\boldsymbol{f}}={\boldsymbol{As}}+{\boldsymbol{N}},\end{eqnarray} \tag{ 1 }$$

where we look for a mixing matrix ${\boldsymbol{A}}$ between several unknown emitting sources ${\boldsymbol{s}}$. Originally motivated by the denoising of all-sky CMB maps, a number of component separation algorithms have been proposed for solving Equation (1):

• 1.  
  parametric methods where the spectral energy distributions (SEDs) of all components are used as a priori knowledge to maximize their spatial distribution entropy or likelihood (Hobson et al. 1998; Stolyarov et al. 2005; Eriksen et al. 2008; Khatri 2015);
• 2.  
  internal linear combination methods that map a specific component as a linear combination of the observed maps with the constraint of minimizing its variance (Eriksen et al. 2004; Remazeilles et al. 2011; Hurier et al. 2013); and
• 3.  
  blind or semi-blind source separation methods that rely on spatial source independence or sparsity priors (Delabrouille et al. 2003; Cardoso et al. 2008; Bobin et al. 2013).

Unlike the image of the last scattering surface, some of the foreground sources are localized objects that may not be uniformly modeled by Gaussian random fields. Therefore, one of the challenges of the component separation algorithms is to avoid any unreal amplification of the very low or null amplitude elements of the source vector, ${\boldsymbol{s}}$. In this context, parametric methods have shown their robustness for mapping the brightest extended foregrounds (see, e.g., Planck Collaboration 2014c), while sparsity-regularized reconstruction approaches have proven their ability to spatially separate the CMB anisotropies from fainter and more localized objects (Bobin et al. 2014).

Among sparse representations, wavelet transforms have long been successfully used to denoise X-ray images of galaxy clusters (Slezak et al. 1994; Vikhlinin et al. 1997; Starck & Pierre 1998), and thus have naturally been proposed for SZ imaging (Pierpaoli et al. 2005; Pires et al. 2006). The isotropy of 2D wavelet functions, however, is unsuitable for detecting and preserving the filamentary structures that are likely to populate the cluster outskirts as well as the the Galactic dust. A solution can be represented by curvelets that are especially designed to perform a sparse representation of linear and curved edges in the images (Candès & Donoho 2000). To combine the local robustness of parametric methods to the imaging capabilities offered by sparse representations, we propose demixing the SZ signal with respect to the CMB and thermal dust via a set of kernel-weighted likelihood maximizations that are directly related to the wavelet transform of the CMB map and to the curvelet transforms of the foreground maps.

The Planck HFI is composed of bolometers array whose frequency channels cover the most prominent part of the thermal SZ (tSZ) energy spectrum. In the HFI frequency range (100–857 GHz), the tSZ component is locally contaminated by radio and infrared point sources, and is mostly mixed up with the CMB temperature anisotropies and with the Galactic extended emissions from the thermal dust and CO lines (Planck Collaboration 2015a). Moreover, most of the time the latter component is negligible at high Galactic latitude, as it is a tracer of Galactic molecular clouds. For these reasons, we reduce the specific intensity expected for each pixel (k, l) of the HFI frequency maps to the sum of three SEDs, associated with the tSZ, thermal dust, and CMB components respectively,

$$\begin{eqnarray}I(k,l,\nu ) & = & R(\nu )\times \left[{A}_{\mathrm{SZ}}(k,l){f}_{\mathrm{SZ}}(\nu )\right.\\ & & +{A}_{\mathrm{dust}}(k,l){f}_{\mathrm{dust}}(\nu )\\ & & +\left.{A}_{\mathrm{CMB}}(k,l)\right],\end{eqnarray} \tag{ 2 }$$

where:

• 1.  
  the Galactic dust spectrum is idealized as a uniform modified blackbody with shape ${f}_{\mathrm{dust}}(\nu )\propto {\nu }^{\beta }{B}_{\nu }(\nu ,T);$
• 2.  
  the tSZ distortion of the CMB spectrum is modeled following the Kompaneets non-relativistic approximation:

  $$\begin{eqnarray}&&{f}_{\mathrm{SZ}}(\nu )=\left[\left(\displaystyle \frac{h\nu }{{kT}}\right)\displaystyle \frac{\mathrm{exp}(h\nu /{kT}+1)}{\mathrm{exp}(h\nu /{kT}-1)}-4\right];\end{eqnarray} \tag{ 3 }$$
• 3.  
  each SED is corrected for the HFI spectral response R(ν).

Note that the correction of each SED for the HFI spectral response includes a unit conversion factor between the HFI 100–353 GHz channels and the 545 and 857 GHz channels, which are calibrated in units of CMB temperature and intensities of a power-law SED, respectively (see details in Planck Collaboration 2014a, 2014b). A color correction is further applied to adapt the Galactic dust SED to the power law used to calibrate the high-energy channels. These corrections are calculated using the Unit conversion and Color Correction package that is provided with the current Planck data release.

To perform a local regression of each component amplitude, we minimize a kernel-Weighted Least-Square (WLS) distance separating the total SED of all components, $I(k,l,\nu ),$ from the energy distribution registered in the vicinity of each pixel of the HFI frequency maps, y(k, l). Derived from the weighted likelihood formalism (see the appendix), the minimized quantity wχ²(k, l) is weighted by the variance of the frequency maps, σ(k, l, ν), and spatially smoothed by a positive kernel with norm 1, w(i, j):

$$\begin{eqnarray}w{\chi }^{2}(k,l) & = & \displaystyle \sum _{i,j,\nu }w(i,j)\\ & & \times \displaystyle \frac{{\left[I(i-k,j-l,\nu )-y(i-k,j-l,\nu )\right]}^{2}}{\sigma {(i-k,j-l,\nu )}^{2}}.\end{eqnarray} \tag{ 4 }$$

The minimization of wχ²(k, l) provides us with an estimate of the searched parameter parameter vector, as follows:

$$\begin{eqnarray}&&\widehat{{{\boldsymbol{A}}}_{w}}(k,l)=\left[\begin{array}{c}\widehat{{{A}_{\mathrm{CMB}}}_{,w}}(k,l)\\ \widehat{{{A}_{\mathrm{dust}}}_{,w}}(k,l)\\ \widehat{{{A}_{\mathrm{SZ}}}_{,w}}(k,l)\end{array}\right]=\underset{{\boldsymbol{A}}}{\mathrm{arg}\;\mathrm{min}}\left(w{\chi }^{2}(k,l)\right).\end{eqnarray} \tag{ 5 }$$

It is undertaken for each pixel of the component maps, $\widehat{{{\boldsymbol{A}}}_{w}}(k,l),$ by means of a Levenberg–Marquardt algorithm (Markwardt 2009).

2.3.1. Image Denoising

Multi-scale transforms perform a sparse representation of the 2D discontinuities that are present in the images, thus enhancing their S/N. The efficiency of the process, however, depends on the level of correlation of the transform basis functions to the shape of the discontinuities. For instance, wavelet transforms best represent isotropic details, while, compared to wavelets, curvelets are better for revealing linear and curved edges. Because CMB features are believed to be nearly isotropic, we chose to analyze the CMB map via its wavelet transform. To analyze the Galactic dust and tSZ images that plausibly hold filaments and edges, we adopted a curvelet transform.

The WLS minimization provides us with a mechanism for extracting the undecimated B₃-spline  wavelet transform (Starck et al. 2007) of all component maps, $A$(k, l). For each wavelet scale a, we split the B₃-spline  wavelet function ψ(a) as the sum of its negative and positive parts, ψ₊(a) and ψ₋(a) (see also Figure 1) and introduce the appropriate weighting kernels in Equation (4). Minimizing ψ₊χ²(k, l) and ψ₋χ²(k, l) provides us with an estimate of the component maps convolved with ψ₊(a) and ψ₋(a), and allows us to derive their wavelet transform, following:

$$\begin{eqnarray}&&{\mathcal{W}}(k,l,a)=\displaystyle \frac{1}{2}\times \left[\widehat{{{\boldsymbol{A}}}_{{\psi }_{+}}}(k,l)+\widehat{{{\boldsymbol{A}}}_{{\psi }_{-}}}(k,l)\right].\end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

In the same way, the variance of the wavelet coefficients can be derived from:

$$\begin{eqnarray}&&d{\mathcal{W}}(k,l,a)=\displaystyle \frac{1}{2}\times \sqrt{\widehat{{{\boldsymbol{\sigma }}}_{A,{{\boldsymbol{\psi }}}_{+}}^{2}}{(k,l)}^{2}+\widehat{{{\boldsymbol{\sigma }}}_{A,{{\boldsymbol{\psi }}}_{-}}^{2}}{(k,l)}^{2}}.\end{eqnarray} \tag{ 7 }$$

To further extract the curvelet coefficients of the component maps, we normalize wavelet coefficients to their variance level, ${\mathcal{W}}^{\prime} ={\mathcal{W}}/d{\mathcal{W}}.$ Subsequently, we compute a so-called curvelet transform of the first generation (Starck et al. 2003) that combines a 2D B₃-spline wavelet transform with a set of ridgelet transforms that are performed within binning blocks of the wavelet detail images. Finally, the denoising of wavelet and curvelet transforms is performed via a thresholding of the coefficients whose absolute value does not exceed the variance. In the case of the CMB wavelet transform, the coefficient variance is straightforwardly defined in Equation (7). On the other hand for the Galactic dust and the tSZ curvelet transforms, the coefficient variance depends on the orientation and it has been preliminarily estimated and tabulated by Monte-Carlo simulations of a Gaussian white noise. In order to work with positive component maps, lower the dynamics of the reconstructed tSZ images, and reduce the relative amplitude of the thresholding artifacts, the quantities that have been effectively denoised and mapped are

$$\begin{eqnarray}&&{{\boldsymbol{A}}}_{o}(k,l)=\left[\begin{array}{c}{{\mathcal{R}}}_{{\mathcal{W}}}\bar{{\mathcal{W}}}[1+{A}_{\mathrm{CMB}}(k,l)]\\ {{\mathcal{R}}}_{{\mathcal{C}}}\bar{{\mathcal{C}}}[{A}_{\mathrm{dust}}(k,l)]\\ {{\mathcal{R}}}_{{\mathcal{C}}}\bar{{\mathcal{C}}}[\mathrm{ln}(1+{A}_{\mathrm{SZ}})(k,l)]\end{array}\right],\end{eqnarray} \tag{ 8 }$$

where $\bar{{\mathcal{W}}}$ and $\bar{{\mathcal{C}}}$ denote the thresholded wavelet and curvelet transforms, respectively, and ${{\mathcal{R}}}_{{\mathcal{C}}}$ and ${{\mathcal{R}}}_{{\mathcal{C}}}$ denote the relative reconstruction operators.

2.3.2. Image Restoration

Because the denoised component maps ${{\boldsymbol{A}}}_{o}(k,l)$ have been estimated from the mixing of frequency maps that were registered at various angular resolutions, they might be altered by aliasing artifacts. In order to take advantage of the angular resolution available at each frequency, we iteratively deconvolve all component maps from the HFI beams by means of a Van Cittert algorithm (van Cittert 1931). More precisely, at each iteration i we associate a set of denoised frequency maps with the component maps ${{\boldsymbol{A}}}_{i}$ following Equation (2). We convolve these frequency maps with their specific beams and restore a new set of point-spread function (PSF)-convolved component maps $P\ast {{\boldsymbol{A}}}_{n}$ via a new combination of WLS estimate and multi-scale denoising. Starting from a first guess ${{\boldsymbol{A}}}_{o}(k,l),$ the Van Cittert iteration yields:

$$\begin{eqnarray}{{\boldsymbol{R}}}_{n}(k,l) & = & {{\boldsymbol{A}}}_{o}(k,l)-[P\ast {{\boldsymbol{A}}}_{n}](k,l)\\ {{\boldsymbol{A}}}_{n+1} & = & {{\boldsymbol{A}}}_{n}+\alpha {{\boldsymbol{R}}}_{n}(k,l).\end{eqnarray} \tag{ 9 }$$

As proposed in Murtagh et al. (1995), the locus of the significant (i.e., nonzero) wavelet and curvelet coefficients of ${{\boldsymbol{A}}}_{o}$ defines a multi-resolution support' that regularizes the algorithm. By reconstructing the residual ${{\boldsymbol{R}}}_{n}(k,l)$ from the wavelet and curvelet coefficients located within this multiresolution support (see also Starck et al. 2002), we prevent any noise amplification during the Van Cittert iteration.

In order to test the restoration algorithm detailed in the previous section, we extracted thermal SZ maps of simulated galaxy clusters. The hydro-simulations employed are carried out by the Smoothed-Particle-Hydrodynamic code GADGET (Springel 2005). The details of the simulations and the physics adopted are presented in Planelles et al. (2014) and Rasia et al. (2014). Here, we briefly list the essential points and address the mentioned papers for a more exhaustive description. The clusters used in the present work come from two Lagrangian regions that were selected from a parent dark matter cosmological box. The regions are re-simulated at a higher-mass resolution, adding baryons either in gaseous or stellar form. The re-simulations treat processes such as radiative cooling; star formation and evolution; kinetic feedback by type Ia and type II supernovae and from asymptotic Giant Branch stars; thermal feedback from active galactic nuclei resulting from gas accretion onto supermassive black holes (see also Ragone-Figueroa et al. 2013). Masses of dark matter and gas particles are equal to m_dm = 8.47 × 10⁸h⁻¹ M_⊙ and m_gas = 1.53 × 10⁸h⁻¹M_⊙, respectively. The adopted Plummer-equivalent softening length for computing the gravitational force is set to  = 5h⁻¹ kpc in comoving units up to z = 2 then it is switched to the same value but in physical units. The minimum SPH smoothing length is 0.5 × .

To prepare a few test cases that are representative of nearby clusters observed by Planck, we select three situations occurring in a rich and dense environment and displaced them at various distances from the observer. The first and second situations are two different time snapshots of the same Lagrangian region, referring to z = 0.25 and z = 0, respectively. The third case corresponds to a different region.

• 1.  
  A massive accreting cluster. With a binding mass of M₂₀₀ = 7 × 10¹⁴h⁻¹M_⊙, this cluster exhibits a rather regular morphology in the innermost regions. However, the cluster periphery is more disturbed. In particular, the object has recently accreted a satellite that is currently lying at a projected distance of 2 Mpc. The smaller object has a temperature of about 2–3 keV, corresponding to a mass of a few 10¹⁴h⁻¹M_⊙. Following its cosmic evolution, we see that it will orbit on the south to the main halo with a projected impact parameter of about 500 kpc, and will lose most of its gas content. In the cosmic time selected, the satellite and main cluster are connected by a tenuous and irregular filament with Compton parameter of y ∼ 10⁻⁵.
• 2.  
  A connected cluster pair. After the merger described above and after having stripped most of the gas of the subclump, the main halo is approached by another massive object of comparable mass. At redshift zero the two systems lie on the same plane of the sky and are about 1.5 Mpc apart. By studying their evolutions, we observed that they are moving toward each other and in the future they will merge. At the current epoch, they are connected by a bright filament with a Compton parameter of y ∼ 3 × 10⁻⁵.
• 3.  
  A colliding cluster system. This system is composed of multiple objects that are all located in the plane of the sky along a large-scale filamentary structure. The main cluster is at the center of the image and has a mass of M₅₀₀ ∼ 10¹⁵M_⊙. It shows a prominent substructure at ∼50 arcmin. This is the residual of a merger with two smaller groups. The violent merger with the main halo had an impact parameter of almost zero and the gas of one of the smallest systems has been completely stripped away. The collision direction is witnessed by a prominent bow shock that is visible as a sharp SZ edge immediately to the west of the merged objects. A third cluster, connected by a filament, is located ahead of the shock outside the field of view. The snapshot shown corresponds to z ∼ 0, thus we cannot clearly see the future of this object. However, by looking at its past motion we presume that it will most likely produce a second major merger.

To show that these test cases may be representative of real cluster configurations, in Figure 2 we compare the y parameter of two of the simulations with actual Planck measurements derived from a Modified Inter Linear Combination Algorithm (MILCA; Hurier et al. 2013; Planck Collaboration 2015b). Figure 2 demonstrates, in particular, that displacing the connected cluster at z = 0.04 yields an angular separation between the two cluster peaks that is comparable to the separation of both components of the cluster pair A399-A401. Moreover, displacing the colliding cluster at z = 0.01 leads to a Comptonization parameter decrement that matches its true analog that is measured across the most prominent shock front in Coma (see also Planck Collaboration 2013a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To emulate HFI observations of our simulated clusters, we choose a few characteristic regions of the sky in terms of thermal SZ S/N and predict the 2D signal to be registered in each HFI channel by adding our tSZ maps to the expected specific intensities of the CMB and Galactic dust components. For each HFI channel and sky region, frequency maps I(k, l, ν) are convolved with the HFI beam and added to a noise realization that matches the noise variance and power spectrum registered to the cluster scales in the raw data. The amplitude of the CMB SED is extracted from the Planck CMB map that is reconstructed using the Spectral Matching Independent Component Analysis (SMICA) component separation algorithm (Planck Collaboration 2014d), while the amplitude of the Galactic dust SED accommodates the dust optical depth mapped at 353 GHz as a result of the Planck all-sky model of thermal dust (Planck Collaboration 2014c). For simplicity, a constant spectral index, β = 1.8, is assumed for the dust SED, consistent with the first all-sky modeling from Planck and IRAS data (Planck Collaboration 2011b). An example of mock HFI frequency maps of the colliding cluster system is shown for (l = 270°, b = −30°) in Figure 3, while a few properties of the Galactic dust and instrumental noise included in the 350 GHz maps are summarized in Table 1 for each simulated sky position. The dust intensity in particular mostly depends on the Galactic latitude, and is shown to dominate the positive side of the thermal SZ signal (f_ν > 217 GHz).

Figures 4–6 present the restoration of simulated SZ maps for a soft thresholding,

$$\begin{eqnarray}&&\bar{{\mathcal{W}}}={\rm{sgn}}({\mathcal{W}}){\left(| {\mathcal{W}}| -\lambda \right)}_{+},\end{eqnarray} \tag{ 10 }$$

of the wavelet and curvelet coefficients at λ = 1.5σ, after 3 Van Cittert iterations performed with α = 0.25 (see Equation (9)). For each map, the restored signal, ${A}_{\mathrm{SZ}},$ is shown together with the S/N_2d:

$$\begin{eqnarray}&&{{\rm{S}}/{\rm{N}}}_{2{\rm{d}}}=10\mathrm{log}\displaystyle \frac{\parallel {A}_{\mathrm{SZ}}\parallel }{\parallel {A}_{\mathrm{SZ}}-\widehat{{A}_{\mathrm{SZ}}}\parallel }.\end{eqnarray} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Effect of the instrumental noise variance.

Figure 7 shows the tSZ maps of our test clusters positioned at low redshift, in a couple of sky regions corresponding to the lowest and highest variance of the instrumental noise in Table 1. In this figure, all characteristic features with a typical Compton parameter of y ≥ 10⁻⁵ are bright enough to be restored. More specifically, both components of the accreting cluster are recovered for each sky region, though the S/N_2d decrement makes it difficult to recover the true shape of the tenuous filament connecting the two clusters at (l = 0°, b = −15°). The extended bow shock escaping from the colliding clusters is also detectable for each sky region, together with the filament connecting these clusters to their third companion.

Effect of the instrumental noise variance and dust intensity.

In Figure 4, the tSZ signal of the connected cluster displaced at z = 0.04 is mapped for the sky positions detailed in Table 1. The cluster system is restored with an S/N_2d larger than 14 for each sky region. In particular, the position of the cluster peaks and the orientation of the connecting filament are preserved in any case. The image S/N_2d turns out to be roughly correlated with the instrumental noise variance regardless of the Galactic latitude, which suggests that it essentially depends on statistics.

Effect of the angular distance.

A few tSZ maps of our test clusters that were displaced at higher redshifts are finally shown on Figure 5. The redshift values of 0.02, 0.04, and 0.07 correspond to two successive decrements of the angular distances by a factor of two. Separated by two PSF radii, both components of the connected and accreting clusters are well-resolved, even at these redshifts. The shock front is also clearly visible at z = 0.02, while its companion filament is detected up to z = 0.04.

In the previous section we showed that ICM anisotropies are accurately restored by our algorithm, provided that they are spatially resolved in each HFI frequency map and located in cluster regions where the Compton parameter, y, exceeds 10⁻⁵. As shown in Figure 4, the situation is more critical for y values that are lower than 10⁻⁵, since the restoration depends on statistics and dust emissivity for a given region of the sky. In practice, a severe detection threshold (3σ or more) would allow us to isolate ICM features down to y values of 5 × 10⁻⁶ in most of the sky regions, in particular when they are extended enough to be analyzed and averaged out on several PSF sizes.

Though not intended to be considered absolute limits for each specific target, these values provide us with a coarse assessment of the potential of the presented algorithm in the 2015 release of Planck data. To predict a number of appropriate targets, we selected a sample of galaxy clusters in the Planck catalog of SZ sources whose r₅₀₀ radius exceeds the PSF size of each HFI frequency map, and whose detection thresholds' S/Ns are higher than 5 in the sense of the Multifrequency Matched Filter detection algorithm detailed in Melin et al. (2006). Figure 7 exhibits the radial profiles of the y parameter expected toward these clusters if one assumes a pressure structure that follows the average profile measured for representative samples of the X-ray cluster population (Arnaud et al. 2010). The cluster radii, r₅₀₀, and integrated-pressure, ${{\rm{Y}}}_{5{r}_{500}},$ follow the scaling relation proposed by the Planck collaboration to constrain cluster masses from SZ measurements (Planck Collaboration 2014e). The radially averaged Compton parameter exceeds 10⁻⁵ for 57 clusters at a distance from the center that exceeds the PSF size of each HFI frequency map. Lowering our selection criterion, at variance, to a y parameter of 7.5 × 10⁻⁶ would yield 81 clusters for which ICM anisotropies would be detectable beyond the PSF size of each HFI frequency map.

Beyond the brightest central regions of clusters, matter filaments also connect the component of cluster systems, and might become detectable at larger radii than expected within isolated, spherically symmetric clusters. A typical sample of targets for such studies is composed of 18 pairs of clusters from the Planck catalog, detected with an S/N threshold that is higher 5, showing an angular separation that is lower than 1° and a redshift interval that is lower than 0.01.