© ESO 2018

1. Introduction

Quantifying dust emission from galaxy clusters is interesting for both astrophysical and cosmological studies. Dust emission from member galaxies is a tracer of the star formation rate (SFR) in dense environments (Alberts et al. 2014, 2016), and the question of intracluster dust embedded in the hot intracluster medium (ICM) concerns stellar feedback and the physical state of the ICM (Montier & Giard 2004). Cluster dust emission also has potentially important ramifications for cosmology from Sunyaev–Zeldovich (SZ) cluster counts because it can contaminate the SZ signal and modify survey selection functions.

One of the first appearances of contaminating dust emission was found by Planck Collaboration Int. XXXVII (2013) when stacking the SZ signal from central halo galaxies. Contamination by dust emission came to dominate the SZ signal when approaching the low-mass group scale. Planck Collaboration XXIII (2016), Planck Collaboration XLIII (2016) examine dust emission by stacking signal in the high frequency Planck maps around massive clusters from the Planck catalogs. The former work separated the dust and SZ signals to conclude that the dust emission evolved with redshift and was more spatially extended than the SZ signal. The authors of the latter work combined IRAS data with Planck observations to measure dust temperature and determine dust masses in cluster systems. These studies extend the work of Montier & Giard (2005) and Giard et al. (2008), who detected dust emission by stacking IRAS maps around clusters.

Because the Planck beam has an angular extent similar to or larger than that of the studied clusters, these observations integrate their total emission. Planck Collaboration XXIII (2016) find that the dust emission could be fully accounted for by cluster member galaxies, in agreement with previous work (Roncarelli et al. 2010), a conclusion further supported by the temperatures of T ∼ 20 K determined by Planck Collaboration XLIII (2016) that are typical of late-type galaxies.

Little attention has yet been given to studying the impact of dust emission on SZ cluster surveys, largely because the level of the emission relative to the SZ signal is poorly known. The large surveys by the Atacama Cosmology Telescope (ACT, Hasselfield et al. 2013), the South Pole Telescope (SPT, Bleem et al. 2015) and the Planck mission (Planck Collaboration XXVII 2016) do not model the effect of dust emission on their selection functions and photometry. Any effect will depend in detail on the observation bands and how they are used in cluster detection.

In this paper, we examine dust emission from massive, intermediate redshift clusters using Planck observations and evaluate its impact on the Planck SZ cluster selection function and photometry. The data are presented in Sect. 2. We proceed by first establishing a baseline model (Sect. 3) for cluster dust emission and fit key model parameters with our Planck measurements (Sect. 4). These parameters are the normalization of the infrared (IR) luminosity-cluster mass relation and the spatial extent of the dust emission. We then simulate Planck observations of clusters with both SZ signal and dust emission to quantify the effect of the dust emission on the Planck SZ selection function and photometry (size and SZ flux) in Sect. 5. We conclude in Sect. 6. Throughout the paper, we adopt the Planck ΛCDM cosmology (TT,TE,EE+lowP+lensing+ext in Table 4 of Planck Collaboration XIII 2016): h = H₀/(100 km s⁻¹ Mpc⁻¹) = 0.6774, Ω_m = 1 − Ω_Λ = 0.3089, Ω_b = 0.0485976, n_s = 0.9667.

2. Data

We used the all-sky maps of the Planck High Frequency Instrument (HFI). From August 2009 to January 2012, HFI observed the sky in six frequency bands centered on 100, 143, 217, 353, 545 and 857 GHz. We used the full mission temperature maps that can be downloaded online¹. In our analysis, we have assumed that the beam for each map is Gaussian with full width at half maximum (FWHM) values of 9.659, 7.220, 4.900, 4.916, 4.675, 4.216 arcmin, respectively, for each band.

We also used the second Planck catalog of SZ sources (PSZ2, Planck Collaboration XXVII 2016), keeping only sources with an assigned redshift (1093 objects), and re-extract their SZ signal using Multifrequency Matched Filters (MMF3, hereafter noted MMF for simplicity; Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016; Melin et al. 2012). The re-extraction gives a positive signal-to-noise ratio (S/N) for 1091 objects, which constitutes the sample that we have adopted throughout this paper. Figure 1 (top) shows the sample redshift distribution.

