Spectral imaging of the thermal Sunyaev-Zel'dovich effect in X-COP galaxy clusters: method and validation

TL;DR

This paper introduces a new method for creating detailed maps of the thermal Sunyaev-Zel'dovich effect in galaxy clusters, revealing complex structures in cluster outskirts with validation against X-ray data.

Contribution

The paper presents a novel locally parametric algorithm for sparse component separation in thermal SZ imaging, validated on Planck data for multiple galaxy clusters.

Findings

Detection of anisotropic features and filamentary structures in cluster outskirts.

Validation of SZ maps through comparison with X-ray observations.

Robustness of the method in low signal-to-noise regimes.

Abstract

The imaging of galaxy clusters through the Sunyaev--Zel'dovich effect is a valuable tool to probe the thermal pressure of the intra-cluster gas, especially in the outermost regions where X-ray observations suffer from photon statistics. For the first time, we produce maps of the Comptonization parameter by applying a locally parametric algorithm for sparse component separation to the latest frequency maps released by Planck. The algorithm takes into account properties of real cluster data through the two-component modelling of the spectral energy density of thermal dust, and the masking of bright point sources. Its robustness has been improved in the low signal-to-noise regime, thanks to the implementation of a deconvolution of Planck beams in the chi-square minimisation of each wavelet coefficient. We applied this procedure to twelve low-redshift galaxy clusters detected by Planck with the highest signal-to-noise ratio, considered in the XMM Cluster Oustkirts Project (X-COP). Our images show the presence of anisotropic features, such as small-scale blobs and filamentary substructures that are located in the outskirts of a number of clusters in the sample. The significance of their detection is established via a bootstrap-based procedure we propose here for the first time. In particular, we present a qualitative comparison with X-ray data for two interesting systems, namely A2029 and RXCJ1825. Our results are in agreement with the features detected in the outskirts of the clusters in the two bands.

Tables

Table 1. Table 1: Basic data of X-COP galaxy clusters  (taken from Ghirardini et al. 2019 ) . Galactic coordinates are taken from the NASA/IPAC Extragalactic Database ( https://ned.ipac.caltech.edu ).

Cluster name	Planck S/N	$z$	$M_{500}$ ($\times {10}^{14}$${\text{M}}_{\odot }$)	$R_{500}$ (kpc)	${\theta }_{500}$ (arcmin)	$l$ (deg)	$b$ (deg)
A2319	30.8	0.0557	7.31	1346	20.8	75.70	13.51
A3266	27.0	0.0589	8.80	1430	21.0	272.13	-40.13
A2142	21.3	0.0909	8.95	1424	14.1	44.22	48.68
A2255	19.4	0.0809	5.26	1196	13.1	93.97	34.95
A2029	19.3	0.0766	8.65	1414	16.3	6.44	50.53
A85	16.9	0.0555	5.65	1235	19.2	115.23	-72.03
A3158	17.2	0.0597	4.26	1123	16.3	265.05	-48.93
A1795	15.0	0.0622	4.63	1153	16.1	33.82	77.18
A644	13.9	0.0704	5.66	1230	15.3	229.93	15.29
A1644	13.2	0.0473	3.48	1054	19.0	304.88	45.45
RXCJ1825‡	13.4	0.0650	4.08	1105	14.8	58.31	18.54
ZW1215‡	12.8	0.0766	7.66	1358	15.7	282.50	65.19

Table 2. Table 2: Frequencies and corresponding beams of Planck HFI (from Planck Collaboration et al. 2016a ) .

Frequency (GHz)	Beam FWHM (arcmin)
100	9.69
143	7.30
217	5.02
353	4.94
545	4.83
857	4.64

Table 3. Table 3: Astrophysical offsets subtracted from each raw HFI map (see text for details).

Frequency (GHz)	Offset
100	$1.47\text{×}{10}^{-5}$ K${}_{\text{CMB}}$
143	$2.32\text{×}{10}^{-5}$ K${}_{\text{CMB}}$
217	$7.02\text{×}{10}^{-5}$ K${}_{\text{CMB}}$
353	$4.12\text{×}{10}^{-4}$ K${}_{\text{CMB}}$
545	$3.41\text{×}{10}^{-1}$ MJy sr-1
857	$5.84\text{×}{10}^{-1}$ MJy sr-1

Equations

$$\text{res}(\nu ,k;s)=D_{\text{HFI}}(\nu ,k)-M(\nu ,k;s),$$

$$M(\nu ,k;s)=\sum _i^{N_s}{f}_i(\nu )s_i(k)+\eta (\nu ,k),$$

$$f_{\text{tSZ}}(\nu )=\frac{h_p\nu }{k_B{T}_{\text{CMB}}}\times \coth \left(\frac{h_p\nu }{2k_B{T}_{\text{CMB}}}\right)-4,$$

$$y=\frac{{\sigma }_T}{m_e{c}^2}{\int }_{\text{los}}p(l)dl,$$

$$\begin{array}{rl}{\text{res}}_{\Psi }(\nu ,k;s)& =\sum _n^{N_{\text{pix}}}{\bar{a}}_{j_0,n}(\nu ;s){\Phi }_{j_0,n}(k)+\\ & +\sum _{j=j_0}^{N_{\text{scales}}}\sum _n^{N_{\text{pix}}}{a}_{j,n}(\nu ;s){\Psi }_{j,n}(k),\end{array}$$

$${\bar{a}}_{j_0,n}(\nu ;s)=\sum _m^{N_{\text{pix}}}\text{res}(\nu ,m;s){\Phi }_{j_0,n}^{\ast }(m);$$

$$a_{j,n}(\nu ;s)=\sum _m^{N_{\text{pix}}}\text{res}(\nu ,m;s){\Psi }_{j,n}^{\ast }(m),$$

$${\chi }_{{\Psi }_{\pm }}^2=\sum _{\nu }^{N_{\nu }}\sum _k^{N_{\text{pix}}}\frac{{\text{res}}_{{\Psi }_{\pm }}^2(\nu ,k;s)}{{\sigma }_{\text{HFI}}^2(\nu ,k)},$$

$$\hat{s}=\frac{1}{2}\left[\underset{s}{\mathrm{argmin}}({\chi }_{{\Psi }_{+}}^2)-\underset{s}{\mathrm{argmin}}({\chi }_{{\Psi }_{-}}^2)\right].$$

$$\begin{array}{rl}{f}_{\text{td}}(\nu )& =f_1\frac{q_1}{q_2}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_1}B(\nu ;T_1)+\\ & +(1-f_1){\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_2}B(\nu ;T_2),\end{array}$$

