THE IMPACT OF CLUSTER STRUCTURE AND DYNAMICAL STATE ON SCATTER IN THE SUNYAEV–ZEL'DOVICH FLUX–MASS RELATION

The simulation described here was performed using FLASH, an Eulerian hydrodynamics plus N-body code which has been applied to a wide range of problems and extensively validated for hydrodynamical (Calder et al. 2002) and cosmological N-body (Heitmann et al. 2005, 2008) applications. We used version 2.4 of FLASH together with the local transform-based multigrid Poisson solver described by Ricker (2008). Because we are concerned in this paper only with the effect of gravity-driven processes on mass-observable relations, the calculation described here did not employ radiative cooling or feedback due to star formation or active galaxies. Here, we give a brief summary and refer the readers to Yang et al. (2009) for details of the numerical methods and merger tree analysis.

The results presented here are based on a FLASH simulation of structure formation in the ΛCDM cosmology within a three-dimensional cubical volume spanning 256h⁻¹ Mpc. Initial conditions were generated for a starting redshift z = 66 using GRAFIC (Bertschinger 2001) with an initial power spectrum generated using CMBFAST (Seljak & Zaldarriaga 1996). The cosmological parameter values used were chosen to be consistent with the third-year WMAP results (Spergel et al. 2007): present-day matter density parameter Ω_m0 = 0.262, present-day baryonic density parameter Ω_b0 = 0.0437, present-day cosmological constant density parameter Ω_Λ0 = 0.738, matter power spectrum normalization σ₈ = 0.74, and Hubble constant h = 0.708 (H₀ = 100h km s⁻¹ Mpc⁻¹). The simulation contains 1024³ dark matter particles with a particle mass m_p = 9.2 × 10⁸h⁻¹ M_☉. The mesh used for the gas dynamics and potential solution was fully refined to 1024³ zones, which corresponds to a comoving zone spacing of 250h⁻¹ kpc. Considering the effect of resolution on the computed abundances of halos of different mass (Lukić et al. 2007), with these parameters we are able to capture all halos containing more than 3150 particles (i.e., total mass 2.9 × 10¹² h⁻¹ M_☉) and 1150 particles (i.e., 1.1 × 10¹²h⁻¹ M_☉) at z = 0 and 1, respectively. The halos are identified using the friends-of-friends (FOF) algorithm. The overdensity mass and radius, M_Δ and R_Δ, are then found by growing spheres around each FOF center until the averaged total density is Δ times the critical density of the universe. The simulation was carried out using 800 processors of the Cray XT4 system at Oak Ridge National Laboratory, requiring a total of 16,500 CPU hours.

The thermal Sunyaev–Zel'dovich effect is caused by inverse Compton scattering of CMB photons off the hot electrons inside galaxy clusters. Assuming that relativistic corrections are small, the resulting distortion of the CMB temperature can be written as ΔT/T_CMB = g_ν(x)y, where $g_\nu (x)=x(\coth (x/2)-4)$ with x = hν/k_BT_CMB, k_B is the Boltzmann constant, ν is the frequency of observation, and T_CMB is the mean CMB temperature at the current epoch. Here, y is the Compton y-parameter, which can be written as

$$\begin{equation} y= \frac{k_B\sigma _T}{m_ec^2}\int n_e(l)T_e(l)dl, \end{equation} \tag{ 1 }$$

where n_e(l) is the electron density profile, σ_T is the Thomson scattering cross section, and the integration is along the line of sight, which is defined to be the x-direction in the simulation box.

Note, however, that the observable is the integral of the temperature distortion over the cluster's projection onto the sky. For a cluster at redshift z, it is given by

$$\begin{equation} Y(M,z) = \frac{1}{d_A(z)^2} \frac{k_B\sigma _T}{m_ec^2}\int n_e(l)T_e(l)dV, \end{equation} \tag{ 2 }$$

where d_A(z) is the angular diameter distance to the cluster, and the integration is over the volume of the cluster. In the literature, sometimes the factor 1/d²_A is omitted from Equation (2). In this study, when we omit the factor 1/d²_A, we denote the temperature distortion as Y(M); otherwise it is denoted as Y(M,z). Note that Y(M,z) is dimensionless, while the units of Y(M) are Mpc². In the following sections in which we investigate the distribution and origin of intrinsic scatter, we adopt Y(M) unless explicitly stated otherwise.