The MMF re-extraction provides size-flux degeneracy curves for each source that we break using an independent X-ray size-flux relation (combination of Eqs. (7) and (9) of Planck Collaboration XX 2014). This allows us to compute the mass proxy of each object, $M_{500}^{\mathrm{Yz}}$ (and the associated cluster size ${\theta }_{500}^{\mathrm{Yz}}$, since the redshift is known). More detail on the method of computing the mass proxy is given in Sect. 7.2.2 of Planck Collaboration XXIX (2014). We thus have 1091 objects with position, redshift z, mass $M_{500}^{\mathrm{Yz}}$ and size ${\theta }_{500}^{\mathrm{Yz}}$.

3. Dust modeling

Our baseline dust model is built on Herschel observations of field and cluster galaxies (Alberts et al. 2014, 2016), and on the model from Cai et al. (2013) for the luminosity functions and spectral energy distributions. The model was developed for the prediction of cluster fluxes for the CORE space mission (Delabrouille et al. 2018; De Zotti et al. 2018). We briefly summarize its main elements.

The total comoving infrared field luminosity density, Ψ_IR(z), is computed using the model from Cai et al. (2013). The model includes three populations of galaxies: the “warm” and “cold” populations dominate at z < 1, the “spheroidal” population at z > 1.5. Using the luminosity functions, Φ_i, given by the model², we compute, for each redshift, the comoving infrared luminosity density, Ψ_i(z), contributed by each galaxy population i:$$\begin{array}{c}\hfill {\Psi }_i(z)=\int {\Phi }_i(logL,z)L\text{d}logL,\end{array}$$(1)

for i = cold, warm, spheroidal. log is the base-10 logarithm. We then compute the total comoving infrared field luminosity density$$\begin{array}{c}\hfill {\Psi }_{\mathrm{IR}}(z)=\sum _i{\Psi }_i(z).\end{array}$$(2)

For z < 1.2, Alberts et al. (2014) measured the ratio of the mean infrared luminosity in clusters to that of galaxies in the field (see their Table 2)$$\begin{array}{c}\hfill f(z)=\frac{L_{\mathrm{cluster}}}{L_{\mathrm{field}}}=\frac{1540}{267}{\text{e}}^{-(0.76-0.42)t_{\mathrm{Gyr}}(z)}=5.77{\text{e}}^{-0.34t_{\mathrm{Gyr}}(z)}\end{array}$$(3)

where t_Gyr(z) is the cosmic time in Gyr.

From Ψ_IR(z) and f(z), we infer that the total infrared luminosity in the sphere of radius R_Δ for a cluster located at z < 1.2 is approximately$$L_{\Delta ,\text{tot},z\hspace{0.17em}<\hspace{0.17em}1.2}\approx {\Psi }_{\text{IR}}(z)\times f(z)\times \frac{{\rho }_{\text{cluster}}}{{\rho }_{\text{mean}}}\times V_{\text{comoving}},\hspace{0.17em}{L}_{\Delta ,\text{tot},z\hspace{0.17em}<\hspace{0.17em}1.2}\approx \frac{{10}^{14}{M}_{\odot }}{{\rho }_{\text{mean}}\hspace{0.17em}(z=0)}\times {\Psi }_{\text{IR}}(z)\times f(z)\times \frac{M_{\Delta }}{{10}^{14}{M}_{\odot }}$$(4)

where $V_{\mathrm{comoving}}={(1+z)}^3\times \frac{4\pi }{3}{R}_{\mathrm{\Delta }}^3$ is the comoving volume enclosed in the sphere of radius R_Δ, ρ_cluster the cluster density, ρ_mean the matter density of the Universe at redshift z, Δ the overdensity with respect to critical density of the Universe at redshift z, and M_Δ the mass enclosed in the sphere of radius R_Δ.

Fixing Δ = 500, we write$$\begin{array}{c}\hfill L_{500,\mathrm{tot},z<1.2}=r_L\times \frac{{10}^{14}{M}_{\odot }}{{\rho }_{\mathrm{mean}}(z=0)}\times {\Psi }_{\mathrm{IR}}(z)\times f(z)\times \frac{M_{500}}{{10}^{14}{M}_{\odot }},\end{array}$$(5)

with r_L a normalization factor that we will determine from the 857 GHz flux of Planck clusters (see Sect. 4). It is expected to be close to unity if Herschel and Planck data are consistent and our model is valid.

For z > 1.2, the factor f(z) does not apply because star formation in clusters matches that of field galaxies (Alberts et al. 2016); thus the IR luminosity reads$$\begin{array}{c}\hfill L_{500,\mathrm{tot},z>1.2}=r_L\times \frac{{10}^{14}{M}_{\odot }}{{\rho }_{\mathrm{mean}}(z=0)}\times {\Psi }_{\mathrm{IR}}(z)\times \frac{M_{500}}{{10}^{14}{M}_{\odot }}·\end{array}$$(6)