$$\begin{array}{rl}{a}_{j,n}(\nu ;\Delta s)& =\sum _m^{N_{\text{pix}}}\text{res}(\nu ,m;\Delta s){\Psi }_{j,n}^{\ast }(m)\\ & =\sum _m^{N_{\text{pix}}}\{D_{\text{HFI}}(\nu ,m)-\mathcal{B}(\nu )\otimes [HM(\nu ,m;\tilde{s}+\Delta s)+\\ & +(1-H)M(\nu ,m;\tilde{s}-\Delta s)]\}{\Psi }_{j,n}^{\ast }(m),\end{array}$$

$${D_{\text{HFI}}}_u(\nu )=D_{\text{HFI,den}}(\nu )+{\eta }_u(\nu ).$$

Full text

1 Introduction

Clusters of galaxies, which are the largest gravitationally bound structures in the universe, can be studied at millimetre wavelengths through observations of the Sunyaev–Zel’dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972). Indeed, as the photons of the cosmic microwave background (CMB) radiation propagate through the hot and ionised intra-cluster medium (ICM hereafter) of galaxy clusters, they get diffused via inverse Compton scattering by the free electrons. The consequence of this interaction is that photon energy is re-distributed; therefore, the observed black body spectrum of the CMB gets distorted. The spectral distortion produced by the random thermal motion of the electrons is known as the thermal SZ (tSZ) effect. A remarkable property is its linear dependence on both electron density and temperature, making it a direct probe of the thermal pressure of the ICM integrated along the observer line of sight. For this reason, the tSZ effect is extremely useful when investigating the outskirts of galaxy clusters, which host complex physical processes associated to the continuous accretion of matter (see e.g. Walker et al. 2019, for a recent review). In this respect, the XMM Cluster Outskirts Project (X-COP, see Eckert et al. 2017), proposes combining the tSZ effect with X-ray data at radii $R_{500}<r<R_{200}$111Whenever the subscript $\Delta=200,500$ is used, we refer to the mass enclosed by a spherical volume of radius $R_{\Delta}$, whose density is $\Delta$ times the critical density of the universe at a given redshift $z$, $\rho_{c}(z)=3H(z)^{2}/(8\pi G)$, being $H(z)$ the Hubble parameter, and $G$ the gravitational constant, which are barely accessible to X-ray spectroscopy because of significant contamination by unresolved astrophysical or instrumental backgrounds.

Another important feature of the tSZ effect, which makes it unique among other observables at microwave frequencies, is its peculiar spectral signature. Indeed, the observed change of the CMB brightness manifests as a decrement at frequencies below the zero-frequency $\nu_{0}=$217\text{,}\mathrm{GHz}$$, and as an increment at higher frequencies. In the non-relativistic approximation (Kompaneets 1957), both the zero-frequency and the overall spectral signature do not depend on cluster physical parameters. Therefore, thanks to this characteristic, it is possible to disentangle the tSZ signal from other astrophysical sources using multi-frequency data.

The separation of different components contributing to an astrophysical signal is an issue that has been extensively addressed in the literature with a variety of techniques, and for diverse purposes. Among the most used methods, there are refinements of the internal linear combination (ILC), which was first proposed by Eriksen et al. (2004). These methods assume there is no correlation between the different physical components, and they are non-parametric, which means that they do not rely on any particular model. The contribution from each astrophysical source is estimated through a simple linear combination of the total signal at all the available frequencies. The weighting coefficients for each frequency band are chosen in such a way that they provide unit response with respect to the component of interest, while minimising the variance of the reconstructed map. The increasing amount of microwave data at high angular resolution, which can be of a few arcmin or arcsec with currently operating instruments (see e.g. Mroczkowski et al. 2019, for a recent review), and the unprecedented frequency coverage provided by the Planck satellite (Planck Collaboration et al. 2011), motivated the development of several modified ILC algorithms tailored to the extraction of the tSZ effect. Relevant examples in this respect are provided by: Remazeilles et al. (2011), Hurier et al. (2013) (who propose the modified internal linear combination algorithm: MILCA), and Hurier et al. (2017). Despite their success, ILC methods for the extraction of the tSZ effect suffer from some drawbacks (see e.g. Bobin et al. 2013, for a discussion). For instance, results can be significantly biased if the different foreground sources are correlated, as well as if they exhibit a spatially variable spectral energy density. Moreover, ILC methods perform better in the case of Gaussian foreground fluctuations, which is not a realistic hypothesis for the tSZ effect and the emission from thermal dust.

In addition, the combination of heterogeneous data sets from ground-based and space-based experiments has been explored through the proposal of more sophisticated techniques, with the aim of exploiting the advantages coming from observations at different angular resolutions and sensitivities. For instance, Remazeilles et al. (2013) propose a method which implements needlet decomposition, while in the recent work by Aghanim et al. (2019) the authors apply suitable matched multi-filters (see e.g. Herranz et al. 2002; Melin et al. 2006; Tarrío et al. 2016, 2018) to data from Planck and from the Atacama Cosmology Telescope (Swetz et al. 2011).

