MEASURING THE REDSHIFT DEPENDENCE OF THE COSMIC MICROWAVE BACKGROUND MONOPOLE TEMPERATURE WITH PLANCK DATA

Our cluster sample contains 623 clusters outside WMAP Kp0 mask. It was created by combining the ROSAT-ESO Flux Limited X-ray catalog (Böhringer et al. 2004) in the southern hemisphere, the extended Brightest Cluster Sample (Ebeling et al. 1998, 2000) in the north, and the Clusters in the Zone of Avoidance (Ebeling et al. 2002; Kocevski et al. 2007) sample along the Galactic plane. All three surveys are X-ray-selected and X-ray flux limited. A detailed description of the creation of the merged catalog is given in Kocevski & Ebeling (2006). The position, flux, X-ray luminosity, and angular extent of the region containing the measured X-ray flux were determined directly from the ROSAT All Sky Survey (RASS). All clusters have spectroscopically measured redshifts. The X-ray temperature was derived from the L_X − T_X relation of White et al. (1997). The central electron densities and core radii were derived by fitting to the RASS data a spherically symmetric isothermal β-model (Cavaliere & Fusco-Femiano 1976) convolved with the RASS point-spread function. The β was fixed at the canonical value of 2/3 to reduce the dependence of the β-model parameters with the choice of radius over which the model is fit. These data are sufficient to compute the Comptonization parameter at the center of the cluster.

Atrio-Barandela et al. (2008) compared the TSZ predicted from the X-ray data with the signal present in WMAP three year data and found it to be in good agreement within the X-ray emitting region, where the β-model is a good description of the electron distribution. In the cluster outskirts, the predicted TSZ signal was systematically higher than the measured value. The latter was consistent with the Komatsu & Seljak (2002) profile, where baryons are in hydrostatic equilibrium within a dark matter halo well described by a Navarro–Frenk–White profile (hereafter NFW, Navarro et al. 1997), as expected in the concordance ΛCDM model. More recently, Nagai et al. (2007) proposed a scaled three-dimensional electron pressure profile p(x) = P_e(r)/P₅₀₀ based on a generalized NFW profile:

$$\begin{equation} p(x) = \frac{P_0}{(c_{500}x)^\gamma [1+(c_{500}x)^\alpha ]^{(\beta -\gamma)/\alpha }}, \end{equation} \tag{ 3 }$$

where (γ, α, β) are the central, intermediate, and outer slopes, c₅₀₀ characterizes the gas concentration, and x = r/R₅₀₀ is the distance from the center of the cluster in units of the radius at which the average density of the cluster is 500 times the critical density. Arnaud et al. (2010) derived an average cluster pressure profile from observations of a sample of 33 local (z < 0.2) clusters, scaled by mass and redshift with

$$\begin{eqnarray} &&[P_0,c_{500},\gamma,\alpha,\beta ] \nonumber\\ &&\qquad= \big[8.403 h_{70}^{-3/2},1.177,0.3081,1.0510,5.4905\big]. \quad \end{eqnarray} \tag{ 4 }$$

Later, Plagge et al. (2010) showed these parameters to be consistent with the SZ measurements of 15 massive X-ray clusters observed with the South Pole Telescope (Plagge et al. 2010). We determine the scale R₅₀₀ using the Böhringer et al. (2007) relation:

$$\begin{eqnarray} R_{500} &=& \frac{(0.753\pm 0.063) h^{-1} \rm { \ Mpc}}{h(z)}\nonumber\\ &&\times\left(\frac{L_X}{10^{44} h^{-2}\rm {\ erg\ s}^{-1}}\right)^{0.228\pm 0.015}. \end{eqnarray} \tag{ 5 }$$