We then compute, for each redshift, the fraction p_i(z) of the total luminosity contributed by the galaxy population i as$$\begin{array}{c}\hfill p_i(z)=\frac{{\Psi }_i(z)}{{\Psi }_{\mathrm{IR}}(z)}·\end{array}$$(7)

The luminosity from each population is then given by$$\begin{array}{c}\hfill L_{500,i}[\nu (1+z)]=p_i(z)\times L_{500,\mathrm{tot}}\times {\mathrm{SED}}_i[\nu (1+z)],\end{array}$$(8)

where ν is the observation frequency and SED_i is the spectral energy distribution of the population i normalized such that ∫ SED_i[ν(1 + z)]dν = 1 (see e.g., Fig. 4 of Cai et al. 2013). The dust flux density for population i is thus$$\begin{array}{c}\hfill S_{500,i}(\nu )=\frac{(1+z)L_{500,i}[\nu (1+z)]}{4\pi D_L^2(z)},\end{array}$$(9)

with D_L(z) the luminosity distance. This model is identical to the model adopted in Sect. 4 of De Zotti et al. (2018), except that we have introduced the normalization factor r_L. Fixing r_L = 1 makes our model strictly identical to De Zotti et al. (2018).

Figure 1 (bottom) shows the predicted dust flux density integrated over the Planck 857 GHz bandwidth, S₅₀₀, as a function of redshift for a cluster of mass M₅₀₀ = 10^14.5 M_⊙ (thick solid line). It is essentially only composed of the warm (dashed blue line) and cold (dotted red line) components for redshifts z < 1. At z > 0.25, the flux density increases with z because of the increase in luminosity with z. At z < 0.25, the luminosity distance dominates the redshift evolution.

Throughout this paper, we will use this model only in the redshift range 0 < z < 1, which is relevant to the Planck cluster catalog. But the model also includes the spheroidal component, which dominates at z > 1.5, and we intend to use it in future work to examine the impact of dust emission on next generation cosmic microwave background (CMB) experiments.

Our model gives global quantities (infrared luminosity), but it does not give any information on the spatial distribution of the dust emission in clusters. In Sect. 4, we use Planck PSZ2 clusters to jointly constrain the model normalization, r_L (see Eq. (5)), and the emission profile.

4. Planck constraints on the normalization and spatial profile of cluster dust emission

We describe the three dimensional dust emission profile with a Generalized Navarro-Frenk-White (GNFW) profile (Nagai et al. 2007):$$\begin{array}{c}\hfill P(r)\propto \frac{1}{{(r/r_{\mathrm{s}})}^{\gamma }{[1+{(r/r_{\mathrm{s}})}^{\alpha }]}^{\frac{\beta -\gamma }{\alpha }}},\end{array}$$(10)

where α, β and γ are, respectively, the intermediate, external and central slopes, r_s = R₅₀₀/c₅₀₀ is the scale radius and c₅₀₀ is the concentration parameter. Dark matter profiles of massive relaxed clusters follow a standard NFW (Navarro et al. 1996) profile with α = γ = 1 and β = 3. For nearby relaxed galaxy clusters, Pointecouteau et al. (2005) have constrained c₂₀₀ = 4.61 ± 0.12, which we convert to c₅₀₀ = 3.03 ± 0.08.

We used two observables, the stacked profile (Sect. 4.1) and the inverse-variance weighted average matched filter flux (Sect. 4.2), to constrain the normalization, r_L, in Eq. (5) and the dust emission profile. The stacked profile does not provide enough information to simultaneously constrain all of the GNFW parameters. We therefore fixed α = γ = 1 and c₅₀₀ = 3, and leave only the external slope, β, free. A value of β larger (or smaller) than three indicates that the profile is steeper (or shallower) than the dark matter profile of massive relaxed clusters. The constraints from the stacked profile and the inverse-variance weighted matched filter can be first compared and then combined (Sect. 4.3).

4.1. Stacked profile

Our first observable is the stacked profile of the PSZ2 clusters (see Sect. 2), which we constructed following the same methodology as Planck Collaboration XXIII (2016). For each cluster, we computed the unweighted mean flux of the 857 GHz map in annuli of width Δθ = 1 arcmin, starting from θ = 0 (the SZ center) to θ = 30 arcmin. We removed the offset using the mean value of the pixels between θ = 30 arcmin and 60 arcmin. We then took the mean of all the profiles in each annulus.