Besides non-parametric algorithms, methods based on a prior modelling of the components and on likelihood maximisation have also been developed (see e.g. Feroz et al. 2009; Khatri 2016). In particular, parametric techniques are robust, especially for the mapping of Galactic thermal dust (see for instance Planck Collaboration et al. 2014a; Meisner & Finkbeiner 2015). On the other hand, algorithms based on sparse representations allow one to exploit some properties of the foreground sources, such as their morphological independence or positive normalisation.

In Bourdin et al. (2015) (B15 hereafter) the authors propose to link sparse imaging to a parametric estimate of foregrounds, via the introduction of a spatially-weighted likelihood function through wavelets and curvelets. Wavelets are basis functions characterised by a scale (or equivalently by its inverse, that is the resolution) and by a position parameter. Wavelet transforms provide sparse representations of signals characterised by a smooth distribution, apart from localised discontinuities. Since the ICM seems to hold such properties, wavelet transforms have been widely used in cluster astronomy, both for X-ray (Slezak et al. 1994; Vikhlinin et al. 1997; Starck & Pierre 1998) and tSZ imaging (Pierpaoli et al. 2005; Pires et al. 2006). Another relevant application for millimetre astronomy is the foreground cleaning of CMB maps proposed by Bobin et al. (2013, 2016), which relies on sparse priors to spatially separate the components. Curvelet transforms are a higher-dimensional generalisation of wavelet transforms. Indeed, curvelet bases are parametrised by two additional quantities, namely direction and elongation (Candès & Donoho 2000). This property makes them particularly suitable to get sparse representations of images featuring edges and anisotropies (Mallat 2008), as in the case of shocks and filamentary structures in the outskirts of galaxy clusters.

For the reasons above, in B15 the authors employ wavelet decomposition in the extraction of the CMB signal and curvelet decomposition to separate the anisotropic signal from the thermal dust and the tSZ effect. The algorithm is tested using mock Planck frequency maps from the set of high-resolution hydrodynamical simulations of galaxy clusters presented in Planelles et al. (2014) and Rasia et al. (2014).

We applied for the first time the spectral imaging algorithm by B15 to the frequency maps from the latest data released by the Planck Collaboration in 2018. Furthermore, we improved the algorithm by implementing an efficient removal of astrophysical contaminations, such as the emission from thermal dust and bright point sources. Most importantly, we enhanced its stability by means of a new deconvolution technique applied to the wavelet coefficients, proposed for the first time in this work. We show the maps of the tSZ effect for a set of twelve massive and nearby clusters of galaxies selected for X-COP. We focus in particular on two interesting cases of possibly interacting cluster couples, and we compare our images to ancillary X-ray data. The significance of the signal is established through a new boostrap-based procedure for the estimate of the tSZ error maps, which we also used to show the overall enhancement of the algorithm, compared to its previous implementation.

The paper is structured as follows. In Sect. 2 we briefly introduce the cluster data set. Section 3 illustrates the details of the algorithm we used for producing maps of the tSZ effect, focussing on its new aspects. In Sect. 4 we present and discuss our results for the X-COP sample. In particular, we make a comparison with the implementation of the algorithm presented in B15, and we detail the cases of cluster A2029 and RXCJ1825. Finally, we summarise our results and outline future perspectives in Sect. 5. Throughout the paper, we assume a $\Lambda$CDM cosmological model with $\Omega_{\Lambda}=0.7$, $\Omega_{m}=0.3$ and $h=0.7$.

2 Data set

We applied our imaging procedure to the sample of galaxy clusters studied in the XMM Cluster Outskirts Project (X-COP). The goal of the project is to unveil the thermodynamic and dynamical properties of the ICM in the outermost cluster regions, combining data in the X-ray and millimetre bands (see Tchernin et al. 2016; Eckert et al. 2017; Ghirardini et al. 2018; Ettori et al. 2019; Eckert et al. 2019; Ghirardini et al. 2019).

The main physical properties of the X-COP clusters are summarised in Table 1.

The sample consists of twelve massive objects located at low to intermediate redshift, namely in the range $0.04<z<0.10$, observed in X-ray by XMM-Newton (Turner et al. 2001; Strüder et al. 2001) and at millimetre wavelengths by Planck. As discussed in Planck Collaboration et al. (2014b), the six bands observed with the High Frequency Instrument (HFI) centred at 100, 143, 217, 353, 545, and 857 GHz are those which provide the highest signal-to-noise ratio (S/N hereafter) for the detection of galaxy clusters. In particular, the clusters in the X-COP sample have been detected with S/N $>12$. Their characteristic angular size $\theta_{500}$, that is the angle subtended by the $R_{500}$ radius of each cluster, is larger than 10 arcmin, so that they are well-resolved by Planck. A list of the HFI beams at all frequencies is reported in Table 2.

We extracted $256\times 256$ pixel maps centred on each target cluster from the raw full-sky maps provided in the latest Planck data release222Planck data are publicly available at: https://pla.esac.esa.int (see Planck Collaboration et al. 2018, for details on HFI). Specifically, we retrieved the temperature maps, the associated variance maps and the jackknife maps (see Planck Collaboration et al. 2016b) at each frequency band. In particular, the latter are used in our procedure to estimate the error, as explained in detail in Sect. 3.3. Sky patches are obtained through a gnomonic re-projection using the HEALPix package (Górski et al. 2005). Each pixel is 1 arcmin in size, so the maps extend over $\approx 4.3^{\circ}$, corresponding to sufficiently large radii (i.e. cluster-centric distances $r>3R_{500}$) to allow the detection of substructures in cluster outskirts.

