© ESO 2018

1 Introduction

Galaxy clusters can be detected from observations at different bands of the electromagnetic spectrum, each of them probing a different component of the cluster. In optical observations we can see the individual galaxies inside the cluster, which contribute to around 1% of the total cluster mass. Clusters are identified in these images as overdensities of galaxies. Clusters can also be detected in X-ray observations where they appear as bright sources with extended emission. In these images we observe the emission of the hot gas of the intracluster medium (ICM), which accounts for 10–15% of the cluster mass. Over the last decade, this gas has also begun to be detected thanks to the characteristic spectral distortion it produces on the cosmic microwave background (CMB) due to Compton scattering of the CMB photons by the ICM electrons. This effect is known as the Sunyaev–Zeldovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972).

State-of-the-art galaxy cluster detection techniques usually rely on the analysis of single-survey observations. However, combining information from different surveys at different wavelengths can potentially improve the detection performance, allowing us to find more distant or less massive clusters. Although multiwavelength, multisurvey detection of clusters was theoretically conceived some years ago (Maturi 2007; Pace et al. 2008), it is a very complex task and, to date, has only been attempted in practice in the pilot study of Schuecker et al. (2004) on X-ray data from the ROSAT All-Sky Survey (RASS; Truemper 1993; Voges et al. 1999) and optical data from the Sloan Digital Sky Survey (SDSS; York et al. 2000). In our previous work (Tarrío et al. 2016), we proposed a new analysis tool based on matched multifrequency filters (MMF) for extracting clusters from SZ and X-ray maps. The method was based on the combination of the classical SZ MMF (Herranz et al. 2002; Melin et al. 2006, 2012) and an analogous single-frequency matched filter developed for X-ray maps. It was shown that combining these two complementary sources of information improved the signal-to-noise ratio (S/N) with respect to SZ-only or X-ray-only cluster extraction, and also provided correct photometry as long as the physical relation between X-ray and SZ emission of the clusters, namely the expected F_X∕Y₅₀₀ relation, was known. The filter was used as an extraction tool to estimate some properties of already detected clusters, but not to detect new clusters in a blind manner, since the position, the size, and the redshift of the clusters were assumed to be known.

In this paper we propose a blind detection method based on the X-ray–SZ filter studied in Tarrío et al. (2016). The goal is to adapt this filter to use it as a blind cluster detection tool, given that we do not know the position, the size, and the redshift of the clusters. As demonstrated in Tarrío et al. (2016), combining X-ray and SZ information increases the cluster S/N with respect to single-map extractions. This gain in S/N will translate directly into a higher detection probability for a given threshold in the S/N. We would also expect to obtain, in principle, a higher purity than when using the classical SZ MMF since the combined version will clean out objects whose emission is far from the expected F_X∕Y₅₀₀ relation. However, even if an object does not follow the expected relation, it could still pass the detection threshold if it has a very strong signal in the X-ray band. This is similar to the strong infrared emissions that were detected with the classical SZ MMF used by the Planck Collaboration VIII (2011); Planck Collaboration XXIX (2014); Planck Collaboration XXVII (2016) even though their spectra did not fit the expected SZ spectrum. This may be the case of some non-cluster X-ray sources, such as active galactic nuclei (AGNs). Since the X-ray filter was designed to be easily compatible with the classical SZ MMF and it is not specifically optimized for X-ray cluster detection, it does not consider the extent of the sources, as other X-ray cluster detection techniques do (Böhringer et al. 2000; Vikhlinin et al. 1998; Pacaud et al. 2006; Ebeling & Wiedenmann 1993; Scharf et al. 1997). As a result, when we add the X-ray information, we will also add false detections produced by non-cluster X-ray sources (mainly AGNs). Therefore, as already noted in Tarrío et al. (2016), the main challenge to be solved when using the proposed X-ray–SZ MMF for blind detection is to obtain a high purity.

The proposed method is applied to observations from the ROSAT All-Sky Survey (RASS) and the Planck survey, the latest full-sky X-ray and SZ surveys available to date. Nevertheless, the proposed joint detection technique is general and also applicable to other surveys, including those from future missions such as e-ROSITA (Merloni et al. 2012), a four-year X-ray survey which is planned to start in 2018 and which will be much deeper than RASS.

The structure of the paper is as follows. Section 2 presents RASS and Planck observations. In Sect. 3 we describe the joint X-ray–SZ detection algorithm. Section 4 presents an evaluation of its performanceusing RASS and Planck maps by comparing its results with other cluster catalogues in the SPT region. Finally, we conclude the paper and discuss ongoing and future research directions in Sect. 5.

Throughout, we adopt a flat Λ CDM cosmological model with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_M = 1 −Ω_Λ = 0.3. We define R₅₀₀ as the radius within which the average density of the cluster is 500 times the critical density of the universe, θ₅₀₀ as the corresponding angular radius, and M₅₀₀ as the mass enclosed within R₅₀₀.

2 Description of the observations

Although the joint algorithm proposed in this paper is a general technique that can be applied, in principle, to any X-ray and SZ surveys of the sky, we have tested it in this paper using all-sky maps from the Planck and RASS surveys. In this section, we briefly describe these observations.

2.1 RASS data description

The ROSAT All-Sky Survey (RASS) is, to date, the only full-sky X-ray survey conducted with an X-ray telescope (Truemper 1993; Voges et al. 1999). The survey data release¹ contains 1378 individual RASS fields in three different bands: TOTAL (0.1–2.4 keV), HARD (0.5–2.0 keV), and SOFT (0.1–0.4 keV). Each field covers an area of 6.4 deg× 6.4 deg (512 × 512 pixels) and has a resolution of 0.75 arcmin/pixel.

In this paper, we use an X-ray all-sky HEALPix map that we built from the HARD band information. This map has a resolution of 0.86 arcmin/pixel (HEALPix resolution closest to the RASS resolution). The details of its construction can be found in Appendix B of Tarrío et al. (2016). Although our approach could optimally include all X-ray bands, we chose to use only the HARD band. This provides a better S/N for the clusters because the SOFT band is dominated by the diffuse X-ray background of the local bubble. This is a common choice in cluster detection surveys based on RASS data, such as REFLEX (Böhringer et al. 2001, 2013).

2.2 Planck data description

Planck is the most recent space mission that was launched to measure the anisotropy of the CMB. It observed the sky in nine frequency bands from 30 to 857 GHz with high sensitivity and angular resolution. The Low Frequency Instrument (LFI) covers the 30, 44, and 70 GHz bands, while the High Frequency Instrument (HFI) covers the 100, 143, 217, 353, 545, and 857 GHz bands.

In this paper, we use only the six temperature channel maps of HFI, which are the same channels used by the Planck Collaboration to produce their cluster catalogues (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). In particular, we used the latest version of these maps; their description can be found in Planck Collaboration VIII (2016). The published full resolution maps have a resolution of 1.72 arcmin/pixel. However, to make them directly compatible with the all-sky X-ray map mentioned above, we up-sampled them to a resolution of 0.86 arcmin/pixel by zero-padding in the spherical harmonics domain.

3 Joint detection of galaxy clusters on X-ray and SZ maps

In this section, we describe the proposed algorithm for the blind detection of galaxy clusters using X-ray and sub-mm maps. The algorithm isbased on the X-ray–SZ extraction method proposed in Tarrío et al. (2016), which extracts the characteristics of a cluster given its known position, size, and redshift. In this paper, we adapted this extraction method to perform a blind detection of clusters, i.e. to discover clusters in the maps without knowing their positions, sizes, or redshifts.

3.1 X-ray–SZ MMF

Let us first briefly recall the joint X-ray–SZ extraction method proposed in Tarrío et al. (2016). This method is based on a matched filter approach and was designed to be compatible with the SZ MMF known as MMF3, described by Melin et al. (2012) and used by the Planck Collaboration VIII (2011); Planck Collaboration XXIX (2014); Planck Collaboration XXVII (2016) to construct their SZ cluster catalogues.

The main idea of the joint extraction algorithm is to consider the X-ray map as an additional SZ map at a given frequency and to introduce it, together with the other SZ maps, in the classical SZ-MMF. In order to do so, the X-ray map needs to be converted into an equivalent SZ map at a reference frequency ν_ref, leveraging the expected F_X∕Y₅₀₀ relation. The details of this conversion are described in Appendix B of Tarrío et al. (2016). Once the X-ray map is expressed in the same units as the SZ maps, we can apply the classical MMF to the complete set of maps (the original N_ν SZ maps obtained at sub-mm frequencies , and an additional SZ map at the reference frequency ν_ref obtained from the X-ray map). The reference frequency ν_ref is just a fiducial value with no effect on the extraction algorithm. In our case, we took ν_ref = 1000 GHz.

The X-ray–SZ MMF presented in Tarrío et al. (2016) is given, in Fourier space², by (1)

with (2)

where is a (N_ν + 1) × 1 column vector whose ith component will filter the map at observation frequency ν_i; is, approximately, the background noise variance after filtering; and P(k) is the noise power spectrum, a (N_ν + 1) × (N_ν + 1) matrix whose ij component is given by , where N_i (k) is the noise map at observation frequency ν_i, which includes instrumental noise and astrophysical sources different from the cluster signal (extragalactic point sources, diffuse Galactic emission, and the primary CMB anisotropy for the SZ maps and X-ray background for the X-ray map). We note that the properties of the noise are not the same in different sky regions; therefore, P (k) has to be calculated locally at each position. Finally, is a (N_ν + 1) × 1 column vector defined as (3)