This stacked profile is shown as the black diamonds in the lefthand panel of Fig. 2. The error bars were obtained as the standard deviation of 10 000 bootstrap realizations. Stacking the profile in angular radii θ mixes different physical scales. This procedure thus introduces correlations between the bins, which can be estimated from the bootstraps. The righthand panel of Fig. 2 shows the correlation matrix, and we see that the data points are indeed strongly correlated (>75%).

We fit the observed profile using a stacked GNFW profile. For each cluster, a GNFW profile is scaled to ${\theta }_{500}^{\mathrm{Yz}}$ and normalized using Eq. (9). We then applied the same averaging procedure as for the data. We fixed all the parameters and let only r_L and β vary. The 68%/95% confidence limits (C.L.) on r_L and β are shown as the solid and dashed blue lines respectively in Fig. 3. The best fit values are given in the first row of Table 1. The normalization, $r_L=1.27\begin{array}{c}+0.34\\ -0.33\end{array}$, is compatible with one, the value determined from Herschel data. The slope parameter, $\beta =1.36\begin{array}{c}+0.39\\ -0.30\end{array}$, is significantly lower than three, indicating that the dust profile is shallower than the matter profile. The fit of the stacked profile is shown as the blue dash-dotted line in the lefthand panel of Fig. 2. We note that it is systematically below the majority of the data points. This is due to the strong correlation between the points, as given by the correlation matrix in the righthand panel of the figure.

4.2. Inverse-variance weighted average matched filter flux

Our second observable is the inverse-variance weighted average matched filter flux of the PSZ2 clusters measured in 857 GHz maps. We extract individual cluster flux and associated error in the 857 GHz map using a single frequency matched filter (Melin et al. 2006). We fix the position to the SZ center and the size to ${\theta }_{500}^{\mathrm{Yz}}$, and we adopt the universal pressure profile from Arnaud et al. (2010). We also perform the flux extraction at the five other Planck HFI frequencies, although we do not use them to constrain r_L and β. Results are shown in Fig. 4 as black diamonds. The error bars are obtained from the standard deviation of 10 000 bootstrap realizations of the inverse-weighted average.

We then fit the GNFW model to the 857 GHz data point by applying the matched filter to the model and averaging as done on the real data. The constraints are shown as green contours in Fig. 3. There is no absolute minimum, since for each value of β we can find r_L which adjusts the average matched filter flux at 857 GHz. The contour is thus a valley.

4.3. Combination

The r_L values preferred by the inverse-variance weighted average matched filter flux are lower than the values preferred by the stacked profile, but the two observations are nevertheless compatible. We used the 10 000 bootstrap realizations to estimate the correlation of the two observables and then combine them into a signal constraint. The resulting 68% and 95% C.L. are shown as the orange and yellow filled contours in Fig. 3. The corresponding best fit is marked with a white cross and is given in the second row of Table 1. The result, r_L = 0.84 ± 0.20 is compatible with one, as for the stacked profile constraint. The value for β is fully driven by the stacked profile because the average matched filter flux does not constrain it.

The best fit for the stacked profile has 18 degrees of freedom (d.o.f.). corresponding to 20 radial bins minus two parameters. The best fit for the combined constraint is 19 d.o.f. corresponding to 20 radial bins plus 1 bin averaged matched filter flux minus two parameters. The reduced χ² is acceptable (χ²/d.o.f. = 1.15) for the profile only fit and shows some small tension between the two measurements for the combined case (χ²/d.o.f. = 1.40).

The dust profile corresponding to the combined best fit is shown as the orange dash-double dotted line in Fig. 2. The inverse-variance weighted matched filter flux for the profile (combined best fit) model is shown as the blue dash-dotted (orange dash double-dotted) line for SZ+dust in Fig. 4, to be compared to the SZ-only signal shown as the black dashed curve. The blue line is significantly (3.9σ) higher than the measurement in the 857 GHz band. In the 857 GHz band, the orange line is in good agreement with the data by construction.

Although the combined fit is performed using 857 GHz data only, the agreement at lower frequencies is good. This demonstrates that the galaxy populations of the De Zotti et al. (2018) model and their spectral energy distributions provide a satisfactory description of the frequency dependance of the dust emission of Planck clusters. The comparison between the spectral energy distributions of galaxies in the De Zotti et al. (2018) model and the spectral energy distribution of Planck clusters determined by Planck Collaboration XLIII (2016) is discussed in Appendix A.2.

We adopt the result of the combined fit (second row of Table 1) as our fiducial dust model. The corresponding model parameters are r_L = 0.84 in Eq. (5) and β = 1.29 in Eq. (10) (with α = γ = 1, c₅₀₀ = 3 fixed).