To test the effect of the cluster profile on the final results, we construct y-maps from the X-ray cluster catalog (1) using the universal pressure profile of Equation (3) with the parameters given in Equation (4) and (2) using the isothermal β = 2/3 model. The Comptonization parameter is computed by integrating the electron pressure profile along the line of sight. Clusters are assumed to be spherically symmetric and extending up to R₂₀₀, the scale where the cluster overdensity reaches 200 times the critical density. To determine the effect of the cluster profile, the central value of y_c is assumed to be the same for both the β-model and the universal pressure profile. Finally, the cluster templates are convolved with the corresponding antenna beams (see Table 1). In Figure 3(a) we show the pressure profile integrated along of line of sight for the β = 2/3 (solid line) and universal pressure (dashed line) profiles convolved with the antenna of the 44 GHz map. Since the latter profile diverges at x = 0, we restrict the integration of Equation (3) to θ  ⩾  10⁻²θ₅₀₀. The cluster is located at z = 0.094 of M₅₀₀ = 2.4  ×  10¹⁴ h⁻¹ M_☉ and R₅₀₀ = 746 h⁻¹ Kpc. For illustration, in Figure 3(b) we plot the value of Comptonization parameter y_c at the center of all the clusters in our proprietary cluster catalog, derived using the measured X-ray parameters, as a function of cluster mass. The solid line represents the linear regression fit to the data. The central Comptonization parameter scales as y_c = 24.5(M₅₀₀/10¹⁴ h⁻¹ M_☉)^1.35.

Zoom In Zoom Out Reset image size

As an alternative, we also use the low-redshift all-sky maps and the associated galaxy cluster catalogs of the hydrodynamic diffuse and kinetic SZ simulations included in the pre-launch Planck Sky Model. The simulations are fully described in Delabrouille et al. (2012). The catalogs contain cluster positions, mass, and radius for an overdensity contrast of 200 times the critical density. The maps contain the integrated SZ signal up to z  ≃  0.25, computed from a combination of full hydrodynamic simulations using the box stacking method described in A. Valente et al. (2012, in preparation). According to this method, the universe around the observer is generated in concentric layers, each with a comoving thickness of 100 h⁻¹ Mpc, using the outputs of hydrodynamic simulations with periodic boundary conditions. The light-cone integrations of the TSZ and KSZ signals are carried out using the formulae in da Silva et al. (2000, 2001). A total of seven layers were constructed, up to z = 0.25. The innermost layer includes the local constrained simulation of Dolag et al. (2005), whereas all the other layers were produced from gas snapshots of the ΛCDM simulation in De Boni et al. (2011). Both these simulations include explicit treatment for gas cooling, heating by UV, star formation, and feedback processes.

The y-map constructed from the X-ray-selected clusters assumes clusters to be spherically symmetric and relaxed while the TSZ and KSZ templates constructed from the hydro-simulation contain clusters with different dynamical state (relaxed, merging systems, etc.), shape, and ellipticity. Also, since the latter are constructed integrating the signal along the line of sight, the projection effects due to low mass clusters and groups are included. Therefore, these templates are very well suited to study the effect of all these systematics and of the KSZ component in the determination of α. For a more realistic comparison, we select 623 clusters from the simulation according to the measured selection function of the X-ray cluster sample. In Figure 4(a) we plot the mass and in Figure 4(b) the redshift distribution of all clusters in our simulation (solid line) that fulfill the selection criteria. The dashed line shows the same distributions of the X-ray clusters. For a better comparison, the histogram of the largest amplitude was normalized to unity. The main difference between the two samples is that there are 22 clusters in our proprietary cluster catalog that have redshifts larger than z = 0.25, the redshift of the last layer constructed from the simulation.

Zoom In Zoom Out Reset image size

The y-maps described above are multiplied by G(ν, α = 0) to generate TSZ templates and convolved with the antenna beam. Since the cluster catalog contains X-ray temperatures for all clusters, in the template generated from the sample of real clusters we included relativistic corrections using the analytic formulae derived by Nozawa et al. (2000) that are accurate for clusters with temperature T_X ⩽ 25 KeV at all the frequencies measured by Planck. A KSZ template was added to the hydrodynamical but not to the X-ray-selected cluster template since their peculiar velocity are not available. Noise maps were constructed assuming the Blue Book noise levels of Table 1. We model the noise as homogeneous and uncorrelated white noise since at the frequencies 44–353 GHz the 1/f is both small and does not affect the angular scales subtended by clusters, ℓ ⩾ 500.