where j(ν_i) is the SZ spectral function at frequency ν_i and and are the convolutions of the cluster 2D spatial profiles ( for the SZ profile and for the X-ray profile) with the point spread function (PSF) of the instruments at the different frequencies. The 2D cluster profiles and are normalized so that their central value is 1. Finally, the constant C is a geometrical factor that accounts for the different shapes of the SZ and X-ray 3D profiles; it is defined in Eq. (25) of Tarrío et al. (2016) as the ratio of the integrated fluxes of the normalized SZ and X-ray 3D profiles up to R₅₀₀. As we can see, the filter is determined by the shape of the cluster signal and by the power spectrum of the noise, hence the name of matched filter.

This matched filter approach relies on the knowledge of the normalized cluster profile. This profile is not known in practice, so we need to approximate it by the theoretical profile that best represents the clusters we want to detect. As in Tarrío et al. (2016), we assume the generalized Navarro–Frenk–White (GNFW) profile (Nagai et al. 2007) given by (4)

with parameters given by (5)

and (6)

for the components corresponding to the original SZ maps and the additional X-ray map, respectively. These parameters come from assuming the 3D pressure profile of Arnaud et al. (2010) and the average gas density profile from Piffaretti et al. (2011), respectively. We note that x = θ∕θ₅₀₀ represents here the 3D distance to the centre of the cluster in θ₅₀₀ units, and θ₅₀₀ relates to the characteristic cluster scale θ_s through the concentration parameter c₅₀₀ (θ_s = θ₅₀₀∕c₅₀₀). The cluster profile is then obtained by numerically integrating these 3D GNFW profiles along the line of sight.

Finally, this cluster profile needs to be convolved by the instrument beams. As in Tarrío et al. (2016), in this paper we use the six highest frequency Planck maps and the X-ray maps of the ROSAT All-Sky Survey. Therefore, we use the same instrument beams as in Tarrío et al. (2016), namely, a Gaussian PSF for the SZ components, with FWHM depending on the frequency, as shown in Table 6 of Planck Collaboration VIII (2016), and a PSF for the X-ray component that was estimated numerically by stacking observations of X-ray point sources from the Bright Source Catalogue (Voges et al. 1999).

Figure 1 shows two examples of the radial profiles of the filter for θ_s = 1 and θ_s = 30 arcmin where we can see the spectral and the spatial weighting introduced by the filter. The filter was computed at a random position, with galactic coordinates 260.356°, −20.332°. We note that the last component of the filter corresponds to the X-ray band, and that its relative amplitude depends on the chosen reference frequency ν_ref.

Fig. 1 Examples of the matched filter for θs = 1 arcmin (a) and θs = 30 arcmin (b). The curves give the radial profiles of the filters, which are symmetric because we have chosen a symmetric cluster template. The filter is normalized so that its maximum amplitude is equal to 1.

3.2 Blind procedure

In Tarrío et al. (2016), the filter described above was proposed as an extraction tool, to estimate the flux of a cluster once we know that there is a cluster at a given position x₀, and if we know its size θ_s and redshift (necessary to convert the X-ray into an equivalent SZ map through the F_X ∕Y₅₀₀ relation). In this section we describe how this method is adapted to become a detection tool.

Given a set of N_ν + 1 maps (SZ + X-ray) of a given region of the sky , where M_ref is the X-ray map already converted into SZ units, the first step to detect new clusters consists of filtering the maps with the filter defined in Eq. (1) as (7)

In this way, we obtain a ŷ-map (filtered map) and a S/N map (ŷ (x)/) with the same size as the observed maps.

We note that to calculate the filter , we first need to estimate the noise power spectrum P(k). This is done in practice from the X-ray and SZ images themselves, assuming that they contain mostly noise. In the case of the X-ray images, this assumption may not be true due to bright X-ray sources with strong signals. Therefore, to minimize this effect, we masked some regions of the X-ray images for the calculation of P (k). In particular, we masked the areas defined in Table 1 of Böhringer et al. (2001) corresponding to the Large and Small Magellanic Clouds, and we also masked the X-ray sources of the ROSAT bright source catalogue (Voges et al. 1999) that have a count rate >0.3 counts/s.

Since the size of the clusters is unknown, we repeat the filtering process using a set of N_s filters with different sizes, covering the expected range of radii. In our case, we vary θ₅₀₀ from 0.94 to 35.31 arcmin, in N_s = 32 steps equally spaced in logarithmic scale. For each size, we obtain a filtered map and a S/N map. The clusters are then detected as peaks in these S/N maps, down to a given threshold.

Finally, for the conversion of the X-ray map into an equivalent SZ map we need to assume a F_X ∕Y₅₀₀ relation, which depends on the redshift. As studied in Tarrío et al. (2016), the assumed F_X ∕Y₅₀₀ relation does not have a big impact on the estimated S/N, which makes the detection robust against possible errors in the assumed relation. For this reason, we have fixed the redshift to a reference value of z_ref = 0.8 and assumed therelation found by the Planck Collaboration I (2012): (8)

where the K-correction can be obtained by interpolating in Table 2 of Tarrío et al. (2016).

To implement the detection procedure in practice, we proceed in two phases.

• 1.

  Producing a preliminary list of candidates. In this first phase, we project the all-sky maps into 504 small 10° × 10° tangential patches, as was done in MMF3 (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). Each patch is filtered by the X-ray–SZ filter (Eq. (1)) using N_s = 32 different sizes, which produces N_s S/N maps. Then, we construct a list with the peaks in these maps that are above a specified S/N threshold q. The procedure is as follows:

• (a)

  We look for the highest peak among all the N_s S/N maps;

• (b)

  If it is above the specified threshold q, we include its position in a preliminary candidate list and mask it in the N_s S/N maps. The size of the mask is defined as the radius at which the value of the filtered template is 1% of its maximum value. The filtered template is the 2D cluster profile corresponding to the size at which the highest peak was found, convolved by the PSF, and filtered by the X-ray–SZ filter with the N_s = 32 different sizes. Thus, the size of the mask is different in each of the N_s S/N maps;

• (c)

  Then werepeat the search until there are no more peaks above the specified threshold.

Finally, we merge the 504 lists into a single preliminary all-sky list of candidates by merging peaks that are close to each other by <10 arcmin, as was done in MMF3 (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016).

• 2.

  Refiningthe list of candidates. In this second phase, we reanalyse each candidate in the preliminary candidate list. This second phase is necessary to better estimate the candidate properties and S/N, since the results from the first phase may not be accurate. This is especially true if the candidate is situated close to a border of the map, since the estimated noise in this case may not be representative of the noise around the candidate. For each candidate we use the following procedure:

• (a)

  We produce a set of N_ν + 1 (SZ + X-ray) 10° × 10° tangential maps centred at the candidate position;

• (b)

  We filter these maps with the different filter sizes, obtaining N_s S/N maps;

• (c)

  We estimate the S/N of the detection by selecting a small circular region around the centre in each of the N_s S/N maps and by searching for local maxima inside this volume. This is necessary because the position of the peak may have changed slightly when centring the tangential maps. Among all the local maxima, we select the one with highest S/N that is not on the border of the circles (to avoid tails of nearby objects);

• (d)

  If this S/N is above a specified threshold q, we add the detection, with its new position and corresponding size, to the final candidate list, otherwise we discard it.

3.3 Determination of the threshold on the joint S/N

An important point of the blind joint detection algorithm is the selection of the threshold q that is applied to the peaks found in the first and second phases. The goal of this threshold is to discard false detections (noise peaks) with a givenconfidence. This can be achieved by setting the probability that a detection is due to a random fluctuation to a sufficiently low value. In the MMF3 method, a fixed threshold is used under the assumption that the noise distribution is Gaussian, so that a fixed S/N threshold leads to a fixed number of noise peak detections. In the joint X-ray–SZ detection, the Gaussian assumption is no longer valid, as explained below, so the threshold must be selected differently.

The probability density function (PDF) of the S/N in the joint filtered maps depends on the noise properties of the observed maps. Due to the Poisson nature of the noise in the X-ray maps, the final PDF of the joint S/N is not Gaussian. Its shape depends on the exposure time of the X-ray map and also on the filter size. In particular, it becomes more long-tailed when the exposure time is low, especially for small filter sizes; in these cases, most of the pixels of the X-ray map contain zero photons, and just a few pixels contain one photon. As a consequence, the average background is very low and the S/N of the filtered map, defined as ŷ (x)/, at the few pixels with one photon can be easily quite high.

Since the PDF of the joint S/N will have different shapes in different regions of the sky, using a fixed S/N threshold to detect cluster candidates everywhere in the sky will produce a different number of false detections, e.g. more detections will appear in low exposure time regions due to single noise pixels with high S/N. To have an approximately constant number of false detections over the whole sky, we need to establish an adaptive threshold that depends on the noise characteristics of each region. Since the PDF of the joint S/N cannot be calculated analytically, we have determined this adaptive threshold numerically by means of Monte Carlo simulations.