5. Impact on Planck cosmological results

Planck cosmological analyses with clusters (in particular, cluster extraction and cosmological constraints from cluster counts) do not take into account dust emission in clusters. The omission of this emission may possibly impact the cluster physical parameter recovery (cluster size and flux) and the survey completeness. In Sect. 5.1, we use our fiducial dust model built in Sect. 4 to study the effect of cluster dust emission on size and flux recovery. In Sect. 5.2, we calculate the effect of cluster dust emission on the Planck completeness and show the impact on cosmological parameter determination.

5.1. Cluster size and flux recovery

The Planck beams (FWHM ranging from 9.6 to 4.2 arcmin between 100 and 857 GHz) are larger than the typical cluster extent (1 arcmin), meaning that Planck provides weak constraints on cluster size. As a direct consequence, blind fluxes are only weakly constrained by the extraction tools. This problem is often referred as the “size-flux degeneracy” (see e.g., introduction of Sect. 7.2 of Planck Collaboration XXIX 2014).

The Planck collaboration noticed that blindly recovered cluster sizes are over estimated on average with respect to the sizes estimated from X-ray observations. This size over-estimation translates into an over-estimation in the blind flux relative to expectations based on the X-rays. For this reason, the Planck collaboration computed cluster flux fixing the size from X-ray measurements (Sect. 7.2.1 of Planck Collaboration XXIX 2014) or adopting a size-flux relation from X-ray as a prior to break the Planck size-flux degeneracy. The latter approach is used to derive the “mass proxy” $M_{500}^{\mathrm{Yz}}$ (Sect. 7.2.2 of Planck Collaboration XXIX 2014) that we use in this paper. However, this size overestimation is not present in millimeter simulations which include SZ as the only cluster emission and for which the simulated SZ profile perfectly matches the profile assumed in the extraction tool (Fig. 8 of Melin et al. 2006). In this section, we examine the over-estimation seen in the Planck data, looking to see if it is related to some additional cluster component, such as dust, or is linked to the profile assumed for the extraction.

We simulate 1091 clusters with the same masses and redshifts as the PSZ2. We inject them randomly into Planck frequency maps, outside the 85% survey mask³ to avoid contamination by Galactic dust and outside a PSZ2 cluster mask⁴ to avoid contamination by real clusters. We model the SZ emission using the universal pressure profile (UPP, Arnaud et al. 2010) or the Planck pressure profile (PlanckPP, Planck Collaboration Int. XI 2013). Although the PlanckPP is consistent with the UPP in the inner cluster regions (R < R₅₀₀), it is significantly more extended to larger radii (R₅₀₀ < R < 3R₅₀₀, see Fig. 4 left of Planck Collaboration Int. XI 2013). Thus, assuming the UPP for extracting clusters well described with a PlanckPP could possibly lead to a size over-estimation.

We then modeled the dust emission using our combined best fit described in Sect. 3. We simulate the four possible combinations (UPP and PlanckPP, with and without dust) and extract cluster size and flux using the MMF. For the cluster size, θ₅₀₀, we use a grid of 32 filter sizes equally spaced on a logarithmic scale and ranging from 0.94 to 35.31 arcmin. We searched for the cluster position as a maximum of S/N in a circle of radius 20 arcmin around the real or injected position and adopt the UPP in the MMF for the four cases. We then compared the recovered size and flux to their input values to see if we could identify the origin of the blind size and flux over-estimation found in the data.

Results are shown in Fig. 5 for the size. The top left panel shows the extraction at the location of the actual PSZ2 clusters. The recovered size ${\theta }_{500}^{\mathrm{blind,data}}$ is weakly constrained and is significantly biased high with respect to the size derived from the mass proxy ${\theta }_{500}^{\mathrm{Yz}}$. The mean (median) of the ratio of the two quantities is 1.35 (1.20). The thickness of the line encapsulates the 68% error on the mean (median) calculated with bootstrap. The inset shows the histogram of the distribution, which peaks above one. The recovered size ${\theta }_{500}^{\mathrm{blind,data}}$ is discretized and corresponds to the values adopted for our grid. One can also notice that the algorithm sometimes fails to recover cluster size and falls onto the grid limits.