Following the convention in the literature, we denote Y_Δ as the SZ flux integrated out to certain overdensity radius R_Δ. In the simulation box, cell-averaged information about the gas density and temperature is stored for each grid cell. Thus, for each cluster the integrated SZ distortion is calculated using

$$\begin{equation} Y_\Delta (M,z)= \frac{1}{d_A(z)^2} \frac{k_B\sigma _T}{m_ec^2}\sum _{i, r_p\le R_\Delta } n_{e,i}T_{e,i}\Delta V_{i}, \end{equation} \tag{ 3 }$$

where the summation is over all the grid cells across the cluster volume within a cylinder of projected radius R_Δ, and ΔV_i is the volume of each grid cell.

In order to explore the influence of cluster formation and merger history on the intrinsic scatter, we construct a merger tree for each cluster in our simulation in the following way. Our simulation generates snapshots that contain particle tags and positions every 100 h⁻¹ Myr beginning at z = 2. For each snapshot, all the groups with more than 10 particles are found using the FOF halo finder with linking length parameter b = 0.2. Between successive outputs at times t = t_n, we find the progenitors at time t_n−1 for all the halos at t_n by tracing the particle tags, which are uniquely assigned to each particle at the beginning of the simulation. For each halo at t_n, we record the masses of its progenitors, the masses they contributed to the halo, and the number of unbound particles. Then, the merger trees are constructed by linking all the progenitors identified in the previous outputs for halos above our halo completeness limit at z = 0. Deriving the mass accretion histories is straightforwardly accomplished by following the mass of the most massive progenitor back in time. Cluster formation time is often defined as the epoch when a cluster exceeds a certain fraction of its final mass. The commonly adopted thresholds include 10%, 25%, 50%, and 70%. In discussion below we present results using the 50% threshold.

To directly quantify the dynamical state of clusters without relying on morphology, we find the time since last merger for each cluster in our simulation. In our analysis, mergers are defined in two ways: the mass-jump definition, in which a merger is present if there is a mass jump larger than some threshold in the halo's assembly history, and the mass-ratio definition, which identifies a merger if the ratio of contributed masses from the first- and second-ranked progenitors is less than a certain value (Cohn & White 2005). To study the variations in cluster observables induced by different types of mergers, we use 1.2 and 1.33 as thresholds for the mass-jump definition and 10:1, 5:1, and 3:1 in the mass-ratio definition. In this paper, by "merging clusters" at a given lookback time we will refer to those identified by at least one of these five merger diagnostics in the preceding 3 Gyr, chosen to be long enough such that mergers with different impact parameters and mass ratios would have returned to virial equilibrium within R₅₀₀ (Poole et al. 2006). The mergers are "major" if the mass jump is larger than 1.2 or if the mass ratio is less than 5:1; "minor" mergers, on the other hand, have mass ratios between 10:1 and 5:1.

One of our main goals is to investigate the possible sources of intrinsic scatter, including the variations due to concentration, departure from hydrostatic equilibrium, merger boosts, and cluster morphology. In order to distinguish the contribution from each source, we construct a set of idealized cluster samples with different assumptions. Starting from a sample with the most possible constraints, we add one source of scatter at a time. By comparing the scatter of the idealized samples and the simulated sample, we can tell whether the scatter can be successfully reconstructed, or yet other sources still need to be found.

To this end, we construct five idealized samples, going from sample A, with the most constraints, to sample E, which includes the most sources of variation. For sample A, only the cluster masses are taken from the simulation. Given the mass, the halo concentration c is computed using the best-fit c–M relation from Shaw et al. (2006). With the mass and halo concentration, the total density and gas density are assigned using the NFW profile (Navarro, Frenk, & White 1995, 1996, hereafter NFW) and a core-softened NFW profile (Sijacki et al. 2007), respectively. In the core-softened NFW profile, the core radius is set to be 0.02 times the virial radius, and the gas fraction, f_g = 0.12, is also fixed. The pressure and temperature profiles are then computed assuming hydrostatic equilibrium (HSE). Finally, dark matter particles and zone-averaged cell quantities are assigned according to the profiles, assuming spherical symmetry. They are stored using the same file format as the simulated clusters, allowing them to be analyzed in the same way as the simulated clusters. Any scatter in the resulting Y–M relation can only be due to particle shot noise, finite cell resolution, and the simulated observation procedure.