The signal in the raw temperature maps is corrected at all HFI frequencies by the astrophysical offsets listed in Table 3, which account for the cosmic infrared background and the tSZ background.

The offset values are computed as described in Planck Collaboration et al. (2014a), from the correlation between HFI data at 857 GHz and data delivered by the Leiden/Argentine/Bonn (LAB) radio survey of Galactic HI (Kalberla et al. 2005). The pixels in the sky that are used to compute this correlation are selected according to the values of both the column density and the velocity of HI clouds. In order to remove the contamination from galaxy clusters, and to account for calibration at high frequencies, we also added two more selection criteria. One is based on masking pixels where the the Comptonization parameter estimated from the maps obtained with the MILCA algorithm (Hurier et al. 2013) exceed the threshold value of ${10}^{-5}$. The second one accounts for a correction of the CMB using SMICA maps (Cardoso et al. 2008) at frequencies between 100 GHz and 353 GHz as templates.

3 Imaging of the tSZ effect

To map the tSZ signal, we used a parametric imaging algorithm featuring wavelet and curvelet decomposition. This procedure was originally proposed in B15 and tested only on mock HFI observations of interacting clusters. The first application to real HFI data motivated substantial adaptations and improvements of its original implementation. In particular, we used the spectral energy density of a double grey body to model the emission from thermal dust, and we corrected for residual contamination from both dust and bright point sources. Another remarkable improvement concerns with the stability of the algorithm, which we obtained through a new deconvolution procedure proposed here for the first time. In this Section we first provide a synthetic description of the original algorithm, then we focus on the aforementioned enhancements.

3.1 The spectral imaging algorithm

The novelty of our imaging method consists of combining the parametric approach to component separation and sparse representations. This is achieved by finding the parameters which minimise a modified chi-square, accounting for the wavelet transform of the residuals between the data and the model map at each frequency $\nu$. Referring to the $k$-th pixel in the maps, the residuals can be simply written as

[TABLE]

where $D_{\textup{HFI}}$ represent the HFI data. The source component maps $s$, which are the parameters we want to estimate, enter the residuals via the model map, $M$. The latter can be written explicitly as

[TABLE]

that is the sum over the temperature anisotropies produced by the $N_{s}$ physical sources to separate, $s_{i}$, multiplied by the corresponding spectral energy density $f_{i}(\nu)$, plus the instrumental noise, $\eta$.

The dominating astrophysical sources in the frequency range covered by Planck HFI, apart from the CMB, are the tSZ effect and the Galactic thermal dust (see for instance Planck Collaboration et al. 2016c). The spectral energy density of the tSZ effect is treated in the non-relativistic case, so that (e.g. Birkinshaw 1999)

[TABLE]

where $h_{p}$ is the Planck constant, $k_{B}$ is the Boltzmann constant, and $T_{\textup{CMB}}=(2.725\pm 0.001)$ K (Mather et al. 1999) is the CMB temperature, which is constant at all frequencies. The amplitude of the tSZ signal is given by the Compton $y$-parameter (see e.g. Rephaeli 1995):

[TABLE]

where the integral is calculated along the line of sight (los), being $\sigma_{T}$ the Thomson cross section, $m_{e}$ the electron mass, $c$ the speed of light, and $p(l)$ the thermal pressure of the ICM. The spectral modelling of Galactic thermal dust (td subscript) is presented in detail in Sect. 3.2, since it represents a major change in the algorithm.

The residuals in Eqn. (1) can be rewritten in terms of their wavelet transform as (see e.g. Mallat 2008)

[TABLE]

where $j$ and $n$ give the dilation and translation of the wavelet basis function $\Psi$, respectively, being $2^{j}$ the wavelet scale. In our case, the wavelet basis function is a B3 spline, and $\Phi$ is the dual scaling function of $\Psi$ at the scale $j_{0}$. This scale corresponds to the approximation level – also called last smooth – which encodes signal information at the lowest resolution. The wavelet coefficients in Eqn. (5) are given by

[TABLE]

for the approximation and detail levels, respectively. In order to ensure normalisation and positivity when weighting the residuals, the wavelet kernel $\Psi$ is split in its positive and negative parts, $\Psi_{+}$ and $\Psi_{-}$ (see fig. 1 of B15 for reference). This yields the minimisation of two separate chi-squares, which can be expressed as the sum across all the frequencies and pixels ($N_{\nu}$ and $N_{\textup{pix}}$, respectively) of such weighted residuals:

[TABLE]

where $\sigma^{2}_{\textup{HFI}}$ is the HFI variance map. Therefore, the estimate of the component maps is given by the half-sum of the results from the two separate minimisations:

[TABLE]

To maximise the efficiency in the recovery of signal anisotropies, we chose a curvelet basis to decompose the tSZ and dust signals (see e.g. Candès & Donoho 2000), and a wavelet basis for the CMB (see B15, for more details). The decomposition has been performed over four scales333The maximum number of scales that can be used for wavelet decomposition is linked to the number of pixels on each side of the image as: $N_{\textup{scales}}=\text{floor}[\log(N_{\textup{pix}})/\log 2]$ - 4. The curvelet transform has been computed following the procedure described in Starck et al. (2003), that is to say trough ridgelet transforms applied to the bidimensional wavelet bands of the B3 spline wavelet transform. Subsequently, we applied a soft thresholding at $1\sigma$ level to the curvelet and wavelet coefficients, in order to keep only the relevant features of the signal. Eventually, the final maps have been obtained by means of a suitable restoration operator, which re-combines the approximation coefficients and the detail coefficients together. Pixels plagued by contaminants, due for instance to a local bad modelling of thermal dust, have been automatically masked. Specifically, we discarded those pixels where chi-square exceeded a given threshold value, depending on the wavelet scale. At the same time, we imposed a condition of regularity on the error of the wavelet coefficients, in order to select regions characterised by a high S/N of the tSZ component.