In particular, we performed an experiment in which we simulated a set of N_ν + 1 = 7 maps emulating in a simple manner the noise properties of Planck and RASS maps.

• The RASS noise map was simulated as a homogeneous Poisson random field, characterized by a given mean value λ (in counts/pixel). This noise represents the instrumental noise and the astrophysical X-ray background (mainly due to diffuse galactic emission and non-resolved point sources). To express this map into X-ray flux units, we assumed an exposuretime of 400 s and a N_H of 2 × 10²⁰ cm⁻² (average values in the SPT region). We note that the simulation results obtained for these values can be converted to results that would be obtained for any other values of exposure time and N_H, as explained in Appendix A, so this choice does not have any implications. Finally, as was done with the real RASS maps, this X-ray flux map is converted into an equivalent SZ map at the reference frequency ν_ref by following the procedure detailed in Appendix B of Tarrío et al. (2016).

• The N_ν = 6 Planck noise maps were simulated as the sum of two independent components: primary anisotropies and white Gaussian noise. First, we used the Planck Sky Model (Delabrouille et al. 2013) to obtain a realization of the CMB for the N_ν = 6 Planck frequencies (100, 143, 217, 353, 545, 857 GHz). Second, for each frequency, we added zero-mean Gaussian random noise with a frequency-dependent variance. In particular, the variance at frequency ν was fixed to the value σ_ν = σ₂₁₇ * [1.66, 0.70, 1.00, 3.12, 19.50, 649.87]³, where σ₂₁₇ is the standard deviation of the Gaussian noise in the 217 GHz map. Therefore, the simulated Planck noise maps are characterized by this single parameter σ₂₁₇.

We repeated the experiment for different values of the mean Poisson level λ and the Gaussian noise level σ₂₁₇, and for each pair of values λ-σ₂₁₇ we used 450 different realizations of the noise maps. At each realization, we changed the Poisson and the Gaussian noises (maintaining their levels), as well as the CMB realization.

Each set of N_ν + 1 = 7 maps was then filtered using the proposed joint filter with the N_s = 32 different sizes, yielding the variance of the X-ray filtered map , the variance of the SZ filtered maps , and the S/N map for each filter size.

From these results, we calculated the average values of and corresponding to each noise level and filter size (see Tables 1 and 2). Finally, we established the joint S/N threshold q_J for a given λ-σ₂₁₇-θ_s triplet as the S/N value for which the fraction of pixels (considering the 450 realizations) with S/N > q_J does not exceed a given false alarm probability P_FA. We used the number of pixels with S/N > q_J as an approximation of the number of detections with S/N > q_J. Due to the iterative blind detection procedure, where each S/N peak is masked after detection (step 1b of the blind procedure described in Sect. 3.2), one detection spans more than 1 pixel. However, the approximation allows a much faster computation and it is also accurate enough, especially for small filter sizes, which is the regime where the Poisson noise peaks become more important and which we need to characterize better.

The value of P_FA serves to select the operational point of the detection method. The higher the P_FA, the lower the threshold q_J and the more candidates we keep, resulting in a catalogue with higher completeness and lower purity. On the contrary, if we want a very pure catalogue at the expense of being less complete, we can choose a small value of P_FA. Table 3 summarizes the S/N thresholds for each combination of noise and some filter sizes calculated for a false alarm rate of P_FA = 3.4 × 10⁻⁶, which corresponds to a cut at 4.5σ in a zero-mean Gaussian distribution.

For reasons of simplicity, we apply this adaptive threshold q_J after we have obtained the final candidate list from the second phase of the blind procedure. The threshold q to be applied in the first and second phases is established to a sufficiently low value so that it does not introduce any different selection effect, i.e. it does not discard any candidate above q_J. In our case, we selected q = 4 for the first and second phases, which is lower than any of the adaptive thresholds shown in Table 3. The adaptive threshold q_J is then used to discard noise detections in the following way.

• 1.

  For each detection in the final candidate list provided at the second phase of the joint blind algorithm, we save the joint S/N of the detection, the corresponding filter size θ_s, the standard deviation of the X-ray filtered map , and the standard deviation of the SZ filtered maps ;

• 2.

  We then calculate the mean Gaussian level that corresponds to the measured by interpolating in the simulation results (Table 2);

• 3.

  Then, we calculate the mean Poisson level that corresponds to the measured . To this end, we need to take into account that the measured corresponds to a map with an exposure time and a N_H that are different from those used in the simulations (t_exp = 400 s and N_H = 2 × 10²⁰ cm⁻²). Therefore, we will first convert the measured into the value that we would have obtained with the values of exposure time and N_H used in the simulations, and then use Table 1 to interpolate the value of λ. The conversion from the measured into its simulation-equivalent counterpart is done using Eq. (A.6). A detailed description of this conversion can be found in Appendix A;

• 4.

  We choose the two simulated values of σ₂₁₇ that are closer to (one above: , one below: ) and the two simulated values of λ that are closer to (λ⁺ and λ⁻). Then, weselect the four simulations corresponding to the filter size θ_s and the four possible combinations of , , λ⁺ and λ⁻;

• 5.

  Then we calculate the S/N threshold for each of the four selected simulations. We cannot take the thresholds in Table 3 directly because they correspond to an X-ray map with t_exp = 400 s and N_H = 2 × 10²⁰ cm⁻² and the real map where we have detected the candidate will have, in general, different characteristics. To correct for this effect, we need to convert the S/N maps obtained in the simulations into equivalent S/N maps corresponding to the exposuretime and N_H of the real map. This is done using Eq. (A.9), as explained in Appendix A. Then, for each of the four selected simulations, we calculate the S/N threshold as the value q_Ji for which the fraction of pixels on the transformed S/N maps with S/N > q_Ji is at most P_FA;

• 6.

  Finally, the threshold q_J to be applied to our detection is obtained via a 2D interpolation between the four values of q_Ji. If the S/N of the detection is above this threshold (S/N > q_J), we keep it in the list since it is not likely to be a noise detection (with confidence 1-P_FA); otherwise, we discard it.

3.4 Catalogue preparation

The blind detection outputs an all-sky catalogue of joint X-ray–SZ detections that may still contain non-cluster objects or detections caused by noise or contamination. As demonstrated in Tarrío et al. (2016), adding the X-ray information to the SZ maps increases the cluster detection probability, allowing us to detect fainter or more distant clusters with respect to the catalogues constructed from purely SZ information. However, introducing the X-ray maps also increases the number of false detections, produced by non-cluster X-ray sources (mainly AGNs) and Poisson noise. Furthermore, the SZ observations contain regions contaminated with infrared emission that may also produce false detections. Thus, a main challenge of the proposed X-ray–SZ blind detection is to obtain a high purity. To this end, the catalogue produced by the blind detection method needs to be cleaned to discard detections in contaminated regions of the sky, in regions with poor statistics, or that correspond to non-cluster objects. In the rest of this section, we introduce two masking procedures to avoid detections in regions with infrared contamination (Sect. 3.4.1) or X-ray poor statistics (Sect. 3.4.2), and a method for discarding real detections corresponding to non-cluster X-ray objects (Sect. 3.4.3).

3.4.1 SZ mask

To avoid SZ contaminated regions, we follow the same procedure used to build the second Planck cluster catalogue (PSZ2) described in Planck Collaboration XXVII (2016), i.e. we discard all the detections inside the PSZ2 survey mask, we also reject detections within 5σ_beam of any SZ compact source of the second Planck catalogue of compact sources (PCCS2; Planck Collaboration XXVI 2016) with S/N > 10 in any of the six HFI Planck channels, and we remove 7 arcmin matches with the Planck cold-clump catalogue (C3PO), or with PCCS2 detections at both 545 GHz and 857 GHz to eliminate spurious infrared detections.

3.4.2 X-ray mask

In the regions of the sky where the X-ray exposure time is very low, the X-ray count-rate map contains few noise pixels with very high count-rate values (typically, one count divided by a very low exposure time) compared to adjacent pixels (with zero counts). These bright pixels may introduce false detections, which could be discarded using the adaptive threshold calculated by numerical simulations. However, since the amount of simulation time required to properly simulate these regions is significant, and given that the X-ray information provided by these regions is very limited, we decided to just set a threshold in the exposure time to mask the low-exposure regions. In the case of RASS, we decided to use a threshold of 100 s, which masks only 5.5% of the sky, i.e. 2300 deg², where 910 deg² have a non-zero exposure time and 1390 deg² have no RASS observations. The overlap of this masked area with the SPT footprint, where the proposed method will be evaluated (see Sect. 4), is 209 deg², 80 deg² with non-zero exposure time and 129 deg² with no RASS observations.

3.4.3 Classification to distinguish clusters from point sources

Some of the objects detected with the blind joint detection method correspond to point sources in the X-ray maps that coincide with a SZ noise peak. Although the estimated size can be used as a criterion to distinguish between a real cluster and a point source, it is difficult to distinguish between a cluster with a small apparent size and a point source.