Sample B is generated using a similar procedure, except that the c–M relation is assumed to have a lognormal distribution with a dispersion of 0.22 (Jing 2000; Bullock et al. 2001; Dolag et al. 2004; Shaw et al. 2006). In other words, the variation in concentration should be the only additional source of Y–M scatter for sample B.

In addition to the variation in concentration, clusters in sample C are allowed to depart from HSE due to pressure support from random gas motions after mergers. For each cluster, we compute the gas velocity dispersion's radial profile from the simulation and include an extra pressure term, P_rand = ρσ², in the HSE equation for computing the thermal pressure. In this way, we are effectively taking into account the incomplete virialization after merger events as another source of scatter.

With sample D, we model the effect of mergers by including not only incomplete relaxation but also the merger boosts due to shock heating (Ricker & Sarazin 2001; Poole et al. 2007). Therefore, we extract the time histories of mass and integrated SZ flux during mergers from the ideal merger simulations in Poole et al. (2007) and apply a boost in Y to each idealized cluster using the actual times since last merger in the simulation and the mass ratios of those mergers. Note that the time evolution of mass from Poole et al. (2007) is measured for overdensity Δ = 500, but the integrated SZ flux is only available for Δ = 2500. Thus the Y–M scatter of sample D should be considered as an upper limit to the effects of merger boosts.

Sample E adds variation in the gas fraction by using the actual gas fraction computed for each simulated cluster. Sample E essentially includes all possible sources of scatter investigated in this paper except the effect of cluster morphology. In Section 4.5, we will compare results from these idealized samples to the simulated clusters (sample S). Because the idealized samples are all constructed under the assumption of spherical symmetry, we derived another sample (sample SS) using gas profiles directly extracted from the simulated clusters, such that it includes all sources of scatter except the morphological effect. These models are summarized in Table 1.

For each simulation output between z = 0 and 1.5, we derive the Y_Δ − M_Δ relation for all clusters with M₅₀₀ ⩾ 2 × 10¹³ M_☉ (619 clusters at z = 0 and 223 clusters at z = 1) and fit it with a power law of the form

$$\begin{equation} \bar{Y}_\Delta =10^{-6} A_{14}(z) \left(\frac{M_\Delta }{10^{14} h^{-1} \,M_\odot } \right)^\alpha, \end{equation} \tag{ 4 }$$

where the normalization at 10¹⁴h⁻¹ M_☉ in units of 10⁻⁶, A₁₄, and slope α are found by fitting the data points using a Levenberg–Marquardt algorithm. An overdensity of Δ = 500 is adopted throughout the paper.

The normalization (relative to the value at z = 0), slope, and scatter as functions of redshift are plotted in Figure 1. As also found in previous adiabatic simulations (da Silva et al. 2004; Motl et al. 2005), the evolution of the normalization and slope is consistent with the self-similar prediction, α = 5/3 and A₁₄(z) ∝ E(z)^2/3, where E(z) = [Ω_m0(1 + z)³ + Ω_Λ0]^1/2, though we cannot rule out the case of no evolution because the small number of high-mass clusters at higher redshift limits the constraining power of the data. We find A₁₄(z = 0) = 5.43 ± 0.47, which is in agreement with the adiabatic runs in Nagai (2006) (A₁₄(z = 0) = 4.99, f_b = 0.14) and Shaw et al. (2008) (A₁₄(z = 0) = 4.07, f_b = 0.11) after taking into account the differences in the baryon fraction used in the simulations (f_b = 0.167 in our simulation).

Zoom In Zoom Out Reset image size

The rms scatter around the best-fit relation is defined for N clusters as

$$\begin{equation} \sigma _{YM} = \left[ \frac{\Sigma _{i=1}^N (\log Y_i-\log \bar{Y}_i)^2}{N-1} \right] ^{1/2}, \end{equation} \tag{ 5 }$$

where Y_i is the measured flux of the ith cluster, and $\bar{Y_i}$ is the flux predicted by the best-fit relation for that cluster. Hereafter, we use the notation

$$\begin{equation} \delta \log Y \equiv \log Y - \log \bar{Y} \end{equation} \tag{ 6 }$$