The top right panel shows the extraction at the location of the injections with the UPP and without dust emission. As already noticed in Melin et al. (2006), although the recovered size ${\theta }_{500}^{\mathrm{blind,\; w/o\; dust}}$ is weakly constrained, it is significantly less biased with respect to the injected size ${\theta }_{500}^{\mathrm{Yz}}$. The mean (median) of the ratio is 1.12 (1.02) and the histogram peaks around unity. Again, the thickness of the line encapsulates the 68% error on the mean (median) calculated with bootstrap. We tested this bootstrap error on the mean (median) by performing ten injections of the PSZ2 clusters and in computing the standard deviation of the mean (median) values across the ten corresponding extractions. The standard deviation across the ten extractions is in very good agreement with the bootstrap error on a single extraction, with a value of 0.02 for both the mean and the median.

The result of adding dust emission to SZ emission in simulated clusters is shown in the bottom left panel for the UPP. The result is essentially identical to the UPP without dust. The mean (median) of the ratio is 1.10 (1.01) and the histogram peaks around unity.

The effect of changing the SZ profile to the PlanckPP is shown in the bottom right panel, without dust emission. The recovered size is overestimated with a mean (median) ratio equal to 1.25 (1.16), close to the value observed in the actual data. The result of including dust with the PlanckPP is not shown: it is almost identical to the PlanckPP without dust as for the UPP.

This test demonstrates that dust emission has no impact on cluster size estimate with the MMF. It also indicates that the size overestimation may find its origin in the profile mismatch between actual clusters and the UPP. Indeed, adopting the PlanckPP in the simulations and extracting clusters using the UPP reproduces the bias observed in the data. The dispersion of the histogram in the bottom right panel (PlanckPP without dust) is smaller than that of the histogram in the top left panel (actual data). This could be due to the dispersion in the actual pressure profiles which is not included in the simulations, the PlanckPP being the average value.

The size overestimation shown in Fig. 5 directly impacts the flux. Thus the flux estimation depends on the profile assumed for the injected model. The results are shown in Appendix B (Fig. B.1). As for the size, the flux is overestimated for the actual PSZ2 and the PlanckPP case (without and with dust), but not overestimated for the UPP (without and with dust).

Finally, we examined the impact of dust emission on the PSZ2 flux estimation when fixing both cluster position and size. Results are shown in Fig. 6. The impact of dust emission on flux estimation is negligible (<1%) for bright clusters (Y_z > 10⁻³ arcmin²). The bias due to dust then increases from <1% to ∼2% with decreasing flux from Y_z = 10⁻³ arcmin² to 5 × 10⁻⁴ arcmin².

5.2. Planck cluster completeness and cosmological constraints

We now investigate the impact of cluster dust emission on the Planck cluster catalog completeness, and then on the measurement of cosmological parameters from cluster counts.

We randomly drew cluster redshifts and masses from the Tinker et al. (2010) mass function, model SZ emission with the UPP and normalize the flux using the Y–M relation from Arnaud et al. (2010). We adjusted the mass bias to 1 − b = 0.65 (see Eq. (7) of Planck Collaboration XX 2014) to match the model counts in our adopted cosmology to the observed counts. We injected the clusters into the Planck frequency maps at random locations outside the 65% cosmological mask⁵ and outside the same PSZ2 cluster mask as in Sect. 5.1 to avoid contamination by real clusters. We considered two cases: with and without inclusion of dust emission in addition to the SZ signal. We then used the MMF algorithm to extract clusters blindly, following the same procedure as for the Planck analyses (Planck Collaboration VIII 2011; Planck Collaboration XXIX 2014; Planck Collaboration XXVII 2016). We perform ten such injections for each case. In order to improve the statistics at high redshift, we also simulated, injected and extracted ten additional independent catalogs containing only clusters at z > 0.5, but with ten times higher density.

We compared the counts for recovered clusters with and without inclusion of dust emission. Results are shown in Fig. 7. The solid black line gives the total number of detected clusters with dust in the 10 + 10 simulations divided by the total number of detected clusters without dust, as a function of redshift. The red band corresponds to the standard deviation of 10 000 bootstraps over the 10 + 10 simulations.

The dust emission significantly impacts the Planck survey completeness over the redshift range [0.3–0.8], with a loss of ∼9% of clusters in the [0.5–0.8] range. The [0–0.3] redshift range is only affected by < 2%, as is the [0.8–0.9] bin. The [0.9–1] bin may present a small excess of detected clusters due to dust (+9%), although the value is not statistically significant (the bootstrap error is 4.2% in this bin). We show in Appendix D that our dust model depends only weakly on cosmological parameters and, for simplicity, we adopt the curve shown in Fig. 7 as the correction factor to apply to predicted counts before being associated with the observed counts in the Planck likelihood.