To construct the final maps for the two types of data (A) and (B) we add different contributions. To make the simulations as realistic as possible we need to consider CMB and foreground residuals, noise inhomogeneities, and deviations from white noise behavior. For instance, when averaged over the cluster extent, CMB residuals do not scale down as N^−1/2_pixel and, together with foreground residuals, they can induce systematic shifts in the Comptonization parameter. The effect will be dependent on both amplitude and shape of the residual power spectrum. CMB residuals dominate over the noise up to ℓ ∼ 2000 on the 70–353 GHz channels (see Leach et al. 2008, Figure 6). Therefore, one can expect that these residuals would dominate over noise inhomogeneities or 1/f contributions and so we will not consider the latter in our simulations. With respect to foreground residuals, they change with frequency, have different origin, different amplitude and are irregularly distribute in the sky. Further, only low-resolution (∼1°) templates are available for synchrotron and free–free emission (Leach et al. 2008) so the behavior of the foregrounds, not just foreground residuals, on scales <1° relevant for the study of clusters, cannot be modeled accurately. However, Galactic emissions are not associated with clusters, their effect on the measured value of α can be analyzed by studying different cluster subsamples in different regions of the sky. Other foregrounds associated with clusters could have a more damaging effect. Diffuse and point-like radio emission can be found in cluster cores. Their physical origins are very different though and both are rare, in particular among the very X-ray luminous clusters that dominate the signal discussed by us here (Brunetti et al. 2007). Specifically, radio halos are found only in merging clusters (Cassano et al. 2010). Radio emission from individual cluster galaxies is usually limited to the brightest cluster galaxy and, opposite to the case of diffuse emission, observed almost exclusively in relaxed systems. Like for other foregrounds, the contribution of both effects could be estimated by analyzing different cluster subsamples, binning by redshift or by galactic latitude. In any case, all these effects can only be accurately modeled once the data become available.

In comparison, CMB residuals represent an extreme case of foreground: they are the same at all frequencies, except by the change in amplitude due to the dilution by the beam, and they are distributed everywhere in the sky. In this context, CMB residuals are the most difficult to separate so they will be the only ones included in our simulations. Since different component separation methods will leave different residuals, we adopted the filter described in Kashlinsky et al. (2009) to create our CMB residual template. This filter is designed to remove the primary CMB fluctuations from the concordance ΛCDM model by minimizing the mean squared deviation of the CMB measurements from noise 〈(δ_CMB − noise)²〉, and it generates CMB residuals that are homogeneously distributed in the sky. In this respect we are being conservative; the true residuals are likely to be smaller away from the galactic plane (see Planck HFI Core Team 2011b, Figure 39), fact that could permit to estimate the effect of the residuals by analyzing cluster subsamples selected according to galactic latitude.

To summarize, we constructed two sets of maps, according to the specific data sets: (A) CMB subtracted using a template, (B) CMB subtracted exactly using the 217 GHz channel. In simulation (A) six maps are constructed, adding to the TSZ template CMB residuals and white noise. The noise is generated according to the parameters given in Table 1. In all the simulations the residual CMB map is kept fixed and only the noise varies. In this form we can quantify the bias introduced by the residuals. We took residual maps with amplitude 〈(ΔT_{CMB, res})²〉^1/2 = 1, 10 μK. For the y-map constructed from real clusters, only the TSZ signal was included. The temperature anisotropy at each pixel is ΔT_A(ν) = y_cG(ν, 0) + ΔT_TSZ ± σ^A_{noise, ν} ± ΔT_CMBres. For the y-map derived from simulations, we also included a KSZ component. In simulation (B) six maps are constructed adding a cosmological CMB signal and noise to the y-map. The 217 GHz map is used to subtract the cosmological CMB and KSZ signals exactly. Therefore, only five different maps are available for the analysis. We checked the final maps had a power spectrum that was a pure white noise, with a slightly larger variance σ^B_{noise, ν}, sum of the original map plus the noise of the degraded 217 GHz map. At each pixel, the temperature anisotropy is ΔT_B(ν) = y_cG_{[ − 217 GHz]}(ν, 0) ± σ^B_{noise, ν}.