for the deviation from the mean relation for each cluster. For each redshift from z = 0 to 1, we compute the rms scatter for 5 mass bins and plot the result in Figure 2 (the data for some of the redshifts are omitted for clarity). The scatter is ∼5%–15%, consistent with previous findings (e.g., Nagai 2006). Moreover, we find that in general the scatter decreases with both mass and redshift. We fit the scatter using the functional form

$$\begin{equation} \sigma (M,z) = A \log M + B \log (1+z) + C, \end{equation} \tag{ 7 }$$

where the best-fit coefficients are A = −7.06 ± 0.28, B = −11.20 ± 0.81, and C = 7.70 ± 0.19. The mass dependence may be due to increasing non-lognormality of the scatter when considering low-mass clusters (see Section 3.3 for details). The redshift evolution may be understood by considering the self-similar model, in which all quantities for collapsed objects can be expressed in terms of the characteristic mass scale, M_⋆ ∝ (1 + z)^−6/(n+3), where n is the spectral index of the scale-free primordial power spectrum, P(k) ∝ kⁿ (Kaiser 1986). Therefore, Equation (7) is equivalent to the expression σ = A'log(M/M_⋆) + B'. Note that Equation (7) is the first attempt in the literature to quantify the scatter using a functional form of mass and redshift. This expression should be useful for future studies that require assumptions about the form of scatter.

Zoom In Zoom Out Reset image size

Since the SZ flux is only proportional to the first power of gas density, it is sensitive to the contribution from low-density gas clumped along the line of sight to but not associated with a cluster. This is in contrast to the X-ray luminosity, which is proportional to the square of gas density. For this reason, the Y–M relation has been found to have a high-scatter tail in the distribution of its scatter (White et al. 2002; Hallman et al. 2007). Since any deviations from the lognormal scatter would bias cosmological constraints based on cluster counts (Shaw et al. 2010a), we would like to examine whether the form of the log scatter can be well approximated by a Gaussian distribution, or whether generalizations of the parameterization need to be considered.

A purely Gaussian distribution can be described exactly using only its mean value μ and variance σ²:

$$\begin{equation} G(x) = \frac{1}{\sigma \sqrt{2\pi }}\exp \left[-\frac{(x-\mu)^2}{2\sigma ^2}\right]. \end{equation} \tag{ 8 }$$

A nearly-Gaussian distribution can be approximated using the Edgeworth expansion (e.g., Bernardeau & Kofman 1995; Blinnikov & Moessner 1998),

$$\begin{equation} \tilde{G}(x) \approx G(x) - \frac{\gamma }{6} \frac{d^3 G}{dx^3} + \frac{\kappa }{24} \frac{d^4 G}{dx^4} + \frac{\gamma ^2}{72} \frac{d^6 G}{dx^6}, \end{equation} \tag{ 9 }$$

which is parameterized by four moments—the mean and the variance describing the Gaussian distribution, plus the skewness (γ) and the kurtosis (κ) describing the deviation from Gaussianity. The skewness is defined as

$$\begin{equation} \gamma = \frac{\langle (x-\mu)^3\rangle }{\sigma ^3} \end{equation} \tag{ 10 }$$

and the kurtosis as

$$\begin{equation} \kappa = \frac{\langle (x-\mu)^4\rangle }{\sigma ^4}-3. \end{equation} \tag{ 11 }$$

We compute the skewness and kurtosis of the Y–M scatter for clusters at z = 0 above different mass thresholds and plot the results in Figure 3. The error bars represent the uncertainty due to finite sample size and are estimated using 10³ Monte Carlo realizations of random sampling from a nearly-Gaussian distribution given the measured skewness and kurtosis as in Equation (9). Note that because of finite sample size the measured skewness and kurtosis would underestimate the intrinsic values of the underlying distribution. This bias is represented by the offset between the data points and the middle points of the error bars.

Zoom In Zoom Out Reset image size