Our ideal aim would be to recognize whether a detection is a real cluster or a point source given the parameters extracted during the filtering process. To check if this is possible, we cross-matched some joint detections with a list of known clusters and known point sources (see details below), and we labelled each of our matching candidates as belonging to one or the other class. Then, we characterized each sample of this labelled list using five features: the estimated S/N, size and flux of the blind joint detection, and the X-ray and SZ components of the S/N: (S/N)_XR and (S/N)_sz. Finally, in order to get an upper bound on the best classification accuracy that can be obtained with a linear classifier, we trained a support vector machine (SVM) classifier with this labelled list.

Using ten-fold cross-validation⁴, we obtained a classification accuracy of 88%, with 5% of the mis-classifications being point sources classified as clusters and 7% being clusters classified as point sources. We noticed that the parameter that plays the most important role in the classification is (S/N)_sz, followed by the estimated size θ₅₀₀. This is logical because we do not expect to find an SZ signal at the position of an AGN, and we expect them to be small (they are point-like sources). A classification considering only these two parameters gives the same performance as the one obtained with the five parameters. Figure 2 shows the detections used for this experiment in the (S/N)_sz–θ₅₀₀ plane, colour-coded according to the type of object with which they are associated. A red line indicates the best linear classification boundary determined by the trained SVM. A simple classification boundary of (S/N)_sz = 2 provides almost the same performance as the complete SVM classifier: 85% correct classifications, with 4% of the point sources classified as clusters and 11% of the clusters classified as point sources. So, for simplicity reasons, we decided to use only (S/N)_sz for the cluster/point source classification. Finally, since for our purity purposes we prefer to have fewer false clusters at the expense of a lower classification accuracy (and thus, lower completeness), we decided to modify the classification threshold to (S/N)_sz = 3, which provides 82% correct classifications, with only 2% of the point sources classified as clusters (and 16% of the clusters classified as point sources).

As mentioned before, the classification is based on the labels obtained by cross-matching some joint detections with a list of known clusters and known point sources. In particular, to construct the list of known clusters we used the MCXC (Piffaretti et al. 2011), ESZ (Planck Collaboration VIII 2011), PSZ1 (Planck Collaboration XXIX 2014), PSZ2 (Planck Collaboration XXVII 2016), SPT (Bleem et al. 2015b), and ACT (Hasselfield et al. 2013) cluster catalogues, and considered only confirmed clusters. Figure 3 shows the distribution of redshift, mass, and size of these clusters. To construct the list of known point sources, we took the ROSAT bright source catalogue (Voges et al. 1999) and we applied the selection criteria of MACS (Ebeling et al. 2001). All the objects in the resulting list were followed up for confirmation by MACS, so by eliminating the objects that match with a known cluster, we were left with a list of X-ray point sources. Since all these point sources belong to the northern hemisphere, the list of joint detections that we used for this test was obtained by running the joint detection algorithm on the northern hemisphere. Then, we then cross-matched our list of detections with the two lists to label our detections as clusters or point sources. This cross-match was done based on distance, with a matching radius of 2 arcmin. It is worth mentioning that this labelling may not be completely accurate, first because the catalogues of known point sources and clusters may not be completely correct, and second because the cross-matching done according to distance may introduce some incorrect matches. Therefore, the classification results reported above just provide a good idea of the real classification performance (with respect to the unknown ground truth). Finally, we want to emphasize that these classification results are based on the selected training dataset, so they cannot be generalized to the problem of distinguishing point sources from any kind of cluster.

Fig. 2 Estimated (S/N)sz and size θ500 of the joint detections in the northern hemisphere that match a known cluster or a known X-ray point source. The continuous red line shows the classification boundary provided by a SVM classifier trained with this dataset. The dashed red line shows the conservative cut that we adopted for discarding point sources.

Fig. 3 Distributions of the redshift (left), mass (middle), and size (right) of the clusters used for the training of the SVM classifier.

Fig. 4 Positions of the joint detections with respect to the RASS exposure map. Blue dots represent the 225 candidates corresponding to a false alarm rate of PFA =3.40 × 10−6. The sky map is colour-coded according to the logarithm of the RASS exposure time.

Fig. 5 Positions of the joint detections with respect to the Planck noise map. Red dots represent the 225 candidates corresponding to a false alarm rate of PFA =3.40 × 10−6. The sky map is colour-coded according to the logarithm of the Planck noise standard deviation map.

3.5 Output parameters and mass estimation

For each detection, the joint algorithm provides its position, the size θ₅₀₀ of the filter that gives the best joint S/N, the corresponding flux Y₅₀₀ and joint S/N, and the SZ and X-ray components of this S/N: (S/N)_sz and (S/N)_XR.

Additionally, the joint method also provides a value for the significance of each detection. This value is calculated from the simulation results described in Sect. 3.3 in the following way:

• 1.

  For each detection, we select the four simulations corresponding to the filter size that are closer to the mean Poisson level λ and mean Gaussian level σ₂₁₇ of the analysed map, and we convert the S/N maps obtained in the simulations into equivalent S/N maps corresponding to the exposure time and N_H of the real map, as in steps 1 to 3 of Sect. 3.3;

• 2.

  Then, we calculate the fraction of pixels in the simulations with a S/N on these transformed S/N maps greater than the joint S/N of the detection. This measures the probability of a false detection;

• 3.

  Finally, we perform a 2D interpolation using these four values to obtain the probability that the detection is due to noise. From this probability, we calculate the value of significance corresponding to a Gaussian distribution;

• 4.

  We note that if the joint S/N is very high, there are no pixels in the simulations with a higher S/N. In these cases, itis not possible to calculate the significance directly, and we use the following expression to estimate it: ). Appendix B explains how this expression was obtained.

Finally, since the size estimation is not very accurate, as occurred for PSZ2 catalogue, the blind detection provides the degeneracy curves Y₅₀₀(θ₅₀₀) and (S/N)_J(θ₅₀₀) for the assumed reference redshift z_ref, which allow us to determine more precisely the size and flux of the cluster given some a priori information about the cluster (e.g. redshift).

Apart from the degeneracy curve Y₅₀₀(θ₅₀₀) corresponding to the reference redshift z_ref, we can also re-extract the degeneracy curves for different redshifts at the position given by the blind detection. Then, if the detection matches a cluster with known redshift, we can interpolate between these degeneracy curves to obtain the curve corresponding to the real redshift of the cluster. This size–flux degeneracy can be further broken using the relation, which relates θ₅₀₀ and Y₅₀₀ when z is known, as explained in Sect. 7.2.2 of Planck Collaboration XXIX (2014). In this way, we can obtain an estimate of the mass M₅₀₀ of the candidate. In Sect. 3.5 we compare the mass estimated following this approach with the published mass for some of the joint detections that match known clusters.

Fig. 6 Distance from the joint position to the position of the closest cluster vs. the distance normalized to the cluster size. Only the objects with known redshift and mass in the considered SZ and X-ray catalogues were taken as clusters.

Fig. 7 M500 estimated from the joint detection for the candidates matching a confirmed cluster vs. the published M500 of the corresponding cluster. Panel a: comparison of the estimated mass to the published SZ mass for the 155 candidates matching an SZ cluster with known mass. Panel b: comparison of the estimated mass to the published MCXC mass for the 72 candidates matching an MCXC cluster. The dotted red line indicates the line of zero intercept and unity slope. The solid blue line indicates the median ratio. The dashed blue lines indicate the interval of ± 2.5σ around the median ratio. Outliers are highlighted with a blue (high estimated mass) or red (low estimated mass) square; see text. Green crosses indicate SPT clusters that are not in Planck catalogues and cyan crosses indicate MCXC clusters that are notin Planck or SPT.

Fig. 8 S/N maps for the two outliers with underestimated mass with respect to the published SZ mass and X-ray mass: PSZ2 G265.21-24.83 (left) and PSZ2 G269.36-47.20 (right). The colour images show the S/N corresponding to the X-ray map and the contours indicate the S/N corresponding to the SZ maps. The black crosses indicate the positions of the joint detection. The angular size of the images is 86.7 arcmin. The distance between the joint detection and the published SZ position is 8.4 and 7.4 arcmin, respectively.

4 Evaluation in SPT area

In this section we present an evaluation of the proposed blind detection method in the region of the sky covered by the SPT survey. This region was selected because it is a wide-area region of the sky (2500 deg²) where we can assume that almost all the massive clusters (M₅₀₀ > 7 × 10¹⁴ M_⊙) up to redshift 1.5 are already known. On the one hand, the SPT survey, which is deeper than the Planck survey, is almost 100% complete at z > 0.25 for clusters with mass M₅₀₀ > 7 × 10¹⁴ M_⊙ (Bleem et al. 2015b) (~90% complete for M₅₀₀ > 6 × 10¹⁴ M_⊙). On the other hand, the X-ray MCXC catalogue (Piffaretti et al. 2011) should include almost all the clusters with mass M₅₀₀ > 5 × 10¹⁴ M_⊙ at z < 0.25 since it contains the REFLEX sample, which is highly complete (>90%) for that redshift-mass range (see Fig. 10). There is however a small redshift range around 0.2–0.3 where some massive clusters could be still unknown. Less massive clusters could also be unknown in a broader redshift range. A comparison of the blind detection results with these catalogues allows us to determine whether the detected candidates are real clusters (purity) and what fraction of real clusters we actually detect (detection efficiency). Nevertheless, we should keep in mind that there could be some clusters in the transition region that are neither in SPT nor in MCXC. Also, other clusters could be missing due to masked regions in the surveys.

