© ESO 2019

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) 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₅₀₀ <  r <  R₂₀₀ ¹, 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 ν₀ = 217 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 I 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 XI 2014; 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; hereafter B15) 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 RXC J1825. Finally, we summarise our results and outline future perspectives in Sect. 5. Throughout the paper, we assume a ΛCDM cosmological model with Ω_Λ = 0.7, Ω_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, 2019; Ghirardini et al. 2018, 2019; Ettori 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 XXIX (2014), 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) 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 θ₅₀₀, that is the angle subtended by the R₅₀₀ 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 × 256 pixel maps centred on each target cluster from the raw full-sky maps provided in the latest Planck data release² (see Planck Collaboration III 2019, for details on HFI). Specifically, we retrieved the temperature maps, the associated variance maps and the jackknife maps (see Planck Collaboration VIII 2016) 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 ≈4.3°, corresponding to sufficiently large radii (i.e. cluster-centric distances r >  3R₅₀₀) 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 XI (2014), 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 Comptonization parameter estimated from the maps obtained with the MILCA algorithm (Hurier et al. 2013) exceed the threshold value of 10⁻⁵. 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 ν. Referring to the kth pixel in the maps, the residuals can be simply written as

(1)

where D_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

(2)

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(ν), plus the instrumental noise, η.

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 X 2016). The spectral energy density of the tSZ effect is treated in the non-relativistic case, so that (e.g. Birkinshaw 1999)

(3)

where h_p is the Planck constant, k_B is the Boltzmann constant, and T_CMB = (2.725 ± 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):

(4)

where the integral is calculated along the line of sight (l.o.s.), being σ_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 Eq. (1) can be rewritten in terms of their wavelet transform as (see e.g. Mallat 2008)

(5)

where j and n give the dilation and translation of the wavelet basis function Ψ, respectively, being 2^j the wavelet scale. In our case, the wavelet basis function is a B₃ spline, and Φ is the dual scaling function of Ψ at the scale j₀. This scale corresponds to the approximation level – also called last smooth – which encodes signal information at the lowest resolution. The wavelet coefficients in Eq. (5) are given by

(6)

(7)

for the approximation and detail levels, respectively. In order to ensure normalisation and positivity when weighting the residuals, the wavelet kernel Ψ is split in its positive and negative parts, Ψ₊ and Ψ₋ (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_ν and N_pix, respectively) of such weighted residuals:

(8)

where 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:

(9)

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 scales³. 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 B₃ spline wavelet transform. Subsequently, we applied a soft thresholding at 1σ 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

(10)