We find that the scatter is non-lognormal with positive skewness and kurtosis when including only massive clusters (M_500,lim ≳ 10¹⁴ M_☉) or when more and more low-mass clusters are included (M_500,lim ≲ 5 × 10¹³ M_☉). Due to the limited number of clusters at the high-mass end, the log scatter there is expected to follow Poisson statistics and deviate from a Gaussian form. For the lower mass range, on the other hand, we find that there is a tail toward positive values in the distribution of scatter, which increases the skewness and kurtosis. By visual inspection of clusters in the tail, we find that these objects happen to have elongated shapes or clumped gas along the line of sight. Since less massive clusters are more likely to be surrounded by gas with mass comparable to their own, this effect becomes more important for lower mass clusters. We will further address this point in Section 4.3. For other redshifts, the level of non-lognormality is also non-negligible, as shown in Figure 4. The median skewness and kurtosis are 1.43 and 4.21, respectively.

Zoom In Zoom Out Reset image size

The concentration of a dark matter halo is usually defined as the ratio between the virial radius and the NFW scale radius, i.e., c ≡ R_vir/R_s. The concentration parameter characterizes the density inside the core region of a halo and reflects the mean density of the universe when the halo collapsed. Thus, halos formed earlier in time tend to be more concentrated (NFW 1997; Wechsler et al. 2002). By testing the correlation of Y–M scatter with scatter in the concentration parameter, we are effectively probing the influence of cluster formation history.

We choose to use the parameter R₂₀₀/R₅₀₀ instead of the original halo concentration parameter c because it has two advantages. The first is that it avoids introducing the uncertainty of fitting an NFW profile, especially for less massive clusters, since the fitting is very sensitive to the grid resolution in the central region of a cluster. Moreover, our analyses involve not only relaxed clusters but also merging ones, for which R₂₀₀/R₅₀₀ is actually better defined than c, since an NFW profile would yield a poor fit.

Figure 5 (left panel) shows a strong positive correlation between scatter in the Y–M and (R₂₀₀/R₅₀₀)–M₅₀₀ relations. The correlation coefficient is 0.64, with a probability of zero given by the Spearman Rank-Order Correlation test (Press et al. 1992, Section 14.6; probability of one means no correlation). To ensure that this result is not biased by the lower mass clusters, whose R₅₀₀ values are close to the resolution of the simulation, we raised the mass threshold to M₅₀₀ ⩾ 10¹⁴ M_☉ and found that the result is robust for these well-resolved systems (shown as the open squares in Figure 5). Note that we correlate with δlog(R₂₀₀/R₅₀₀) instead of the raw value of R₂₀₀/R₅₀₀ because the latter is a function of cluster mass. By doing so we exclude the effect of different cluster masses, focusing on the variation in halo concentrations. R₂₀₀/R₅₀₀ is a monotonically decreasing function of the halo concentration parameter (see Yang et al. 2009 for derivation). Therefore, for clusters with similar masses, more concentrated clusters tend to lie under the mean Y–M relation, while the "puffier" clusters tend to scatter high.

Zoom In Zoom Out Reset image size

Since halo concentration is related to cluster formation time, we can test the above trend by checking the correlation of Y–M scatter with the formation times derived from the mass assembly histories of our simulated clusters. The formation time here is defined as the time when the cluster first exceeds half of its final mass. As expected, we find that clusters that formed earlier (thus with higher concentrations) tend to scatter low (see the right panel of Figure 5). The correlation is not as tight as the one with the halo concentrations. This is due to the fact that the correlation between the halo concentration and the cluster formation time itself has a very large scatter, and also that the variation in halo concentrations cannot be fully accounted for by the variation in cluster formation time (Neto et al. 2007). But the direction of the correlation with cluster formation time is consistent with the correlation with halo concentration.

To explain the correlation between the Y–M scatter and the concentration, recall the virial theorem for the simplest case of an isolated system: 2K + U = 0, where K and U are the total kinetic and gravitational binding energies of the system, respectively. In general, one can write

$$\begin{equation} \frac{k_B T}{\mu m_p} \propto \frac{{\it GM}}{R}, \end{equation} \tag{ 12 }$$

where μ is the mean molecular mass of the gas, m_p is the mass of a proton, and T, M, and R are the virial temperature, mass, and radius of the system, respectively. Together with the definitions of the SZ flux (Equation (2)) and $M={4\over 3}\pi R^3 \bar{\rho }$, one can derive

$$\begin{equation} Y \propto M_{\rm{gas}} T \propto f_g M^{5/3}. \end{equation} \tag{ 13 }$$