In this method we take the ratio of the temperature anisotropy at two different frequencies. Maps are brought to a common resolution before taking the ratio. Since the ratio of two Gaussian random variables is not a Gaussian, we estimate α using the approach introduced by Luzzi et al. (2009). If the average temperature anisotropy at two different frequencies i = 1, 2 at one cluster location j are denoted by δ_j(ν_i) = 〈ΔT_{j, I}(ν_i)〉, and they are measured with Gaussian errors σ_i then for each cluster j, the probability P_j(R_I) that the ratio R_I(ν₁, ν₂) = δ₁/δ₂ of two temperature anisotropies is in the range (r, r + dr) is given by

$$\begin{eqnarray} &&P_j(R_I(\nu _1,\nu _2,\alpha))=\frac{1}{2\pi \sigma _1\sigma _2}\nonumber\\ &&\quad\times\int _{-\infty }^{\infty }\!\! x\exp \bigg({-}\bigg[\frac{(x-\delta _1)^2}{2\sigma _1^2}+ \frac{(xR_I(\nu _1,\nu _2,\alpha)-\delta _2)^2}{2\sigma _2^2}\bigg]\bigg)dx.\nonumber\\ \end{eqnarray} \tag{ 6 }$$

The most likely value of α is then given by the maximum of the likelihood function:

$$\begin{equation} \log {\cal L}=\sum _{\nu _1,\nu _2}\sum _{j=1}^{N_{cl}} \log [P_j(\nu _1,\nu _2,\alpha)]. \end{equation} \tag{ 7 }$$

Luzzi et al. (2009) demonstrate that the ratio of the two distributions is biased; the ratio is dominated by the error on the denominator. For this reason, we use the denominator to be the measurement with the smallest fractional errors in order to minimize the bias, i.e., the 100 and 143 GHz channels. In total we have nine different ratios.

Alternatively, we can fit the TSZ signal of each cluster to the spectral dependence of Figures 1(b) and 2(b). In this case, errors are Gaussian distributed and the likelihood function is given by

$$\begin{equation} -2\log {\cal L}=\sum _\nu \sum _{j=1}^{N_{cl}} \left[\frac{\langle \Delta T_I(\nu)\rangle -\bar{y}_cG_I(\nu,\alpha)}{\sigma ^I_{{\rm noise},\nu,j}}\right]^2, \end{equation} \tag{ 8 }$$

where σ^I_{noise, ν, i} is the error at each cluster location that includes noise and CMB residuals (simulation (A)) or just noise (simulation (B)).

In this method, since we do not measure G(ν, α) directly but $\Delta T(\hat{n})=T_0y_cG(\nu,\alpha)$, we need an independent determination of y_c using X-ray data, introducing another complication: different frequencies have different resolutions and cluster anisotropies are diluted differently by the antenna beam. The amplitude of the effect depends on the cluster profile and angular extent but does not depend on the scaling of the TSZ signal with redshift, G(ν, α). To show this effect, in Figure 5 we represent the effect for two different clusters. Figure 5(a) corresponds to a cluster at redshift z = 0.218, with mass M₅₀₀ = 3.64 × 10¹⁴ h⁻¹ M_☉ and size 94 while in Figure 5(b) the cluster is located at z = 0.058 with mass M₅₀₀ = 7.7 × 10¹⁴ h⁻¹ M_☉ and size 42'. The dilution effect is most noticeable at 44 GHz since this channel has the smallest resolution.

Zoom In Zoom Out Reset image size