where the dimensionless constant factors f₁ and q₁/q₂ refer to the relative contribution from the coldest component at temperature T₁ and the hottest component at temperature T₂. The β₁ and β₂ parameters give the slopes of the two different power laws, while B(ν; T₁) and B(ν; T₂) are the corresponding Planck functions describing the black body spectra. In order to get the best-fit parameters of the model described by Eq. (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₅₀₀. The only spatially-variable parameter is the temperature T₂, 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 Eq. (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 scales⁴, 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, ā_j0,n (see also Eq. (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₀. 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 XXVIII 2014; Planck Collaboration XXVI 2016, 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 Eq. (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 Eq. (1), and we calculated their wavelet coefficients as follows (see also Eq. (7) for a comparison):

(11)

(the same holds for the approximation level). We rewrite the parameter as s = s̃ ± Δs in Eq. (11), to highlight its spatial variation, Δ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, ℬ(ν), 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.

Fig. 1. Effect of beam-weighting of wavelet coefficients. a: step-like spatial variation of s (solid black line). The filled black area represents its convolution with the instrumental beam (dashed red line). b: effective spatial variation of s. The blue and red lines represent the absolute value of the positive and negative parts of the wavelet kernel, respectively. The envelopes of the shaded areas coincide with the convolution of the instrumental beam with the support of the positive and negative parts of the kernel.

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_tot sets of HFI maps. More specifically, we followed the steps detailed below.

Firstly, we generated N_tot = 100 Monte Carlo realisations of the noise at each HFI frequency, η_u(ν), with u = 1, …, N_tot referring to the uth extraction. We chose the value of N_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_HFI(ν), 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σ. 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, η_u(ν). Lastly, the mock data have been obtained as the summation of the denoised maps, D_HFI, den(ν), and the uth noise realisation, that is

(12)

We used such N_tot synthetic data sets as input to the imaging algorithm to obtain the following vectors of maps: , , . The standard deviations of all the bootstrap sets, namely std(s_tSZ), std(s_td), and std(s_CMB), represent our error estimates for each component. In the following, we label as σ_y = std(s_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/σ_y, for a given minimum value of σ_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 ≳ R₅₀₀ from the cluster centre. The average minimum values of the Compton y-parameter we can recover from the faintest resolved structures are 2.0 × 10⁻⁶ and 4.2 × 10⁻⁶, with a minimum significance of 3σ_y and 5σ_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σ detection is y = 5 × 10⁻⁶. 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).

Fig. 2. Compton y-parameter of X-COP clusters from our spectral imaging. The solid and dashed cyan circles on top of the images mark the R500 and R200 radii, respectively. The colour scale is logarithmic and contours are logarithmically spaced by a factor of , starting from 2 × 10−6. The J2000 Equatorial coordinates of the cluster centres are reported in the bottom right corner of each map.

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, 2019; Eckert 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σ.

Figure 3 shows the maps of the Compton y-parameter from the old and new versions of the algorithm (see Figs. 3a and b, 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 ≈ 6 × 10⁻⁶. 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.

Fig. 3. Region centred on cluster A2319. a: Compton y-parameter reconstructed with Van Cittert deconvolution. b: Compton y-parameter reconstructed with the wavelet coefficient-wise deconvolution. The colour scale is logarithmic and contours are logarithmically spaced by a factor of . c: raw HFI map at 857 GHz.

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 α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 ≲ 7.5 × 10⁻⁶. 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 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_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⁻⁶. 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 ≳ 2R₅₀₀ 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 ≈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.

Fig. 4. Compton y-parameter profiles of cluster A2319. The solid black curve of a is superimposed as a dashed red line to the cuts in b for comparison purposes (see text for details).

4.2. The cases of A2029 and RXCJ1825

Clusters A2029 and RXC J1825 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 Int. VIII (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σ_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/σ_y in Fig. 6a, this elongated excess of signal is significant to better than 5σ_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₂₀₀, 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 Int. VIII (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.

Fig. 5. X-ray and tSZ signals of clusters A2029 and RXCJ2815. The maps show the vignetting-corrected and background-subtracted surface brightness in the X-ray energy band 0.5–2.5 keV. The superimposed contours represent the Compton y-parameter from our algorithm. The solid and dashed cyan circles on top of the images are drawn at R500 and R200, respectively. The shaded circles in the bottom left corner of the maps represent a 5 arcmin beam. a: cluster A2029 mapped by ROSAT/PSPC. The cyan cross indicates the position of the neighbour cluster A2033. b: cluster RXC J1825 mapped by XMM-Newton. The cyan cross marks the position of cluster CIZA J1824, while the green arrow identifies a significant elongation which may be associated to stripped gas from a nearby group of galaxies.

Fig. 6. Maps of y/σy from a bootstrap run on clusters A2029 and RXC J1825. The overlapped contours represent the tSZ signal.

4.2.2. RXC J1825

In Fig. 5b, we show the X-ray surface brightness of RXC J1825 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 (hereafter CIZAJ 1824), 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σ_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₅₀₀, that contributes at a level of y ≈ 5 × 10⁻⁵, 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 σ_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 σ_y/y within a 10 arcmin-wide band, for the clusters A2029 and RXC J1825. 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 ≃ 2R₅₀₀ for both clusters.

Fig. 7. Profiles of σy/y of clusters A2029 and RXC J1825. The purple dashed lines correspond to an effective S/N of 3 (σy/y = 1/3, see text for details).

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₅₀₀. 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₅₀₀, 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.

Acknowledgments

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.

References

All Tables

All Figures

Fig. 1. Effect of beam-weighting of wavelet coefficients. a: step-like spatial variation of s (solid black line). The filled black area represents its convolution with the instrumental beam (dashed red line). b: effective spatial variation of s. The blue and red lines represent the absolute value of the positive and negative parts of the wavelet kernel, respectively. The envelopes of the shaded areas coincide with the convolution of the instrumental beam with the support of the positive and negative parts of the kernel.

In the text

	Fig. 2. Compton y-parameter of X-COP clusters from our spectral imaging. The solid and dashed cyan circles on top of the images mark the R500 and R200 radii, respectively. The colour scale is logarithmic and contours are logarithmically spaced by a factor of , starting from 2 × 10−6. The J2000 Equatorial coordinates of the cluster centres are reported in the bottom right corner of each map.
In the text

	Fig. 3. Region centred on cluster A2319. a: Compton y-parameter reconstructed with Van Cittert deconvolution. b: Compton y-parameter reconstructed with the wavelet coefficient-wise deconvolution. The colour scale is logarithmic and contours are logarithmically spaced by a factor of . c: raw HFI map at 857 GHz.
In the text

	Fig. 4. Compton y-parameter profiles of cluster A2319. The solid black curve of a is superimposed as a dashed red line to the cuts in b for comparison purposes (see text for details).
In the text

Fig. 5. X-ray and tSZ signals of clusters A2029 and RXCJ2815. The maps show the vignetting-corrected and background-subtracted surface brightness in the X-ray energy band 0.5–2.5 keV. The superimposed contours represent the Compton y-parameter from our algorithm. The solid and dashed cyan circles on top of the images are drawn at R500 and R200, respectively. The shaded circles in the bottom left corner of the maps represent a 5 arcmin beam. a: cluster A2029 mapped by ROSAT/PSPC. The cyan cross indicates the position of the neighbour cluster A2033. b: cluster RXC J1825 mapped by XMM-Newton. The cyan cross marks the position of cluster CIZA J1824, while the green arrow identifies a significant elongation which may be associated to stripped gas from a nearby group of galaxies.

In the text

	Fig. 6. Maps of y/σy from a bootstrap run on clusters A2029 and RXC J1825. The overlapped contours represent the tSZ signal.
In the text

	Fig. 7. Profiles of σy/y of clusters A2029 and RXC J1825. The purple dashed lines correspond to an effective S/N of 3 (σy/y = 1/3, see text for details).
In the text