Note that the above relations are for virial quantities of a cluster as a whole, but mass-observable relations are often measured using a certain aperture size, R_Δ. The relation between the virial and overdensity quantities depends on individual cluster profiles, which are determined by the halo concentration and how the gas is distributed on top of the dark matter potential (e.g., the equation of state of gas). Therefore, the normalization, and thus the scatter, of the Y_Δ–M_Δ relation is a function of halo concentration and gas properties.

It is also important to note that, from the above derivation, the direction of the correlation between the scatter and concentration is dependent on the gas properties. Our result shows that less concentrated clusters tend to scatter high, i.e., have higher pressure than clusters of similar masses. Observationally, Comerford et al. (2010) have also found a similar anti-correlation between the X-ray temperature–mass scatter and strong lensing concentration. This may be attributed to the fact that less concentrated clusters have larger scale radii, and hence when comparing with clusters of the same M_Δ or the same aperture size R_Δ, their ratios R_Δ/R_s are smaller, which means the observable integrated within R_Δ would be greater. However, different simulations can have different directions of correlation depending on the input gas physics. For example, Shaw et al. (2008) also found a correlation between the scatter and concentration, but in the opposite direction. The difference may be due to different gas physics included in their models. In principle, by assuming a particular gas model the constant of proportionality in Equation (13) can be computed exactly. Ascasibar et al. (2006) have done this exercise assuming a polytropic equation of state for the gas. According to their calculation, the coefficient in the T–M relation (as in Equation (12); inverse of the Y_MT in their Equation (17)) decreases with concentration for a polytropic index of γ_p = 5/3. As γ_p decreases, the dependence becomes weaker and then the direction is reversed. Since including extra baryonic physics effectively works to decrease γ_p (e.g., for fixed mass and concentration, both changes yield a shallower temperature profile; see also Figures 2 and 3 in Ostriker et al. 2005), this may explain why the dependence on concentration can be different between models with different input gas physics. Note, however, that in reality the situation can be even more complicated because a constant γ_p may not be valid for all gas models (Kay et al. 2004).

The strong correlation in Figure 5 suggests that the variation in halo concentrations contributes a significant amount of the intrinsic scatter in the Y–M relation. Using this strong correlation it is possible to adjust for the dependence of SZ flux on cluster concentrations. We use δlog(R₂₀₀/R₅₀₀) for each cluster to calculate its expected δlog Y from the best-fit relation, (δlog Y)_exp = 7.167 × δlog(R₂₀₀/R₅₀₀). We then subtract this (δlog Y)_exp from the measured SZ flux to obtain a corrected flux,

$$\begin{equation} (\delta \log Y)_{\rm corr} = \delta \log Y - (\delta \log Y)_{\rm exp}. \end{equation} \tag{ 14 }$$

The Y–M relations before and after correcting for concentration are shown in Figure 6. We find that after removing the effect of halo concentration, the rms scatter decreases from 12.07% to 7.34% (i.e. by 38.9%). This method was proposed by Yang et al. (2009) to tighten the X-ray temperature–mass relation and has been successfully applied to observed strong lensing clusters (Comerford et al. 2010). In addition to strong lensing, the NFW concentration can also be measured via weak lensing, X-ray emission, etc. (Comerford & Natarajan 2007; Mandelbaum et al. 2008; Buote et al. 2007, and references therein), although one has to be careful about systematics of each method and when combining different measurements. In fact it can be generalized to any observable other than concentration. That is, if there exists any variable X whose mass scatter δlog X is known and correlates with δlog Y, then given the best-fit correlation (δlog Y)_exp = α × δlog X, the Y–M scatter can be reduced in a similar way using Equation (14) to remove the effect of X from the scatter. Therefore, this method can be a powerful way to reduce the observed mass-observable scatter and obtain better mass estimates.

Zoom In Zoom Out Reset image size

Another possible origin of the Y–M scatter is a cluster dynamical state. Cluster mergers are among the most energetic events in the universe. Shock heating and departure from hydrostatic equilibrium during mergers can drive clusters away from the mean scaling relations. Ideal merger simulations (Ricker & Sarazin 2001; Poole et al. 2007) have shown that the effect of shock heating is to boost the SZ and X-ray observables to values a few times higher than the pre-merger values. Studies that combine the amount of boosting predicted by these simulations with extended Press–Schechter merger trees (Randall et al. 2002; Wik et al. 2008) show that the boosting effect can bias estimates of cosmological parameters such as σ₈ and Ω_m. The other effect of mergers is departure from hydrostatic equilibrium. Before the gas within a merger is completely virialized, the pressure support from random gas motions can contribute ∼10%–20% of its thermal pressure (Rasia et al. 2006; Lau et al. 2009). Thus, the thermal pressure and hence the SZ flux of unrelaxed systems is expected to be smaller than that of relaxed systems of similar masses. These previous studies are primarily based on small cluster samples. Therefore, our aim is to investigate how merger events statistically influence the cluster scaling relations.