We reran the Markov chains for the N(z) likelihood, correcting the completeness from the effects of the dust, over the full redshift range [0–1], and on the two distincts ranges [0–0.2] and [0.2–1], reproducing what was done in Planck Collaboration XXIV (2016). The results are shown in the lefthand panel of Fig. 8, presented in the same format as Fig. 7 of Planck Collaboration XXIV (2016) to ease comparison.

Over the entire redshift range [0–1] and also the high redshift range [0.2–1], the change in the contours is negligible. There is a change, on the other hand, for the low redshift range [0–0.2]; in particular, the Ω_m posterior loses its bimodality. This could be due to a lack of convergence of the original chains or to the fact that the low redshift likelihood is unstable.

To decide between these two possibilities, we reran the original Planck likelihood (i.e., without any dust correction to the completeness) to a higher level of convergence. The results are shown in the righthand panel of Fig. 8. The bimodality of the Ω_m posterior remains. We note that the low redshift contours in the righthand panel of Fig. 8 and Fig. 7 of Planck Collaboration XXIV (2016) differ slightly. Specifically, the maximum of the Ω_m posterior is now the high Ω_m solution, while it was the low Ω_m solution in Planck Collaboration XXIV (2016). We conclude that the low-z likelihood is somewhat unstable. This is supported by the change in contour shape with increasing convergence, as just noted, and also by the fact that the dust correction in the first two redshift bins is small (0.994 for 0 < z < 0.1 and 0.988 for 0.1 < z < 0.2).

Despite this change in the low z likelihood, the two panels of Fig. 8 are remarkably similar. This demonstrates that taking dust contamination into account in the analysis does not significantly change the preferred cosmological parameters, and does not ease the tension with the primary CMB.

6. Conclusion

We have modeled dust emission in galaxy clusters at millimeter wavelengths using the model by De Zotti et al. (2018), which we augmented by stacking PSZ2 clusters. The model now gives the shape of the dust profile and a normalization for the dust emission based on the Planck 857 GHz channel. We used this model to simulate clusters that we injected into the Planck maps. We then assessed the impact of dust emission on Planck cluster results, finding that:

• Dust emission is not responsible for the cluster size over-estimation seen in the real data.

• The size over-estimation is plausibly caused by a mismatch, in the external regions, between the true cluster pressure profiles and the UPP adopted in the cluster extraction tool.

• When fixing cluster size and position, dust emission biases Planck cluster flux measurements low at only the 1–2% level.

• Dust emission impacts the completeness of the cluster cosmology catalog over the redshift range [0.3–0.8], with a maximum loss of ∼9% of clusters between z = 0.5 and z = 0.8.

• This cluster loss has a negligible effect on cosmological parameter estimation. Taking dust contamination into account in the Planck cluster cosmology analysis does not help to ease the tension with the primary CMB.

Our constraints on the cluster dust emission model are general, and the calibrated model can be used to evaluate the impact of dust emission on other SZ surveys. This is of particular interest for the highly sensitive next generation SZ cluster surveys such as the Advanced ACT (De Bernardis et al. 2016), SPT-3G (Benson et al. 2014), the Simons Observatory⁶, CMB-S4 (Abazajian et al. 2016), CORE (Delabrouille et al. 2018; Melin et al. 2018), and PICO (Hanany & Inflation Probe Mission Study Team 2018).

---

Acknowledgments

The authors would like to thank the anonymous referee for useful comments that helped to clarify some important aspects of this work. JBM would like to thank G. Hurier for the detailed discussion about Planck Collaboration XXIII (2016). A portion of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. GDZ gratefully acknowledges financial support from ASI/INAF agreement n. 2014-024-R.1 and from the agreement ASI/Physics Department of the university of Roma–Tor Vergata n. 2016-24-H.0. Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package.

References

Appendix A: Comparison with published results

In this appendix, we compare our augmented De Zotti et al. (2018) model to previously published results. We discuss the comparison with the PSZ2 dust emission estimate from Planck Collaboration XXIII (2016) in Appendix A.1, and comparison between the dust SED from Planck Collaboration XLIII (2016) with the SED from Cai et al. (2013; used in the De Zotti et al. 2018 model) in Appendix A.2.

A.1. Comparison to PSZ2 dust emission from Planck Collaboration XXIII (2016)

In Fig. A.1, we reproduce Fig. 7 from Planck Collaboration XXIII (2016). Instead of displaying the total flux density of the stack, we prefer to show the average value, so our y-axis must be multiplied by 1091(number of clusters in the analysis) to match the y-axis from Planck Collaboration XXIII (2016). The central values (black diamonds) are in good agreement with the values found in Planck Collaboration XXIII (2016; displayed in our figure as red filled circles and shifted by +10 GHz for clarity), but our error bars are much larger. This could plausibly be due to the different methods used to estimate errors.