3.2 Improvements and new features

The three main changes we implemented with respect to the version of the imaging algorithm presented in B15 are detailed in the following subsections. Two of them concern with a more realistic treatment of astrophysical contaminants. The third one consists of a new deconvolution technique proposed here for the first time.

3.2.1 Modelling of thermal dust

Following Meisner & Finkbeiner (2015) (see also Bourdin et al. 2017), we modelled the thermal dust as a double grey body, by assuming two populations of dust grains, instead of the idealised case of a single grey body spectrum. Indeed, the latter provides an accurate representation of the thermal emission from Galactic dust only at frequencies higher than 353 GHz. The spectral energy density we set for this component is therefore

[TABLE]

where the dimensionless constant factors $f_{1}$ and $q_{1}/q_{2}$ refer to the relative contribution from the coldest component at temperature $T_{1}$ and the hottest component at temperature $T_{2}$. The $\beta_{1}$ and $\beta_{2}$ parameters give the slopes of the two different power laws, while $B(\nu;T_{1})$ and $B(\nu;T_{2})$ are the corresponding Planck functions describing the black body spectra. In order to get the best-fit parameters of the model described by Eqn. (10), we calculated an independent fit to the dust component through a Monte Carlo Markov Chain sampling. To treat only the signal from the dust, we limited this fit to the pixels in the frequency maps that are located sufficiently far from the cluster, at radial distances from the centre larger than $5R_{500}$. The only spatially-variable parameter is the temperature $T_{2}$, which is fixed a priori to the value determined by a joint fit to IRAS and Planck data, as detailed in Meisner & Finkbeiner (2015). From this fit we obtain maps of the dust component at all frequencies, which we plug in the model maps of Eqn. (2).

3.2.2 Removal of contamination from thermal dust and point sources

Despite the two grey body model provides an accurate description of the thermal dust component, our maps of the Compton $y$-parameter turned out to be affected by contamination either from diffuse signal from thermal dust on large spatial scales444We notice, however, a significant reduction of this contamination when comparing maps from Planck 2015 data release to those from the 2018 data release, or from compact point sources, mostly at the frequency of 857 GHz. To solve these issues, we proceeded as described in the following.

We reconstructed the tSZ images without accounting for the contribution from the 857 GHz frequency in the last approximation of its discrete wavelet transform, $\bar{a}_{j_{0},n}$ (see also Eqn. (5)). On the other hand, we took advantage of all the HFI frequency maps to compute the detail coefficients of the wavelet transform, $a_{j,n}$, supposedly insensitive to image gradients below the last approximation scale, $j_{0}$. This allowed us to suppress unphysical large-scale gradients in the approximation coefficients of some cluster images. These gradients could be possibly related to a residual contamination of the thermal dust signal by the large-scale variance of other astrophysical sources or of the instrumental noise.

To remove the contamination from point sources, we used the masks in the second Planck catalogue of compact sources (PCCS2) (see Planck Collaboration et al. 2014c, 2016d, for details). These objects are dusty galaxies and synchrotron emitting in the HFI range at high and low frequencies, respectively. Specifically, we combined the gnomonic projection of all the masks at the six HFI frequencies into a unique mask, after smoothing them to the common resolution of 18 arcmin. To remove the contribution from such sources, we multiplied this mask by the wavelet kernel used in the chi-square minimisation. We applied this procedure only to the clusters in our sample which are significantly affected by bright point sources at far infrared frequencies, namely A3266, A85, and ZW1215.

3.2.3 Deconvolution procedure

As can be seen from Eqn. (8), the weighted chi-square mixes the different angular resolutions of the data. Therefore, a deconvolution must be implemented to restore the signal. The original version of the algorithm features a Van Cittert deconvolution (Van Cittert 1931), regularised via successive projections in the significant curvelet domain. In other words, the resolution of the model map is adjusted to the resolution of the data map through iterative corrections of the curvelet coefficients. However, several tests revealed that this technique hampered the reliability of the signal reconstructed in low signal regimes, because of the amplification of pixel-sized diverging artefacts across the iterations.

To avoid this problem, we implemented a new “wavelet coefficient-wise” deconvolution of the Planck beams. Specifically, we deconvolved the amplitude of the single wavelet coefficients, instead of deconvolving the final image. To this end, we took the full expression of the residuals as in Eqn. (1), and we calculated their wavelet coefficients as follows (see also Eqn. (7) for a comparison):

[TABLE]

(the same holds for the approximation level). We rewrite the parameter as $s=\tilde{s}\pm\Delta s$ in Eqn. (11), to highlight its spatial variation, $\Delta s$. Instead of splitting the wavelet kernel in two components, we used the Heaviside step function, $H$, that coincides with the positive part of the kernel, while $(1-H)$ coincides with its negative part. This operation is equivalent to taking the absolute value of the wavelet function. To match the angular resolutions of the data and the model, we convolved the model maps at each frequency with the Planck beam, $\mathcal{B}(\nu)$, which can be approximated with a Gaussian having the FWHM listed in Table 2. The effect of such a convolution is that of correlating the contributions from the negative and the positive components of the kernel (see also Fig. 1). Thus, the final estimate of the parameters takes into account the results from the single decompositions. Besides allowing us to deconvolve spatial variations of the model maps, the weighting by the absolute value of the wavelet kernel within a single chi-square minimisation increases the support of the kernel itself, as can be seen in Fig. 1b. This is equivalent to enlarging the size of the sky region where we compute the chi-square, thus enhancing the S/N.