If the scatter were dominated by the boosting effect of cluster mergers as described above, then one would expect to find merging clusters to preferentially lie above the mean relation. However, if during mergers the departure from hydrostatic equilibrium due to non-thermal pressure support were dominant, mergers would tend to scatter low. In order to see which effect is more prominent, we correlate the Y–M scatter with the time since last merger (Figure 7, left panel). Substructure measures such as centroid offset (Mohr et al. 1995) and power ratios (Buote & Tsai 1995, 1996) are often used to quantify departures from equilibrium in clusters. We compute the centroid offset and power ratios for the simulated clusters using the same definition as in Yang et al. (2009). The right panel of Figure 7 shows the correlation with one of the power ratios, P₂/P₀. Based on the Spearman Rank-Order Correlation test (Press et al. 1992, Section 14.6), both correlations are weak (with a small correlation coefficient) but significant (with a high probability), in the direction that more disturbed clusters tend to scatter low. This implies that during mergers the incomplete virialization may be the more important factor in driving the scatter than the shock heating effect. However, the fact that these two effects operate in opposite directions may be the reason why there is not a clear trend with merger activities. Moreover, a number of factors can dilute the shock boosting effect, such as the small chance of finding mergers in progress, capturing the shocks within an overdensity radius at the right projection and at the right moment during a merger's transient boost, and the fact that merging clusters tend to move along the scaling relations because their masses also increase at the same time their observables increase, as also found by Wik et al. (2008) and Kravtsov et al. (2006); see Yang et al. (2009) for an extensive discussion.

Zoom In Zoom Out Reset image size

Do mergers bias the Y–M relation due to incomplete virialization? To answer this question, we plot the normalized distributions of the Y–M scatter for relaxed and merging clusters in Figure 8 and use the Wilcoxon Rank-Sum (R-S) test and the F-variance (F-V) test to see whether these two populations have significantly different mean values or variances, respectively. A value smaller than 0.05 (for a significance level of 5%) returned by these tests is commonly adopted to indicate a significant difference between two populations. We find that their mean values do not differ significantly, but mergers have a wider distribution compared to relaxed clusters (with significance 0.014). Therefore, although mergers do not tend to bias the scaling relation, they do have a greater amount of scatter than relaxed clusters. If merging (relaxed) clusters are chosen to be those within the quartile with the highest (lowest) substructure measures, similar trends with significant probabilities are found for 9 out of 21 substructure measures (P₂/P₀, P₃/P₀, and the centroid offset W measured from different viewing directions and varying aperture sizes). We find that including only the merging clusters would result in ∼15%–45% greater scatter than when only relaxed clusters are taken into account, consistent with previous findings (Shaw et al. 2008). Note however that separating mergers from relaxed clusters does not reduce the skewness or kurtosis of the scatter distribution, which suggests that the non-lognormality has causes other than mergers (see Section 4.3).

Zoom In Zoom Out Reset image size

Despite the spherical symmetry that theoretical models usually assume, both simulations (Warren et al. 1992; Cole & Lacey 1996; Jing & Suto 2002; Bailin & Steinmetz 2005; Kasun & Evrard 2005) and observations (e.g., Basilakos et al. 2000) have shown that clusters are triaxial (or elliptical when projected) rather than simple spheres, even for relaxed clusters. How the gas is distributed in the cluster potential well should vary depending on the axes ratios of the cluster. Moreover, viewing a triaxial cluster from different angles should also yield different observed quantities integrated along the line of sight. Both these factors can contribute to the Y–M scatter.