For clusters drawn from a simulation, its size, ellipticity, and profile are well known. For such clusters, the deconvolution factor F can be determined exactly by comparing the average Comptonization parameter before, 〈y_c〉, and after, 〈y_c*B(ν)〉, convolving with the antenna beam B(ν): F = 〈y_c〉/〈y_c*B(ν)〉. This factor would be different for resolved and unresolved clusters and would depend on the cluster profile and redshift. For illustration, in Figure 6(a) we represent F for a subset of 110 clusters in the mass range M₅₀₀ = 5–6 × 10¹⁴ h⁻¹ M_☉. In Figure 6(b), we plot the deconvolution factors for all clusters in our simulation with masses M₅₀₀  ⩾  10¹⁵ h⁻¹ M_☉. Solid black circles correspond to the 353 GHz frequency and open squares to 44 GHz. All clusters are resolved at 353 GHz. At 44 GHz clusters with redshift z  ⩾  0.08 are unresolved. For clarity, the resolved clusters are not shown. The solid straight lines correspond to the linear regression fit, F_lin, to the deconvolution factor, F, different for each cluster mass range and channel. Arrows indicate the deconvolution factor of the clusters plotted in Figures 5(a) and (b).

Zoom In Zoom Out Reset image size

If ΔF is the rms dispersion of the true deconvolution values F around F_lin for each mass bin and antenna, then F = F_lin ± ΔF. We used the y-map computed with clusters drawn from a numerical simulation to compute the deconvolution factors F_lin and its uncertainty ΔF in three mass bins of equal number of clusters: M₅₀₀ ⩽ 2 × 10¹⁴ h⁻¹ M_☉, M₅₀₀ = 2–3.6 × 10¹⁴ h⁻¹ M_☉, and M₅₀₀ ⩾ ×3.6 × 10¹⁴ h⁻¹ M_☉. These deconvolution factors were later used to deconvolve the effect of the beam on the sample of X-ray-selected clusters. If for a real cluster we use F_lin instead of the (unknown) true factor F, the deconvolved signal (ΔT_TSZ*B)F_lin would differ from the true signal ΔT_TSZ by an amount (ΔT_TSZ*B)ΔF. This uncertainty is uncorrelated with the instrumental noise at the cluster location and can be included in the likelihood analysis of Equation (8) by adding it in quadrature: σ²_{tot, i} = σ_{noise, i}² + [(ΔT_TSZ*B)ΔF]². On the negative side, the deconvolution coefficient does not scale linearly with redshift, and F_lin underestimates the true deconvolution factor F especially at high redshifts, potentially biasing our estimation of α. The effect would be more noticeable at 44 GHz, since the linear factor would underestimate the correction at small redshifts and overestimate it at higher redshifts (see the discussion of Figure 8(a) below).

In Figure 7 we present the likelihood function of a single simulation randomly selected of our ensemble of 1000 simulations. In Figure 7(a) we represent the likelihood for clusters within the three redshift bins given above, marginalized over cluster mass, and nine frequency ratios. Dashed, dot-dashed, and solid lines correspond to the lower, intermediate, and high redshift bins. For the latter, α = −0.052  ±  0.011. In Figure 7(b) we present the results binning clusters according to mass. Dashed, dot-dashed, and solid correspond to lower, intermediate, and high mass bins. As expected, the most massive clusters dominate the likelihood. The value α in this single simulation measured from the bin containing the most massive clusters is α = −0.028  ±  0.013. The full likelihood, including all redshifts and all masses, is represented in Figure 7(c) (dashed line). We also plot the distribution of the values of α measured in the thousand simulations (histogram). The thick solid line is a Gaussian fit to the histogram. For the single simulation, we obtain α = −0.048  ±  0.010 while the average over simulations is $\bar{\alpha }=-0.045,\; \sigma _{\bar{\alpha }}=0.010$. By keeping fixed the template of CMB residuals we verified that the effect of CMB residuals is to bias our estimate of α since the rms dispersion in α equals the error in α in one single realization. In this particular case, α is biased to negative values, but for different realizations of the residuals, the bias could be positive with equal probability. If the CMB residual is as large as ΔT_CMBres = 10 μK, then the bias rises to $\bar{\alpha }=0.6$. In the ideal case that the CMB template removes the cosmological signal exactly and no CMB residuals are present in the data, the result would be $\bar{\alpha }=-0.003\pm 0.011$; the error bar is similar but the bias becomes negligible.