We estimated our errors from the standard deviation of a bootstrap resampling of the sample of 1091 clusters, while Planck Collaboration XXIII (2016) computes the standard deviation at 1000 random locations on the sky. This second method does not capture the intrinsic variation of dust emission across the cluster population. Therefore, our fitted models (blue and orange lines) are fully consistent with our data points and error bars, but are significantly below the model proposed in Planck Collaboration XXIII (2016) and shown as the red short dashed line. The black dashed line shows the contribution of the SZ signal only.

The difference between Fig. 4 and Fig. A.1 comes from the difference in the measurement and averaging procedure. For Fig. 4, the signal extraction is performed on individual clusters within an area of radius 5 × R₅₀₀ using matched filters and assuming the profile from Arnaud et al. (2010). With this template, the flux within 5 × R₅₀₀ is then converted to the flux within a sphere of radius R₅₀₀. Individual cluster fluxes are combined using an inverse-variance weighted average. For Fig. A.1, on the other hand, the signal is obtained from raw stacked maps and the error bars determined via bootstrap resampling. The flux in this case is estimated within a 20 arcmin radius aperture.

A.2. Comparison of the SEDs from Planck Collaboration XLIII (2016) and Cai et al. (2013)

In Fig. A.2, we compare the SED determined on PSZ2 clusters in Planck Collaboration XLIII (2016) to the SED from Cai et al. (2013) used in the De Zotti et al. (2018) model. The frequency dependence of the warm and cold SEDs of Cai et al. (2013) are similar to the SED from the PSZ2 over the range covered by Planck (100–857 GHz, or equivalently from 350 to 3000 μm i.e., log λ between 2.54 and 3.48). This explains why our model, although only adjusted to the 857 GHz data, also provides a good match at lower Planck frequencies in Fig. 4.

Appendix B: Blind flux estimation of PSZ2 clusters

In this appendix, we compare the SZ signal extracted blindly to the injected signal for simulated PSZ2 clusters. Fig. B.1 is equivalent to Fig. 5, but for the flux instead of the size. The figure shows that the over-estimation of the blind flux is not due to the dust emission in clusters. It likely finds its origin in the mismatch between the profile assumed for cluster extraction (UPP) and the actual cluster profile (closer to the PlanckPP), as already shown in Fig. 5 for cluster size and discussed in Sect. 5.1.

Appendix C: Stacked maps

We show Planck maps stacked at the PSZ2 (Fig. C.1) positions and at random positions (Fig. C.2). In Fig. C.1, the negative part of the SZ effect is clearly visible at 100 and 143 GHz, and the dust contribution mixed with the increment of the SZ emission is visible at frequencies above 217 GHz. No significant emission is found in Fig. C.2.

Appendix D: Cosmology dependence of the effect of dust emission on Planck completeness

The effect of dust emission on the Planck completeness may depend on cosmology, in particular because the dust model fit performed in Sect. 4 may depend on the assumed cosmological parameters. Expressing the completeness as a function of redshift introduces an additional dependance on cosmology. It is technically feasible to express it as a function of cluster flux and size, as in the Planck analyses, to avoid this latter cosmological dependance, but this would require significant additional computing time to run more simulations to build the two dimensional quantity. The major difficulty would be to assess the dependence on cosmology: we would need to Monte Carlo the whole analysis chain (fit for the dust model, simulations, injections, extractions) on each set of cosmological parameters, which would require some unmanageable computing time. Running the full analysis takes about two weeks for a single cosmology.

In order to test the cosmology dependance of the effect of dust on Planck completeness, we thus performed a second full analysis and changed the value for Ω_m to 0.4 and Ω_Λ to 0.6, while keeping the other parameters fixed to the Planck ΛCDM cosmology. This model is located far from our fiducial Planck ΛCDM cosmology in the 95% C.L. region of the Planck cluster cosmological constraints (Fig. 7 of Planck Collaboration XXIV 2016). The impact of the adopted model on completeness is shown in Fig. D.1 and the ratio between the black line of this figure and the one from Fig. 7 is shown in the inset. The change in the effect of dust emission between the two sets of cosmological parameters is weak (<4% in the full Planck cluster redshift range [0−1] as shown in inset). We thus adopted the curve from Fig. 7 to correct our predicted cluster counts for all sets of cosmological parameters.

All Tables

All Figures