3.3 Procedure for error assessment

To associate a statistical error to the estimate of the component maps, we used a bootstrap procedure that allowed us to simulate $N_{\textup{tot}}$ sets of HFI maps. More specifically, we followed the steps detailed below.

Firstly, we generated $N_{\textup{tot}}=100$ Monte Carlo realisations of the noise at each HFI frequency, $\eta_{u}(\nu)$, with $u=1,\dots,N_{\textup{tot}}$ referring to the $u$-th extraction. We chose the value of $N_{\textup{tot}}$ as a trade-off between the computational time needed to perform wavelet and curvelet transforms and the statistical significance. We assumed the instrumental noise to be Gaussian and spatially correlated, and we imposed the noise maps to have the same power spectrum as the jackknife maps. Subsequently, the HFI raw data at each frequency, $D_{\textup{HFI}}(\nu)$, have been denoised through a simple wavelet-based procedure, by calculating the wavelet transform on three scales, and soft-thresholding the coefficients at $1.5\sigma$. This procedure relies on spatially-variant thresholds for each wavelet band, which have been preliminarily inferred by computing the variance of the coefficients across wavelet transforms of the noise maps, $\eta_{u}(\nu)$. Lastly, the mock data have been obtained as the summation of the denoised maps, $D_{\textup{HFI,den}}(\nu)$, and the $u$-th noise realisation, that is

[TABLE]

We used such $N_{\textup{tot}}$ synthetic data sets as input to the imaging algorithm to obtain the following vectors of maps: $\bm{s}_{\textup{tSZ}}=(s_{\textup{tSZ}}^{1},\dots,s_{\textup{tSZ}}^{N_{\textup{tot}}})$, $\bm{s}_{\textup{td}}=(s_{\textup{td}}^{1},\dots,s_{\textup{td}}^{N_{\textup{tot}}})$, $\bm{s}_{\textup{CMB}}=(s_{\textup{CMB}}^{1},\dots,s_{\textup{CMB}}^{N_{\textup{tot}}})$. The standard deviations of all the bootstrap sets, namely $\text{std}(\bm{s}_{\textup{tSZ}})$, $\text{std}(\bm{s}_{\textup{td}})$, and $\text{std}(\bm{s}_{\textup{CMB}})$, represent our error estimates for each component. In the following, we label as $\sigma_{y}=\text{std}(\bm{s}_{\textup{tSZ}})$ the error for the tSZ signal only. The level of significance of the anisotropies (such as blobs and filaments) we find in the tSZ images of each cluster can be assessed by computing the ratio $y/\sigma_{y}$, for a given minimum value of $\sigma_{y}$, which represents an ‘effective’ S/N.

4 Results and discussion

We show in Fig. 2 the zoomed maps of the Compton $y$-parameter we obtained from the application of the procedure described in Sect. 3 to the HFI maps of the X-COP clusters.

In the majority of cases, the algorithm highlights the presence of blobs or extended filamentary structures at radial distances $r\gtrsim R_{500}$ from the cluster centre. The average minimum values of the Compton $y$-parameter we can recover from the faintest resolved structures are $2.0\text{\times}{10}^{-6}$ and $4.2\text{\times}{10}^{-6}$, with a minimum significance of $3\sigma_{y}$ and $5\sigma_{y}$, respectively. This represents an improvement of a factor of 2.5 in sensitivity with respect to the version of the algorithm presented in B15, where the minimum signal for a $3\sigma$ detection is $y=$5\text{\times}{10}^{-6}$$. It is important to stress that, thanks to the deconvolution, the algorithm is capable to recovering the elongated structure of the signal in the cluster central regions (see e.g. the cases of A2142 and A1644).

A particular cluster case which is worth mentioning is that of A2319, which is known to be a complex merging system that behaves as an outlier with respect to the other clusters in the sample (see Ghirardini et al. 2018; Eckert et al. 2019; Ghirardini et al. 2019). The tSZ map of this object shows a number of blobs located in the virial region, which may be due to either substructures or clumpy ICM patches. The detailed study of these features and of their impact on the thermodynamic properties of this system will be addressed in a forthcoming dedicated paper.

In the following, we first perform a comparison with the results from the version of the algorithm presented in B15 for cluster A2319. Subsequently, we briefly discuss a couple of two possibly interacting systems, namely A2029 and RXCJ1825.

4.1 Comparison with the previous version of the algorithm

We illustrate here the improvements introduced by the adaptation of our procedure to real cluster data, and by the wavelet coefficient-wise deconvolution. In particular, in order to highlight the impact of removing the highest frequency from the approximation coefficients (see Sect. 3.2.2), we adopted a double grey body model for thermal dust also in the procedure by B15. To perform a consistent comparison, we applied the two versions of the algorithm to the same HFI maps. The wavelet decomposition has been performed here over three scales, and coefficients have been thresholded at $1\sigma$.

Fig. 3 shows the maps of the Compton $y$-parameter from the old and new versions of the algorithm (see Figs. 3a and 3b, respectively), in the case of cluster A2319.

This object is particularly suitable to show the effect of removing the 857 GHz channel from the last smooth coefficients. Indeed, in Fig. 3a it is possible to see a large-scale gradient to the left-hand side of the image at $y\approx$6\text{\times}{10}^{-6}$$. Such a contaminant is likely due to a residual signal from thermal dust, as can be seen from the comparison with the HFI map at 857 GHz, shown in Fig. 3c. This residual signal is significantly removed with the new adaptations, as demonstrated in the map in Fig. 3b.