The ratio method is biased both by the error in the denominator (Luzzi et al. 2009) and also by CMB and foreground residuals present in the data. For this method to provide an unbiased estimate of α, the residual contribution due to CMB and astrophysical foregrounds must be smaller than ∼1 μK. The bias can be estimated by analyzing clusters selected by galactic latitude, since, as indicated, CMB and foreground residuals are largest close to the galactic plane (see Planck HFI Core Team 2011b). Also, we checked that ratios using lower resolution maps have less statistical power to constrain α. In retrospect, this justifies neglecting the 30, 545, and 857 GHz channels, with large noise levels or (1/f) contributions, that would complicate the analysis without adding more information.

The results of the fit method using the 217 GHz map to remove the intrinsic CMB are presented in Figure 8. In panel (a) we plot the likelihood function of a single simulation for three different frequencies: 44 GHz (dashed), 100 GHz (solid), and 343 GHz (dot-dashed line). The figure shows that the 100 GHz channel is the most restrictive and the 44 GHz channel is the less restrictive of the three. The latter is also biased since the deconvolution factor is the largest of all the Planck channels (see Figure 5) and is not as well represented by a linear approximation with a small correction, especially at low masses (see Figure 6(a)). The biased introduced does not alter significatively the final results since the final likelihood is also dominated by the channels that have the highest resolution and lowest noise, in this case the 100 and 143 GHz channels. In Figure 8(b) we represent the likelihood for the three mass bins given above and marginalized over frequencies; dashed, dot-dashed, and solid lines correspond to the low, intermediate, and high mass intervals. The signal is dominated by the most massive clusters that, on a flux-limited sample, are, on average, at high redshift than the lower and intermediate mass samples. For this particular realization, the value estimated from the subsample of most massive clusters is α = 0.020  ±  0.018, adiabatic evolution being only marginally outside the 1σ level.

In Figure 8(c) we represent the histograms of 1000 simulations together with their linear fits for cluster templates constructed using the universal pressure profile (dot-dashed line) and the β = 2/3 (solid line). The mean and rms dispersion of the estimated values are 〈α〉 = −0.011  ±  0.018 for the universal profile and 〈α〉 = 0.009  ±  0.008 for the β-model profile. When all cluster properties are identical, the TSZ integrated over the cluster extent will be larger for the β-model than for the universal profile (see Figure 3(a)) so it must constrain α better, as shown. We also carried out 1000 simulations using method (A) with the hydrodynamical template, which contains the KSZ component and CMB residuals of 1 μK amplitude. The result was 〈α〉 = −0.008  ±  0.015, identical to the result with method (B) above. Therefore, in the fit method is not so important to have maps free of residuals as, for example, in the ratio method. Using the 217 GHz map to remove the intrinsic CMB signal alters the TSZ frequency dependence, but the TSZ signal is still strongly dependent with redshift (see Figures 1(b) and 2(b)) what is not the case in the ratio method (compare Figures 1(a) and 2(a)) and an extra error of ∼1 μK, even if correlated across frequencies, is not so significant.

Let us remark that the rms dispersion of α on 1000 simulations, σ_α, is very similar to the error on α in one single realization, both in the ratio and in the fit method, indicating that our pipelines are efficient. The results obtained using y-maps constructed with clusters drawn from a hydrodynamical simulation or from a catalog of X-ray-selected clusters show no significant differences. In the hydro-simulation, the y-map integrates the SZ signal up to z  ≃  0.25 and contains all the projection effects up to that redshift, not included in the X-ray-selected cluster template; then, we can conclude that projection effects play no significant role. This can be understood in the light of the results presented in Figure 8(b); the full likelihood is dominated by the most massive clusters for which projection effects are not significant (see A. Valente et al. 2012, in preparation for a discussion on this point).