We ran the blind joint detection algorithm on the SPT footprint⁵ and obtained 2767 detections in the second phase (using q = 4). Then, we applied the cuts in the Planck S/N (S/N_sz > 3), the RASS exposure time (t_exp > 100 s) and the joint S/N (S/N_J > q_J), and we applied the SZ cleaning procedure described in Sect. 3.4.1. If we choose a false alarm rate of P_FA = 3.40 × 10⁻⁶ to calculate the joint S/N threshold q_J to apply to the 2767 detections, we are left with 225 candidates. Table C.1 summarizes the properties of these candidates. If we decrease this false alarm probability we get fewer candidates; for example, for P_FA = 2.04 × 10⁻⁷ (equivalent to a 5σ cut in a Gaussian distribution) we get 185 candidates and for P_FA = 1.90 × 10⁻⁸ (equivalent to a 5.5σ cut in a Gaussian distribution) we get 165 candidates.

The detection area is not covered homogeneously: there are slightly more candidates in the regions where the RASS exposure time is higher and where the Planck noise is lower. This is expected, since in those regions both surveys are deeper. This effect is shown in Figs. 4 and 5.

Fig. 9 Distance between the SPT and joint positions (blue) and between the SPT and PSZ2 positions (red). The median and mean values of each set of distances are represented by solid and dashed lines, respectively. Upper panel: histograms. Lower panel: corresponding cumulative distributions.

Fig. 10 Massand redshift of the clusters in the MMF3 catalogue and of the clusters detected with the proposed method for (a) PFA = 3.40 × 10−6 and (b) PFA = 2.04 × 10−7. Open circles represent the MMF3 confirmed clusters in the considered region, while filled symbols represent the joint detections colour-coded according to the associated cluster. Yellow filled circles represent joint detections matching confirmed MMF3 clusters, green filled circles represent joint detections matching other confirmed Planck clusters (not MMF3), red squares represent joint detections matching confirmed SPT clusters not detected by Planck, and blue diamonds represent joint detections matching confirmed MCXC clusters that do not match any of the previously mentioned catalogues. The blue solid line shows the REFLEX detection limit, calculated from the REFLEX flux limit and the LX − M500 relation presented in Piffaretti et al. (2011). It corresponds to a completeness of at least 90% (Böhringer et al. 2001). The green solid line shows the Planck mass limit for the SPT zone at 20% completeness.

4.1 Crossmatch with other cluster catalogues

To estimate the purity and the detection efficiency of these catalogues, we cross-matched all the candidates with various published catalogues of clusters. In particular, we took several SZ-selected catalogues covering the considered region, namely the three Planck catalogues: ESZ (Planck Collaboration VIII 2011), PSZ1 (Planck Collaboration XXIX 2014), and PSZ2 (Planck Collaboration XXVII 2016), the SPT catalogue (Bleem et al. 2015b), and the ACT catalogue (Hasselfield et al. 2013). It is worth mentioning that a subsample of the PSZ2 catalogue, namely the MMF3 sub-catalogue, is especially interesting for us since the proposed joint detection method is based on the MMF3 detection method. The SPT and ACT surveys are deeper than the Planck survey, so these catalogues contain additional clusters that were not detected by Planck. We also took as reference the X-ray selected MCXC catalogue (Piffaretti et al. 2011). This is a metacatalogue of X-ray detected clusters that was constructed from publicly available cluster catalogues of two kinds: RASS-based catalogues, obtained from the RASS survey data, and serendipitous catalogues, based on deeper pointed X-ray observations. Finally, we also considered the optically selected Abell catalogue (Abell et al. 1989). We did not use the Zwicky and redMaPPer catalogues since they do not contain clusters in the considered region.

To decide whether our detections match these previously known candidates, we first determined the closest cluster to each of our detections. To this end, we selected only the objects in the considered SZ and X-ray catalogues with known redshift and mass (i.e. confirmed clusters). Figure 6 shows a scatter plot of the absolute distance versus the relative distance (in terms of θ₅₀₀) between theassociated objects. We observe two types of associations: those with a small distance in absolute and in relative terms, and those with a long distance in absolute and in relative terms. The first group of points represents true detections of clusters, whereas the second group corresponds to the detections that are randomly distributed with respect to the considered known clusters. From this observation, we decided to use the following association rule: if the distance is <10 arcmin the detection is considered as associated with a known cluster, otherwise the detection is considered as not associated with a known cluster. We show in Sect. 4.2 that the resulting associations are valid, since the masses of the detected objects and the associated clusters agree. Furthermore, given that the considered catalogues also contain objects without redshift and mass information, we decided not to introduce an additional criterion based on the relative distance, which can be only calculated if the θ₅₀₀ of the object isknown. This association rule is very simple, but has the advantage that it can be applied to all the candidates in the considered catalogues.

After the cross-match of our candidates with these published catalogues, we found that 187 of the 225 detections corresponding to P_FA = 3.40 × 10⁻⁶ match a previously known confirmed cluster within a distance of 10 arcmin. This corresponds to a purity >83.1%. For the case of P_FA = 2.04 × 10⁻⁷ we found that 166 of the 185 detections match a previously known confirmed cluster, whereas for P_FA = 1.90 × 10⁻⁸ there are 151 matches out of 165 detections. This corresponds to a purity >89.7% and >91.5%, respectively, which is higher than before, as expected, since we decreased the false detection probability. Table 4 shows more details about the number of candidates matching the different cluster catalogues that we considered.

In this context, and for the rest of this section, we have defined “purity” as the percentage of joint detections that are associated with a confirmed cluster in these published catalogues. It is important to keep in mind that these values of purity are just rough estimations since our simple association rule could introduce a few incorrect associations. Furthermore, the candidates without a match may also be real clusters that have not been detected or included in the published catalogues (objects in masked regions, objects in a mass-redshift region where the considered surveys are not complete, etc.). Therefore, the value of purity with respect to previously known confirmed clusters can be considered an approximate lower limit.

On the other hand, we define “detection efficiency” as the percentage of candidates in the considered published catalogues that are detectedby our joint algorithm. This magnitude is related to the completeness. To calculate it, we have cross-matched all the previously known clusters in the considered region (SPT region with RASS exposure time >100 s, outside the PSZ2 masked region) withour detections. Table 4 shows the results from this cross-match for different values of P_FA. A higher P_FA allows us to recover more clusters, but at the expense of a lower purity, as seen at the bottom line of Table 4. As expected, the proposed method is able to recover most of the MMF3 objects, given that it is based on the same data and a similar approach. Amore detailed comparison with respect to MMF3 can be found in Sect. 4.4. The recovery rate of PSZ2 and Planck clusters is also high. On the contrary, the method recovers only a small fraction of SPT clusters not detected by Planck, which was foreseen, since the SPT survey data is deeper. Finally, 68.9% of the RASS clusters situated in the considered region are recovered. Given that the proposed method also uses RASS observations, this value may seem low; however, it is not due to the detection algorithm, but to the additional cut we imposed to discard possible X-ray point sources. A more detailed comparison with respect to RASS clusters can be found in Sect. 4.5.

4.2 Mass comparison

Following the procedure described in Sect. 3.5 and using the relation proposed in Planck Collaboration XX (2014), we estimated the mass M₅₀₀ for the detections matching confirmed SZ or MCXC clusters. There are a total of 163 detections matching a confirmed SZ cluster (155 with known mass) and 72 detections matching MCXC clusters. Figure 7 shows the relation between the estimated mass and the published mass for the corresponding clusters.

The comparison with respect to the published SZ mass (Fig. 7a) shows that the estimated mass closely follows the published mass; the median ratio is 0.98, very close to 1. We identified four outliers that are at more than 2.5σ from the median ratio, two with overestimated masses and two with underestimated masses. Figure 7a also shows that the SPT clusters that were not detected by Planck have, on average, a higher mass ratio than those detected by Planck. This behaviour can be explained by the Malmquist bias.

The two outliers with overestimated mass correspond to clusters PSZ2 G252.99-56.09 (also RXC J0317.9-4414, ABELL 3112) and PSZ2 G348.46-64.83 (also SPT-CLJ2313-4243, RXC J2313.9-4244, and ACO S 1101). Both clusters are known strong cool-cores according to the classification of Hudson et al. (2010), so the assumed F_X ∕Y₅₀₀ relation used in the detection does not represent these clusters accurately. Since our mass estimation is obtained from the combinationof X-ray and SZ information, its value compared to the SZ mass is boosted due to the high X-ray luminosity. We checked that the estimated mass using only the SZ information agrees with the published mass, which supports this explanation.

The twooutliers with underestimated masses correspond to clusters PSZ2 G265.21-24.83 (also RXC J0631.3-5610) and PSZ2 G269.36-47.20 (also RXC J0346.1-5702 and ABELL 3164). They can be justified by the high distance between the joint detection and the published position, which is 8.4 and 7.4 arcmin, respectively. This implies that the SZ signal at the detected position is not at its peak value, which explains the low value obtained for the mass. The reason for this distance is that the detection is centred on the X-ray peak, while the X-ray emission is not coincident with the SZ emission, as can be seen in Fig. 8. We note that in both cases the distance normalized by the cluster size is <1, so the association can be still considered correct.

