Dust in galaxy clusters: Modeling at millimeter wavelengths and impact on Planck cluster cosmology

We have examined dust emission in galaxy clusters at millimeter wavelengths using the Planck $857 \, {\rm GHz}$ map to constrain the model based on Herschel observations that was used in studies for the Cosmic ORigins Explorer (CORE) mission concept. By stacking the emission from Planck-detected clusters, we estimated the normalization of the infrared luminosity versus mass relation and constrained the spatial profile of the dust emission. We used this newly constrained model to simulate clusters that we inject into Planck frequency maps. The comparison between clusters extracted using these gas+dust simulations and the basic gas-only simulations allows us to assess the impact of cluster dust emission on Planck results. In particular, we determined the impact on cluster parameter recovery (size, flux) and on Planck cluster cosmology results (survey completeness, determination of cosmological parameters). We show that dust emission has a negligible effect on the recovery of individual cluster parameters for the Planck mission, but that it impacts the cluster catalog completeness, reducing the number of detections in the redshift range [0.3-0.8] by up to $\sim 9\%$. Correcting for this incompleteness in the cosmological analysis has a negligible effect on cosmological parameter measurements: in particular, it does not ease the tension between Planck cluster and primary cosmic microwave background cosmologies.


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(Alberts et al. , 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 (2013b) 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 (2016b) and Planck Collaboration (2016e) 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 (2016b) 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 (2016e) 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 2016d) 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 2016a): h = H 0 /(100 kms −1 Mpc −1 ) = 0.6774, Ω m = 1 − Ω Λ = 0.3089, Ω b = 0.0485976, n s = 0.9667.

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 from http://pla.esac.esa. int/pla/. 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 2016d), 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 2011; Melin et al. 2012;Planck Collaboration 2014b, 2016d. 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 2014a). This allows us to compute the mass proxy of each object, M Yz 500 (and the associated cluster size θ Yz 500 , since the redshift is known). More detail on the method of computing the mass proxy is given in Sec. 7.2.2 of Planck Collaboration (2014b). We thus have 1091 objects with position, redshift z, mass M Yz 500 and size θ Yz 500 .

Dust modeling
Our baseline dust model is built on Herschel observations of field and cluster galaxies (Alberts et al. 2014(Alberts et al. , 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. 2017;De Zotti et al. 2016). 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 1 , we compute, for each redshift, the comoving infrared luminosity density, Ψ i (z), contributed by each galaxy population i: for i = cold, warm, spheroidal. log is the base-10 logarithm. We then compute the total comoving infrared field luminosity density 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) 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 where V comoving = (1 + z) 3 × 4π 3 R 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 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 We then compute, for each redshift, the fraction p i (z) of the total luminosity contributed by the galaxy population i as The luminosity from each population is then given by 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 with D L (z) the luminosity distance. This model is identical to the model adopted in Sect Figure 1 (bottom) shows the predicted dust flux density integrated over the Planck 857 GHz bandwidth, S 500 , as a function of redshift for a cluster of mass M 500 = 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.

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): where α, β and γ are, respectively, the intermediate, external and central slopes, r s = R 500 /c 500 is the scale radius and c 500 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 200 = 4.61 ± 0.12, which we convert to c 500 = 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 500 = 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).

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 (2016b). 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 θ Yz 500 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 +0.34 −0.33 , is compatible with one, the value determined from Herschel data. The slope parameter, β = 1.36 +0.39 −0.30 , 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.   Table 1: Best fit values for the normalization, r L , of the infrared L 500,tot − M 500 relation and for the external slope, β, of the spatial profile of the dust emission. Errors are 68% C.L. Fig. 3: Contours at 68% and 95% C.L. on the normalization, r L , of the infrared L 500,tot − M 500 relation and on the external slope, β, of the spatial profile of the dust emission. Constraints are obtained from the stacked profile (blue) and from the inversevariance weighted average matched filter flux (green), both in the Planck 857 GHz band. The combined constraint is shown as filled orange and yellow contours. The blue and white crosses shows the best value for the profile and combined fits, respectively.

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 θ Yz 500 , 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.

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 χ 2 is acceptable (χ 2 /d.o. f.=1.15) for the profile only fit and shows some small tension between the two measurements for the combined case (χ 2 /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 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 500 = 3 fixed).

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.

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 2014b).
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 2014b) 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 Yz 500 (Sect. 7.2.2 of Planck Collaboration 2014b) 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 overestimation 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 2 to avoid contamination by Galactic dust and outside a PSZ2 cluster mask 3 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 2013a). Although the PlanckPP is consistent with the UPP in the inner cluster regions (R < R 500 ), it is significantly more extended to larger radii (R 500 < R < 3R 500 , see Fig. 4 left of Planck Collaboration 2013a). 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, θ 500 , 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 θ blind,data 500 is weakly constrained and is significantly biased high with respect to the size derived from the mass proxy θ Yz 500 . 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 θ blind,data 500 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 θ blind, w/o dust 500 is weakly constrained, it is significantly less biased with respect to the injected size θ Yz 500 . 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 −3 arcmin 2 ). The bias due to dust then increases from < 1% to ∼ 2% with decreasing flux from Y z = 10 −3 arcmin 2 to 5 × 10 −4 arcmin 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 2014a) 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 4 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 2011, 2014b, 2016d. 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 (2016c). The results are shown in the lefthand panel of Fig. 8, presented in the same format as Fig. 7 of Planck Collaboration (2016c) 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 The value for the mean (median) is 1.12 (1.02). There is no strong overestimation of the size as on the actual data. Bottom left: Same as top right but adding the dust component based on the best combined fit (white cross in Fig. 3). The mean (median) is 1.10 (1.01). The impact of the dust component on the size estimation is negligible. Bottom right: Same as top right but using the Planck pressure profile instead of the UPP to simulate clusters. The mean (median) is 1.25 (1.16). The blind sizes are overestimated as for the Planck data, although the histogram in the inset is less dispersed. The dotted line in all four panels is the equality line.
Ω m posterior remains. We note that the low redshift contours in the righthand panel of Fig. 8 and Fig. 7 of Planck Collaboration (2016c) 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 (2016c). 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.