Most notably, Fig. 3 highlights how the two deconvolution techniques have a different impact on the final maps. Indeed, the result in Fig. 3a is produced with the Van Cittert deconvolution using three iterations, and with the convergence parameter $\alpha$ set to 0.25 (see eq. (9) of B15). It can be seen that this map shows a number of pixel-sized outliers, which plague regions where the signal is $y\lesssim$7.5\text{\times}{10}^{-6}$$. Such artefacts are due to a diverging amplification of the curvelet coefficients, as a consequence of the iterative nature of this deconvolution algorithm. The map in Fig. 3b shows instead that the result from the new wavelet coefficient-wise deconvolution is cleaner, and the mildly ellipsoidal shape of the signal in the centre is still well recovered. It is worth noting that both techniques are sensitive to the basis functions used to project the deconvolved signal. To be specific, Van Cittert method is a regularisation applied to curvelet coefficients, therefore it is particularly suitable to highlight ellipsoidal and elongated features. The wavelet coefficient-wise deconvolution, on the contrary, is applied to wavelet coefficients, thus it provides a better localisation and detection of peaks. This can be seen in in Fig. 3b, where the signal in the cluster centre is recovered more efficiently, with respect to Fig. 3a. Another important element that improves peak localisation is the larger kernel support which enters the chi-square minimisation in the new deconvolution (see Sect. 3.2.3 and Fig. 1). Indeed, this feature allows the reconstruction of the signal with a higher S/N, since in this case the wavelet coefficients are computed and deconvolved on a larger sky region.

To further illustrate the differences between the two deconvolutions, we show as solid black lines in Fig. 4 the profiles of the Compton $y$-parameter. Each profile has been extracted as the average vertical cut within a 10 arcmin-wide band, passing through the centre of the corresponding tSZ image in Fig. 3. The superimposed light grey lines represent the same profiles from each of the $N_{\textup{tot}}=100$ different maps of the Comptonization parameter obtained with the bootstrap procedure (see Sect. 3.3). It can be seen from Fig. 4b that the non-iterative nature of the new deconvolution improves significantly the stability in low-signal regimes, allowing one to reach a minimum level of $y$ of the order of ${10}^{-6}$. On the other hand, the diverging artefacts which hampered the reliability of the detection of substructures in the outskirts with the Van Cittert deconvolution, are clearly visible in several realisations of the signal in Fig. 4a. Moreover, the dispersion of the bootstrap profiles at radii $r\gtrsim 2R_{500}$ is on average 50 per cent lower with the new procedure, resulting into a lower error in the reconstructed signal in these regions.

The different characterisation of the central peak is also clear from Fig. 4b, where the value recovered with the Van Cittert deconvolution is $\approx 23$ per cent lower with respect to the value obtained with the new one based on wavelet coefficients. Such a difference is due to the joint effect of the larger S/N provided by the absolute value of the wavelet kernel, and of the limited number of iterations (set to 3) used in the Van Cittert deconvolution. The latter has been arbitrarily chosen as a trade-off between a reliable recovery of the peak, and the avoidance of divergences. Therefore, it has a non-negligible impact on the final results; on the contrary, the new deconvolution does not rely on any arbitrary parameter. From the comparisons reported above, we can conclude that the new version of the algorithm does provide a more reliable and stable reconstruction of the Comptonization parameter.

4.2 The cases of A2029 and RXCJ1825

Clusters A2029 and RXCJ1825 are two interesting targets for testing the algorithm. Among the X-COP clusters, they have been detected by Planck with intermediate (19.3) and low (13.4) S/N, respectively. A2029 has been widely studied in X-rays (see e.g. Lewis et al. 2002; Clarke et al. 2004; Bourdin & Mazzotta 2008; Walker et al. 2012); on the contrary, RXCJ1825 has been poorly investigated in this band since its discovery (Ebeling et al. 2002). Both these clusters may be interacting with two known neighbouring systems, as suggested in Planck Collaboration et al. (2013). Nevertheless, given the low significance of the data, no further analysis on the tSZ signal from possible connecting filaments has been explored in their work.

4.2.1 A2029

We show in Fig. 5a the contours from our map of the tSZ effect, superimposed to the X-ray surface brightness of A2029. We report the data from the ROSAT satellite, which allows the detection of the X-ray emission at larger radii than those probed by Chandra or XMM-Newton, thanks to its low particle background (Trümper 1982; Vikhlinin et al. 1999). The raw data have been denoised via the adaptive smoothing technique (Eckert et al. 2012).

It can be seen that the tSZ signal in the central region matches fairly well with the X-ray surface brightness. Thanks to the deconvolution, we are able to detect the tSZ emission corresponding to the neighbouring cluster A2033 with a significance of $8\sigma_{y}$, which together with A2029 belongs to a small supercluster (Einasto et al. 2001). The X-ray and the tSZ peaks of A2033 show an offset which is, in any case, smaller than the best resolution of 5 arcmin provided by Planck. Our tSZ map clearly highlights an elongated projected structure that connects the two clusters. As testified by the image showing $y/\sigma_{y}$ in Fig. 6a, this elongated excess of signal is significant to better than 5$\sigma_{y}$. X-ray images do also show such an elongated morphology in the peripheral regions to the north-east, pointing to A2033. This suggested a possible ongoing merger between the two objects (see e.g. Eckert et al. 2012; Walker et al. 2012). However, recent analyses of the reconstructed density field through gravitational lensing, indicate that this signal is likely due to the gas in the overlapping outskirts of the two clusters at $R_{200}$, rather than to a filament connecting them (Gonzalez et al. 2018). Our tSZ imaging does represent an improvement with respect to the map shown in fig. 1 of Planck Collaboration et al. (2013), where the emission from A2033 is barely detected, as well as the signal between the two clusters. However, we cannot favour any of the two hypothesised processes. On the other hand, our maps may help in discriminating among possible scenarios proposed to explain the observed excess, for instance in future works combining microwave and X-ray data to model the three-dimensional thermodynamic properties of the ICM in each component of the system.