The comparison with respect to the published MCXC mass (Fig. 7b) shows that the ratio between theestimated mass and the published mass is >1, with a median value of 1.19. The same value is found for the ratio between the published SZ mass and the published MCXC mass for the same clusters. This behaviour was also observed by the Planck Collaboration XXVII (2016) when they compared the SZ mass and the X-ray luminosity of common PSZ2-MCXC objects. We identified two outliers that are at more than 2.5σ from the median ratio. They coincide with the two outliers with underestimated mass with respect to published SZ mass (Fig. 7a).

This mass comparison indicates that the joint extraction provides in general a good mass proxy when the redshift is known. The main sources of bias in the mass estimation are the presence of a cool core, which tends to overestimate the mass, and an offset between the X-ray and SZ peaks, which tends to give underestimated masses.

The good agreement between the estimated and the published masses also indicates that the 10 arcmin association rule is appropriate.

4.3 Position accuracy

Since the proposed method combines Planck maps with RASS observations, which have better position accuracy due to the smaller beam, we expect the positions provided by the joint detection method to be more accurate than those provided by Planck. To assess this accuracy, we took as reference the positions given in the SPT catalogue, which are more accurate than Planck positions. Then, we selected the joint detections that match clusters detected both by SPT and PSZ2 and calculated the distance between the joint position and the SPT position. Finally, we compared this distance to the distance between the SPT and the PSZ2 position for the same clusters. Figure 9 summarizes this comparison. On average, the joint position is closer to the SPT position than the PSZ2 position is. Therefore, we can conclude that the joint detection method introduces a gain in the position determination with respect to Planck, thanks to the use of the X-ray information.

4.4 Performance comparison with the MMF3 method

Since the proposed joint detection method is built as an extension of the MMF3 detection method, we expect it to have a better performancethan that of its predecessor.

As shown in Table 4, the proposed method is able to recover most of the MMF3 candidates for P_FA = 3.40 × 10⁻⁶ and only eight MMF3 candidates are missing, three of which are confirmed clusters. The three missing clusters were initially detected (in the second phase), but then were discarded because the joint S/N was not high enough. The five MMF3 unconfirmed clusters that do not appear in our candidate list are missing for several reasons: two of them are below the initial threshold q = 4 used to include S/N peaks in the list, one has a joint S/N lower than the corresponding threshold q_J, and two were discarded because they have a (S/N)_sz < 3. These five objects belong to the PSZ2 catalogue, but they are not detected by SPT or ACT, and they do not have any known external counterpart. The recovery rate for MMF3 confirmed clusters is much higher than that of MMF3 candidates, which indicates that thanks to the combination with the X-ray information, our joint detection is able to clean the MMF3 catalogue of non-cluster objects.

Even though the proposed method misses a small fraction of the MMF3 clusters, it detects other previously known clusters that are missed by MMF3 (see Table 4). In particular, for P_FA = 3.40 × 10⁻⁶, it detects 16 additional Planck clusters, 32 SPT clusters that were not detected by Planck, and 9 MCXC clusters that were not detected by Planck or SPT. The overall effect is an improvement of the purity-detection efficiency performance with respect to the reference method MMF3. A comparison of the two methods can be seen in Fig. 10, which shows the MMF3 clusters and the joint detections in the mass-redshift plane for P_FA = 3.40 × 10⁻⁶ and P_FA = 2.04 × 10⁻⁷. The MMF3 clusters are represented as black open circles, whereas the joint detections are represented as coloured symbols. For P_FA = 3.40 × 10⁻⁶, the proposed method is able to recover almost all the MMF3 clusters while detecting at the same time additional clusters down to a mass of 2.6 × 10¹⁴ M_⊙ at redshift 0.5. For P_FA = 2.04 × 10⁻⁷, the proposed method recovers fewer clusters due to the increased purity.

Figures 11 and 12 show two examples of additional clusters detected by the proposed method thanks to the combination of SZ and X-ray information. Figure 11 shows SPT-CLJ0351-4109, a cluster at z = 0.68 with M₅₀₀ = 4.26 × 10¹⁴ M_⊙ detected by SPT, but not detected by Planck. The SZ S/N obtained from Planck observations is too low to pass the Planck detection threshold. However, the presence of some X-ray photons at the same position (11 photons within a 4 arcmin-radius circle, compared to 3.5 photons expected from the average background level) boosts the joint S/N so that the cluster is detected. The red smooth region on the right side of the middle panel of Fig. 11 is due to a negative ripple introduced by the filter around a very strong X-ray source (X-ray S/N of 75.4) situated to the right of the cluster, at a distance of 39 arcmin (outside the region represented here). Figure 12 shows RXC J0211.4-4017, a cluster at z = 0.1 with M₅₀₀ = 1.65 × 10¹⁴ M_⊙ included in the MCXC catalogue, but not detected by Planck or SPT. In this case, the presence of a strong X-ray signal at the same position of a faint SZ signal allows the detection of this cluster.

A directcomparison of the purity-detection efficiency performance of the joint detection method and MMF3 can be seen in Fig. 13. The purity and the detection efficiency are calculated with respect to all confirmed clusters from Planck, SPT, and MCXC catalogues in the considered region, thus, they are both rough estimations. Nevertheless, they serve as indicators to compare our method with the reference MMF3. The figure shows different operational points of both methods. For MMF3, the operational point is chosen through the S/N threshold. For the nominal Planck catalogue, this threshold is set to 4.5, but different thresholds can be used, producing catalogues with different purity and detection efficiency (thus, different completeness). The proposed joint method can be tuned by changing the false alarm probability that is used to calculate the joint S/N threshold. This figure shows that our detection method outperforms MMF3 in the sense that it can simultaneously achieve higher purity and higher detection efficiency if the operational point is chosen appropriately.

Fig. 11 S/N maps for the joint detection that matches SPT-CLJ0351-4109, an SPT cluster not detected by Planck situated at z = 0.68 with M500 = 4.26 × 1014 M⊙. The three colour images show the S/N corresponding to the SZ filtered maps (left), the X-ray filtered map (middle), and the joint filtered maps (right). The filter size is 0.8 arcmin, which provides the best S/N for this detection. The angular size of the images is 68.7 arcmin. The RASS exposure time at the position of the detection is 545 s.

Fig. 12 S/N maps for the joint detection that matches RXC J0211.4-4017, a MCXC cluster not detected by Planck or SPT situated at z = 0.1 with M500 = 1.65 × 1014 M⊙. The three colour images show the S/N corresponding to the SZ filtered maps (left), the X-ray filtered map (middle) and the joint filtered maps (right). The filter size is 0.8 arcmin, which corresponds to the one that provides the best S/N for this detection. The angular size of the images is 68.7 arcmin. The RASS exposure time at the position of the detection is 946 s.

Fig. 13 Comparison of the performance of the proposed method (green, blue, cyan, and magenta circles) and the MMF3 method (red squares). The horizontal axis shows the percentage of detections which match a real cluster within a 10 arcmin radius. It is an indicator of the purity of the methods. The vertical axis shows the number of real clusters detected (which match a detection within a 10 arcmin radius). It is related to the detection efficiency of the methods. We assume thatthe real clusters in the considered region are all the clusters from the SPT, MCXC, and Planck catalogues with known mass and redshift (confirmed). The different points in the curves correspond to different operational points of the detection algorithms. For the MMF3 case we have represented the results for S/N thresholds of 4.00, 4.25, 4.50, 4.75, and 5.00, increasing from left to right. For the proposed algorithm, we have represented the results for PFA =3.40 × 10−6, PFA =2.04 × 10−7 and PFA = 1.90 × 10−8, decreasing from left to right for each (S/N)sz cut.

4.5 Comparison with RASS clusters

Since the proposed joint detection method uses RASS observations, it is interesting to check whether it is able to recover known clusters that have been detected using the same observations. Table 4 shows that we detect 71 of the 103 RASS clusters situated in the considered region (SPT area with RASS exposure time >100 s, outside the PSZ2 masked region), which corresponds to 68.9%. Most of the RASS clusters that we do not recover (30 of the 32) were in fact included in the list of detections provided by the second phase of the algorithm, but 29 were discarded because their (S/N)_sz was lower than 3 and one was discarded because the RASS exposure time at the detection’s position is lower than 100 s. There are just two RASS clusters that were not originally detected by the joint algorithm because their joint S/N does not reach the threshold of q = 4: RXC J0040.1-5607 and RXC J2326.7-5242, which are quite faint both in X-ray and SZ. To cross-check these results, we used the MMF3 method of Planck Collaboration XXVII (2016) to extract the S/N of RASS clusters from Planck maps. We found that the 2 undetected clusters and the 29 clusters that were discarded due to a low (S/N)_sz have a very low S/N, which supports our results. Therefore, we can conclude that the joint detection method is able to recover almost all the RASS clusters, as expected, but we discard some of them later in order to maintain a high purity by eliminating possible AGN detections with a threshold in (S/N)_sz, which has a similar effect to a mass cut at each redshift. Figure 14 illustrates this comparison by showing the RASS clusters and the joint detections in the mass-redshift plane.