Conclusion
We have modeled dust emission in galaxy clusters at millimeter wavelengths using the model by De Zotti et al. (2016), 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 overestimation 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 to 2% level. , when fixing position and size for clusters simulated with and without dust as a function of injected SZ flux Y z . The red line is the raw mean value. The impact of dust emission is negligible for bright clusters (Y z > 10 −3 arcmin 2 ) and increases to ∼ 2% with decreasing flux down to Y z = 5 × 10 −4 arcmin 2 .  (Fig 7). Shifts in cosmological parameters (black curves) are negligible with respect to the case when dust is not taken into account (right panel). Right: Cosmological parameters from the N(z) Planck likelihood without any dust correction to the completeness. This figure was obtained with the same likelihood as the original analysis ( Fig. 7 of Planck Collaboration 2016c), but the convergence of the chains is higher.
We estimated our errors from the standard deviation of a bootstrap resampling of the sample of 1091 clusters, while Planck Collaboration (2016b) 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 (2016b) 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 500 using matched filters and assuming the profile from Arnaud et al. (2010). With this template, the . The y-axis shows the average flux density of a cluster, so it needs to be multiplied by 1091 (the number of clusters in the analysis) to give the total flux density of the stack. Our data points are shown as black diamonds, and the points from Planck Collaboration (2016b) as red filled circles (shifted by +10 GHz for clarity). Our models are shown in blue and orange, and the model from Planck Collaboration (2016b) is shown in red. The black dashed line gives the contribution from the SZ signal alone. flux within 5 × R 500 is then converted to the flux within a sphere of radius R 500 . 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.  Cai et al. (2013). We have added the best fit spectrum from Planck Collaboration (2016e) as the solid black line. It has been normalized so that the warm SED and the Planck SED have the same integrated luminosity (0.348 L ) between 100 and 3000 GHz (or equivalently between 100 and 3000 µm i.e., log λ between 2 and 3.48).

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 2016c). 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.  The maps are 2 × 2 deg 2 in units of Jy/arcmin 2 . We have adopted the same color scales as in Fig. C.1. The white circle is located at the stack center and is 20 arcmin in radius.