4.2.2 RXCJ1825

In Fig. 5b, we show the X-ray surface brightness of RXCJ1825 from XMM-Newton, featuring wavelet denoising. The tSZ contours from our algorithm are superimposed to the image. As in the case of A2029, the tSZ emission follows pretty well the X-ray surface brightness. In particular, we notice that the peaks of the signal in the two bands coincide within 1 arcmin. The tSZ signal shows a significant elongation towards the neighbouring cluster CIZA J1824.1+3029 (CIZAJ1824 hereafter), even if the peak corresponding to this object is not clearly visible. A recent analysis based on the kinematics of their member galaxies, suggests that RXCJ1825 and CIZAJ1824 are in a pre-merger state (Girardi et al. 2019). Thus, also in this case, our maps may shed a new light on the hypothesis of a possible interaction in future analyses. Interestingly, our tSZ map shows a $14\sigma_{y}$-significant elongation to the south-west of the cluster, highlighted with a green arrow in Fig. 5b, which agrees with the excess in the X-ray surface brightness detected with XMM-Newton. Such an emission may be due to gas stripped from the past interaction between RXCJ1825 and a small, disrupted group of galaxies, which have been detected at the same redshift in the optical band (Clavico et al. 2019). This scenario is also supported by the recent finding of a radio halo extending in the same direction of the X-ray elongation (Botteon et al. 2019). Another feature we notice in the signal is the presence of a third structure located within $R_{500}$, that contributes at a level of $y\approx$5\text{\times}{10}^{-5}$$, and which is not detected in X-rays. Nevertheless, the significance of this detection is 35 per cent lower with respect to the signal in the centre, as demonstrated by the map in Fig. 6b, constructed from a bootstrap run. Therefore, its detection is likely due to some localised irregularity in the instrumental noise in the raw HFI data.

In order to further quantify the significance of the signal reconstructed with our algorithm as a function of the radius, we computed the profiles of the ratio $\sigma_{y}/y$ that is the inverse of the effective S/N. Specifically, we show in Fig. 7 the average of the vertical cuts passing through the centre of the maps of $\sigma_{y}/y$ within a 10 arcmin-wide band, for the clusters A2029 and RXCJ1825. It can be seen that our algorithm provides a reconstruction of the signal with an effective S/N $>3$ (marked with a dashed purple line in the plots), up to radii $r\simeq 2R_{500}$ for both clusters.

5 Conclusions

We presented maps of the Comptonization parameter of the twelve nearby, massive galaxy clusters constituting the X-COP sample (Eckert et al. 2017). We used an improved version of the spectral imaging algorithm for parametric component separation proposed in B15, featuring wavelet and curvelet decomposition. In particular, the enhancements we introduced are: (i) a double grey body spectral energy density to model the emission from thermal dust; (ii) the cleaning of residual contamination from dust and compact point sources; (iii) a new wavelet coefficient-wise deconvolution via beam-weighting in the calculation of the wavelet transform. We also illustrated a new method to estimate the error on the reconstructed tSZ signal based on bootstrap extractions of Planck noise maps. With this method, we showed it is possible to detect the signal from diffuse filaments and small substructures located beyond $R_{500}$. Indeed, we detected interesting features in the outskirts of the majority of the objects in the sample under study. In particular, we highlighted diffuse bridges in significant excess with respect to the background level, which connect the brightest objects of the cluster systems A2029 and RXCJ1825. These findings are consistent with ancillary surface brightness maps in the X-ray band. When compared to its original implementation, our new deconvolution technique provided more stability, an enhanced localisation of the central peak and a more precise reconstruction of the signal in the outskirts. In particular, the test case of cluster A2319 showed a reduction of contamination from outliers at radii $r>2R_{500}$, and about 50 per cent lower values of the bootstrap error in the same region.

We plan to use the tSZ maps presented in this work to isolate and mask the signal from substructures located in cluster outskirts, which will be the topic of a forthcoming paper. Indeed, this signal may bias-high the reconstructed radial profiles of thermal pressure in the virial region, with non-negligible consequences on the parameters of the profile, and on estimates of cluster masses relying on the assumption of hydrostatic equilibrium.

Cluster A2319, which is the one detected by Planck with the highest S/N, will be the subject of a dedicated upcoming work, in which we will highlight in detail the interesting features characterising the tSZ map that we were able to reconstruct by employing sparse representations.

Acknowledgements.

The authors thank the anonymous referee for useful comments on the manuscript. They acknowledge the Planck Collaboration for the extraordinary legacy they delivered to the scientific community. They also acknowledge the use of the image analysis routines of the Interactive Sparse astronomical data Analysis Package (ISAP) developed in the CosmoStat laboratory at CEA Saclay. This work makes also an extensive use of python, particularly of the astropy and matplotlib libraries. The authors thank Sabrina De Grandi for helpful suggestions. A.S.B. is grateful to Rocco D’Agostino for useful discussions, and she acknowledges funding from Sapienza Università di Roma - Progetti per Avvio alla Ricerca Anno 2018, prot. AR11816430D9FBEA. A.S.B., H.B. and P.M. acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.0., from ASI Grant 2016-24-H.0., and from “Tor Vergata” Grant “Mission: Sustainability” EnClOS (E81I18000130005). S.E. acknowledges financial contribution from the contracts ASI 2015-046-R.0 and ASI-INAF n.2017-14-H.0. M.G. is supported by the Lyman Spitzer Jr. Fellowship (Princeton University) and by NASA Chandra grants GO7-18121X and GO8-19104X.