Fig. 14 Mass and redshift of the clusters in the RASS catalogue. Open circles represent RASS clusters in the considered region; red filled circles represent RASS clusters that are detected by the joint detection algorithm; cyan filled circles represent RASS clusters that were detected in the second phase of the algorithm, but discarded due to a low (S/N)sz; and the green-filled circle represents the RASS cluster that was detected in the second phase of the algorithm, but discarded due to a low exposure time.

4.6 New candidates

As mentioned before, 193 of the 225 detections corresponding to P_FA = 3.40 × 10⁻⁶ match with a real cluster within a distance of 10 arcmin. This means that our catalogue contains 32 candidates that are not known clusters. They could be either false detections due to noise or to X-ray point sources, or true clusters not detected before. Table 5 summarizes the coordinates and some additional information of these 32 new candidates.

Since weset a false alarm probability of P_FA = 3.40 × 10⁻⁶ to calculate the joint S/N threshold, we expect to have 1.7 pixels in each filtered patch with (S/N)_J greater than the threshold. This means that in the whole SPT area we expect to have 42 pixels above the threshold, producing at most 42 false detections due to noise fluctuations. This number is close to the number of new candidates, so can we expect most of them to be false detections. On the other hand, for P_FA = 2.04 × 10⁻⁷ and P_FA = 1.90 × 10⁻⁸ we would expect at most 2.5 and 0.2 false detections due to noise fluctuations in the SPT footprint, respectively. However, the number of new candidates for these two false alarm probabilities is 14 and 10, respectively (see Table 4). We expect then that most of these candidates are real detections (either clusters or other objects like X-ray point sources). We indicate in Table 5 which of the 32 candidates corresponding to P_FA = 3.40 × 10⁻⁶ are also candidates for P_FA = 2.04 × 10⁻⁷ and P_FA = 1.90 × 10⁻⁸.

We searched for archival X-ray observations covering these 32 positions. Three of the candidates (1, 8, and 12) were observed by Swift (Burrows et al. 2005) and two more (14 and 32) by XMM-Newton. In the three Swift observations, there is evidence that the candidates are real clusters, since we can see an extended emission, as shown in Fig. 15. On the contrary, the two XMM-Newton observations show that the candidates are false detections: candidate 14 is a point source, while candidate 32 is just a noise fluctuation. As shown in the Table 5, the joint S/N of candidate 32 is just above the threshold.

We also looked around these 32 positions for other known galaxy clusters or groups in the NED⁶ and SIMBAD⁷ databases. Table 6 shows all the clusters and groups found closer than 10 arcmin to our candidates. Most of the joint candidates do not have any NED or SIMBAD clusters close to them. For 13 candidates, we found some close-by objects, but in most of the cases, they do not seem to be associated with the candidate since their separation is too big. Only candidates 12 and 21 might be associated with clusters: the SPT cluster SPT-CL J0438-4907 and the optical cluster LCS-CL J051723-5325.5, respectively.

We clarify that we did not match candidate 12 with SPT-CL J0438-4907 before (Sect. 4.1) because this cluster is not included in the published SPT catalogue of Bleem et al. (2015b) since its significance is lower than 4.5, which was the limit for the published catalogue. However, it is detected at lower significance and confirmed with opticalobservations in Saro et al. (2015). The presence of this cluster at only 0.7 arcmin from candidate 12 and the Swift observation presented in Fig. 15c are strong indicators that this candidate is a real cluster. Moreover, the mass obtained from the joint extraction assuming the redshift of the SPT cluster is M₅₀₀ = 3.04 × 10¹⁴ M_⊙, very close to the mass published by Saro et al. (2015) for the SPT cluster (M₅₀₀ = (3.13 ± 0.81) × 10¹⁴ M_⊙).

Regarding candidate 21, the estimated mass assuming the redshift of the nearby optical cluster is M₅₀₀ = 3.38 × 10¹⁴ M_⊙. If we apply the richness-mass relation of Rozo et al. (2015), we get an estimated richness of λ_e = 63.5. According to Bleem et al. (2015a), the optical cluster LCS-CL J051723-5325.5 has a richness of λ = 29.2, which differs from λ_e by 2.92 σ_lnλ. Given the large scatter of the richness-mass relation, it is reasonable to associate candidate 21 with cluster LCS-CL J051723-5325.5.

Fig. 15 Swift observations of candidates 1 (a), 8 (b), and 12 (c). A red circle with 2 arcmin radius centred at the candidate position wasadded for visual reference.

5 Conclusions

In this paper we have proposed a galaxy cluster detection method based on matched multifrequency filters (MMF) that combines X-ray and SZ observations. This method builds on the previously proposed joint X-ray–SZ extraction method and allows us to blindly detect clusters, i.e. finding new clusters without knowing their position, size, or redshift, by searching on SZ and X-ray maps simultaneously. It can be seen as an evolution of the MMF3 detection method, one of the MMF methods used to detect clusters from Planck observations, that incorporates X-ray observations to improve the detection performance.

The main challenge to solve was to obtain a high purity, since the addition of the X-ray information increases the cluster detection probability, but also the number of false detections, produced by AGNs and Poisson noise. To deal with the Poisson noise correctly, we proposed an adaptive S/N threshold to keep or discard detections depending on the noise characteristics of the region. To discard AGN detections, we proposed an additional classification according to the SZ part of the S/N.

The proposed method is tested using data from the RASS and Planck surveys and evaluated by comparing the detection results with existing cluster catalogues in the area of the sky covered by the SPT survey. We have shown that, thanks to the addition of the X-ray information, the method is able to simultaneously achieve better purity, better detection efficiency, and better position accuracy than its predecessor, the MMF3 detection method.

We have also shown that if the redshift of a cluster is known by any other means, the joint detection allows a good estimation of its mass. Some bias may appear in the presence of a cool core (overestimated mass) or if there is an offset between the X-ray and SZ peaks (underestimated mass).

Finally, we have produced a catalogue of candidates in the SPT region composed of 225 objects, with 32 new objects that are not included in other SZ or X-ray cluster catalogues. We have found, using Swift observations, that three of these new objects are probably real clusters. This supports the fact that the proposed method can be used to find new clusters.

In future work we will run the joint detection on all the sky using Planck and RASS maps and we will provide the last and deepest all-sky cluster catalogue before the e-ROSITA mission.

Acknowledgements

This research is based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. This research has made use of the ROSAT all-sky survey data which have been processed at MPE. The authors acknowledge the use of the HEALPix package (Górski et al. 2005). The authors would like to thank Iacopo Bartalucci for his help in the preparation of the Swift images. The authors thank an anonymous referee for the constructive comments. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement n° 340519.

Appendix A: Matching of simulation results with real maps parameters

As explained in Sect. 3.3, the X-ray noise maps used in the Monte Carlo simulations were simulated as homogeneous Poisson random fields, characterized by a given mean value λ (in counts). To express these count maps in ΔT∕T units, as was done with real RASS count maps, we need to assume an exposure time map and a N_H map, and then apply the conversion procedure described in Appendix B of Tarrío et al. (2016), which can be summarized as follows: (A.1)

In this expression c(N_H) represents the factor that converts the count rate into X-ray flux and it depends on the N_H map; the expected F_X ∕Y₅₀₀ relation is used to convert the X-ray flux into equivalent Y₅₀₀ integrated flux and depends on the reference redshift; is the HEALPix pixel area and g(ν_ref) is the factor that converts from y units into ΔT∕T_CMB units, which depends on the reference frequency assumed for the map (1000 GHz in our case).

For the Monte Carlo simulations of this paper, we assumed a constant exposure time t_exp = 400 s and a constant N_H = 2 × 10²⁰ cm⁻² (average values in the SPT region). If other values were used, the resulting X-ray maps in Δ T∕T units would only differ from the ones obtained with these values by a constant factor a. This allows us to convert some of the simulation results obtained for the reference t_exp and N_H values (in particular and S/N) into the results corresponding to any other value of t_exp and/or N_H. In the following, we explain how this conversion is done.

Let M_X and be two X-ray maps in Δ T∕T units, calculated from the same count map using different values of t_exp and N_H. From Eq. (A.1), we have the following relation between the two maps: (A.2)

where (A.3)

Let and be two multifrequency noise maps whose SZ components are equal and whose X-ray components differ by a constant factor a.

Considering that the noise in the X-ray map and the SZ maps is uncorrelated, we can write the noise power spectrum of M as (A.4)

where P_sz is the noise power spectrum of the SZ maps N_sz, P_X is the noise power spectrum of the X-ray map N_X, and 0_n×m denotes a vector with n rows and m columns whose elements are all equal to 0. The noise power spectrum of M′ can be decomposed into P_sz and in an analogous way. Using the definition of the noise power spectrum (see Sect. 3.1), it is immediate to show that .

Using the definition of the variance of the filtered maps (Eq. (2)) and applying (A.4), we can also decompose the variance of the filtered maps into an SZ and an X-ray component as (A.5)

where F_sz and F_X are the SZ and X-ray components of (Eq. (3)). From this expression it is easy to show that (A.6)