In order to explore the impact of morphology, for each simulated cluster we find the orientation of the principal axes by diagonalizing the moment of inertia tensor, I_αβ = m_iΣ_ir^α_ir^β_i, where the summation is over all the particles and cells in the cluster, m_i is the mass of a particle or a gas cell, and r^α_i is the x, y, or z component of the distance from the cluster center of mass. The lengths of the major, intermediate, and minor axes, denoted as a, b, and c, are found by finding the intercepts of the axes with the isodensity surface having overdensity Δ = 200. The angles θ_α are defined to be the angles between the major axis and the directions of projection (α = x, y, z; note that the x-direction is the projection used for all the analyses in this paper).

Motivated by our results in Section 3.3 that the non-lognormality may be due to clusters that happen to have elongated shapes aligned with the viewing direction, we invented a measure, acos θ_x, to trace cluster morphology along the line of sight. Figure 9 shows the Y–M scatter versus the scatter in the acos θ_x–M relation. The positive correlation indicates that clusters that are more elongated along the line of sight preferentially have higher Y–M scatter, which is expected because the SZ flux is roughly proportional to the column density of the gas (see Equation (1)). Given the best-fit relation, δlog Y = 0.076 × δlog(acos θ_x), we can again reduce the scatter by applying a correction as in Equation (14). By doing so we find that the scatter is reduced from 12.07% to 11.63% (i.e., by 3.6%).

Zoom In Zoom Out Reset image size

We further divide the cluster sample in half using the values of acos θ_x and plot the normalized distributions of the Y–M scatter in Figure 10. We find that the scatter distribution of the clusters with larger acos θ_x is non-lognormal (γ = 1.04, κ = 2.63), while the skewness and kurtosis of the remaining population are greatly reduced (γ = 0.29, κ = −0.08). Therefore, the non-lognormality is indeed caused by the clusters with more elongated shapes along the line of sight.

Zoom In Zoom Out Reset image size

Since the SZ flux is obtained by integrating along the line of sight within a projected radius, it is subject to contamination by gas that lies along the same line of sight, which causes the large number of high-scatter objects in the Y–M relation found in simulations that include light cones (White et al. 2002; Hallman et al. 2007). These outliers and the outliers due to morphology discussed in the previous section can both drive the non-lognormality of the scatter. Since our simulated observations only include isolated clusters and thus do not take the projection effect into account, the skewness and kurtosis estimated from our simulation may be considered to be underestimates of the true values. In this case, it is even more important to adopt the higher moments in the parameterization of the scatter in order to get unbiased cosmological constraints.

Figure 11 shows the percentage contribution of scatter for each idealized cluster sample with respect to the simulated sample. The cluster sample at the bottom has the most constraints on and least freedom in the model parameters, and the assumptions are loosened one at a time from bottom to top (see Table 1 for a summary of model descriptions). In general, the values are independent of mass, that is, the processes shape the scaling relation in a self-similar way, as expected in the absence of additional baryonic physics. The only exception is the bottom curve for which only the masses of clusters are assigned. In principle, this sample should have zero scatter if the resolution of gas cells were infinite, but in reality the finite resolution introduces a nonzero scatter which becomes bigger as more low-mass clusters are included. Because the idealized samples are all constructed under the assumption of spherical symmetry, we derived sample SS (second line from the top in the figure) using gas profiles directly extracted from the simulated clusters, such that it includes all sources of scatter except the morphological effect. In other words, the difference between the simulated clusters and the spherically smoothed clusters is solely due to the variation in cluster morphology, which is ∼10%.

Zoom In Zoom Out Reset image size

From the differences between the subsequent samples, we are able to isolate the contribution of each effect to the total scatter: the variation in halo concentration contributes ∼10%–20% (difference between A and B), the departure from hydrostatic equilibrium results in ∼10%–15% (between B and C), merger boosts add another ∼30%–60% (between C and D), the variation in gas fractions introduces ∼0%–10% (between D and E), and the rest (between D and SS) due to other unaccounted for effects is ∼0%–30%. Note that the contribution from variation in concentration quoted here is estimated using sample B, which assumes spherical symmetry and hydrostatic equilibrium. However, in reality, both changing cluster morphology and including random gas motions (Lau et al. 2009) would further modify the distribution of gas and hence alter the measured value of concentration. In other words, the scatter is driven not only by the variation of concentration in sample B, but also by that due to cluster morphology and departure from HSE. Therefore, the total effect of concentration, as suggested by the correction for concentration in Section 4.1 (i.e., ∼40%), would be more appropriately accounted for also by considering the contributions from morphology and random gas motions.