If we filter M and M′ with the proposed joint filter (Eq. (1)), the S/N of the filtered maps will be given by (A.7)

and (A.8)

From these two expressions, and taking into account (A.5) and (A.6), we can obtain the relation between the S/N of the two noise maps as (A.9)

where (A.10)

and (A.11)

Appendix B: Significance estimation

As explained in Sect. 3.5, the significance of a detection is calculated with the aid of the S/N maps obtained from the numerical simulations described in Sect. 3.3. In particular, it is obtained by measuring the number of pixels in the simulations that have a S/N higher than the S/N of the detection. Given that the number of simulated pixels is finite, when the joint S/N of the detection is very high, there may be no pixels satisfying this condition. In these cases, it is not possible to calculate the significance directly. Since for each combination of mean Poisson level λ and mean Gaussian level σ₂₁₇ we created 450 random realizations of the noise maps, and since each map has 700 × 700 pixels, the minimum false alarm probability that we can determine is 4.5 × 10⁻⁹. This means that when we do not find any pixels above a given S/N, we can only affirm that the significance will be higher than 5.75.

To find a way to estimate the significance in these cases, we analysed the second phase detections corresponding to P_FA = 3.40 × 10⁻⁶. We calculated the significance of the detections with (S/N)_sz > 3 outside the SZ mask and we found a good correlation between the significance and the difference between the joint S/N and the threshold q_J, as shown in Fig. B.1. Therefore, we decided to use a linear extrapolation to estimate the values of significance as a function of (S/N)_J − q_J for the detections whose calculated significance is >5.75. By definition, the significance corresponding to (S/N)_J = q_J is 4.5 for P_FA = 3.40 × 10⁻⁶. Thus, the linear extrapolation was found by fitting a line with a fixed intercept of 4.5 to the points in Fig. B.1, excluding the outliers. The final expression we obtained is (B.1)

Fig. B.1 Significance as a function of (S/N)J − qJ for the second phase detections corresponding to PFA = 3.40 × 10−6, with (S/N)sz > 3 and outside the SZ mask. Red circles indicate the detections that pass all the other cuts (texp > 100 s and (S/N)J > qJ), while black circles indicate the detections that do not pass them. The solid blue line is the best linear fit with a fixed intercept of 4.5.

Appendix C: List of candidates

References

All Tables

All Figures

	Fig. 1 Examples of the matched filter for θs = 1 arcmin (a) and θs = 30 arcmin (b). The curves give the radial profiles of the filters, which are symmetric because we have chosen a symmetric cluster template. The filter is normalized so that its maximum amplitude is equal to 1.
In the text

	Fig. 2 Estimated (S/N)sz and size θ500 of the joint detections in the northern hemisphere that match a known cluster or a known X-ray point source. The continuous red line shows the classification boundary provided by a SVM classifier trained with this dataset. The dashed red line shows the conservative cut that we adopted for discarding point sources.
In the text

	Fig. 3 Distributions of the redshift (left), mass (middle), and size (right) of the clusters used for the training of the SVM classifier.
In the text

	Fig. 4 Positions of the joint detections with respect to the RASS exposure map. Blue dots represent the 225 candidates corresponding to a false alarm rate of PFA =3.40 × 10−6. The sky map is colour-coded according to the logarithm of the RASS exposure time.
In the text

	Fig. 5 Positions of the joint detections with respect to the Planck noise map. Red dots represent the 225 candidates corresponding to a false alarm rate of PFA =3.40 × 10−6. The sky map is colour-coded according to the logarithm of the Planck noise standard deviation map.
In the text

	Fig. 6 Distance from the joint position to the position of the closest cluster vs. the distance normalized to the cluster size. Only the objects with known redshift and mass in the considered SZ and X-ray catalogues were taken as clusters.
In the text

Fig. 7 M500 estimated from the joint detection for the candidates matching a confirmed cluster vs. the published M500 of the corresponding cluster. Panel a: comparison of the estimated mass to the published SZ mass for the 155 candidates matching an SZ cluster with known mass. Panel b: comparison of the estimated mass to the published MCXC mass for the 72 candidates matching an MCXC cluster. The dotted red line indicates the line of zero intercept and unity slope. The solid blue line indicates the median ratio. The dashed blue lines indicate the interval of ± 2.5σ around the median ratio. Outliers are highlighted with a blue (high estimated mass) or red (low estimated mass) square; see text. Green crosses indicate SPT clusters that are not in Planck catalogues and cyan crosses indicate MCXC clusters that are notin Planck or SPT.

In the text

Fig. 8 S/N maps for the two outliers with underestimated mass with respect to the published SZ mass and X-ray mass: PSZ2 G265.21-24.83 (left) and PSZ2 G269.36-47.20 (right). The colour images show the S/N corresponding to the X-ray map and the contours indicate the S/N corresponding to the SZ maps. The black crosses indicate the positions of the joint detection. The angular size of the images is 86.7 arcmin. The distance between the joint detection and the published SZ position is 8.4 and 7.4 arcmin, respectively.

In the text

	Fig. 9 Distance between the SPT and joint positions (blue) and between the SPT and PSZ2 positions (red). The median and mean values of each set of distances are represented by solid and dashed lines, respectively. Upper panel: histograms. Lower panel: corresponding cumulative distributions.
In the text

Fig. 10 Massand redshift of the clusters in the MMF3 catalogue and of the clusters detected with the proposed method for (a) PFA = 3.40 × 10−6 and (b) PFA = 2.04 × 10−7. Open circles represent the MMF3 confirmed clusters in the considered region, while filled symbols represent the joint detections colour-coded according to the associated cluster. Yellow filled circles represent joint detections matching confirmed MMF3 clusters, green filled circles represent joint detections matching other confirmed Planck clusters (not MMF3), red squares represent joint detections matching confirmed SPT clusters not detected by Planck, and blue diamonds represent joint detections matching confirmed MCXC clusters that do not match any of the previously mentioned catalogues. The blue solid line shows the REFLEX detection limit, calculated from the REFLEX flux limit and the LX − M500 relation presented in Piffaretti et al. (2011). It corresponds to a completeness of at least 90% (Böhringer et al. 2001). The green solid line shows the Planck mass limit for the SPT zone at 20% completeness.

In the text

Fig. 11 S/N maps for the joint detection that matches SPT-CLJ0351-4109, an SPT cluster not detected by Planck situated at z = 0.68 with M500 = 4.26 × 1014 M⊙. The three colour images show the S/N corresponding to the SZ filtered maps (left), the X-ray filtered map (middle), and the joint filtered maps (right). The filter size is 0.8 arcmin, which provides the best S/N for this detection. The angular size of the images is 68.7 arcmin. The RASS exposure time at the position of the detection is 545 s.

In the text

Fig. 12 S/N maps for the joint detection that matches RXC J0211.4-4017, a MCXC cluster not detected by Planck or SPT situated at z = 0.1 with M500 = 1.65 × 1014 M⊙. The three colour images show the S/N corresponding to the SZ filtered maps (left), the X-ray filtered map (middle) and the joint filtered maps (right). The filter size is 0.8 arcmin, which corresponds to the one that provides the best S/N for this detection. The angular size of the images is 68.7 arcmin. The RASS exposure time at the position of the detection is 946 s.

In the text

Fig. 13 Comparison of the performance of the proposed method (green, blue, cyan, and magenta circles) and the MMF3 method (red squares). The horizontal axis shows the percentage of detections which match a real cluster within a 10 arcmin radius. It is an indicator of the purity of the methods. The vertical axis shows the number of real clusters detected (which match a detection within a 10 arcmin radius). It is related to the detection efficiency of the methods. We assume thatthe real clusters in the considered region are all the clusters from the SPT, MCXC, and Planck catalogues with known mass and redshift (confirmed). The different points in the curves correspond to different operational points of the detection algorithms. For the MMF3 case we have represented the results for S/N thresholds of 4.00, 4.25, 4.50, 4.75, and 5.00, increasing from left to right. For the proposed algorithm, we have represented the results for PFA =3.40 × 10−6, PFA =2.04 × 10−7 and PFA = 1.90 × 10−8, decreasing from left to right for each (S/N)sz cut.

In the text

Fig. 14 Mass and redshift of the clusters in the RASS catalogue. Open circles represent RASS clusters in the considered region; red filled circles represent RASS clusters that are detected by the joint detection algorithm; cyan filled circles represent RASS clusters that were detected in the second phase of the algorithm, but discarded due to a low (S/N)sz; and the green-filled circle represents the RASS cluster that was detected in the second phase of the algorithm, but discarded due to a low exposure time.

In the text

	Fig. 15 Swift observations of candidates 1 (a), 8 (b), and 12 (c). A red circle with 2 arcmin radius centred at the candidate position wasadded for visual reference.
In the text

	Fig. B.1 Significance as a function of (S/N)J − qJ for the second phase detections corresponding to PFA = 3.40 × 10−6, with (S/N)sz > 3 and outside the SZ mask. Red circles indicate the detections that pass all the other cuts (texp > 100 s and (S/N)J > qJ), while black circles indicate the detections that do not pass them. The solid blue line is the best linear fit with a fixed intercept of 4.5.
In the text