A LINE-OF-SIGHT GALAXY CLUSTER COLLISION: SIMULATED X-RAY OBSERVATIONS

We performed our simulations using FLASH, a parallel hydrodynamics/N-body astrophysical simulation code developed at the Center for Astrophysical Thermonuclear Flashes at the University of Chicago (Fryxell et al. 2000). FLASH uses adaptive mesh refinement (AMR), a technique that places higher-resolution elements of the grid only where they are needed. In our case, we are interested in capturing sharp ICM features like shocks and "cold fronts" accurately, as well as resolving the inner cores of the cluster DM halos. As such it is particularly important to be able to resolve the grid adequately in these regions. AMR allows us to do so without needing to have the whole grid at the same resolution. FLASH solves the Euler equations of hydrodynamics using the Piecewise–Parabolic Method (PPM) of Colella & Woodward (1984), which is ideally suited for capturing shocks. FLASH also includes an N-body module which uses the particle-mesh method to solve for the forces on gravitating particles. The gravitational potential is computed using a multigrid solver included with FLASH (Ricker 2008).

The cluster DM halos are initialized using the Navarro-Frenk-White profile (NFW; Navarro et al. 1997) density profile, where we follow the method of Kazantzidis et al. (2006). The NFW functional form is used for r ⩽ r₂₀₀:

$$\begin{equation} \rho _{\rm DM}(r) = {\rho _{\rm s} \over {r/r_{\rm s}{(1+r/r_{\rm s})}^2}},\ r \le r_{200}. \end{equation} \tag{ 1 }$$

Outside of r₂₀₀, to keep the mass of the halo finite and to avoid an unphysical distribution function, we implement an exponential cutoff that turns off the profile on a scale r_decay, a free parameter which we set to 0.1r₂₀₀:

$$\begin{equation} \rho _{\rm DM}(r) = {\rho _{\rm s} \over {c_{200}(1+c_{200})^2}}{\left({r \over r_{200}}\right)^\kappa }{\rm exp}{\left(-{r-r_{200} \over r_{\rm decay}}\right)},\ r > r_{200}. \end{equation} \tag{ 2 }$$

Here r₂₀₀ is the radius within which the mean mass density is 200 times the cosmic mean, r_s is the NFW scale radius, and c₂₀₀ ≡ r₂₀₀/r_s is the NFW concentration parameter. We require that at r₂₀₀ the profile and its logarithmic slope be continuous. This is achieved via the parameter κ:

$$\begin{equation} \kappa = -{(1+3c_{200}) \over (1+c_{200})} + {r_{200} \over r_{\rm decay}}. \end{equation} \tag{ 3 }$$

Truncating the cluster NFW profiles at r₂₀₀ is unphysical. In real DM halos, the density profile is expected to remain NFW-like out to where the density of the halo reaches the average matter density of the universe (Prada et al. 2006; Tavio et al. 2008). The single cluster tests we have performed (as well as those performed in other studies, e.g., Ricker & Sarazin 2001) demonstrate that the truncated density profile is stable for the length of the run, with minor steepening of the profile near the cutoff radius.

The truncation also reduces the halo mass by approximately 25%. Assuming the clusters fell in from infinity, the larger halo masses would accelerate the halos to an impact velocity approximately 12% larger than the value we use. After the collision, the corresponding deceleration due to gravity would also be higher. The effect of the larger halo masses on the deceleration due to dynamical friction (f_dyn ∝ ρM2/v2) is more complicated to determine, since although the density in the outskirts of the halos is higher (thus increasing the dynamical friction), the velocity of clusters would also be higher due to the increased acceleration (thus decreasing it). The higher infall velocity of the clusters and the increased mass outside r₂₀₀ might also lead to some increase in the dissipation of the shock energy, and therefore in the heating of the gas in the outskirts. However, in general any such effects are inaccessible to current X-ray observations, due to the low densities in the cluster outskirts. We conclude that the effects due to truncation of the cluster density profile are not likely to be significant for the purposes of this study.

To initialize the gas density we choose a Burkert profile:

$$\begin{equation} \rho _{\rm gas}(r) = {\rho _{\rm c} \over {(1 + ({r/r_{\rm c}})^2)(1 + {r/r_{\rm c}})}}\,, \end{equation} \tag{ 4 }$$

which originally (Burkert 1995) was given as a fitting function for DM profiles of dwarf galaxies but is also a good fit to the gas density profiles of non-cooling-core clusters (A. Kravtsov 2007, private communication). The more traditional β-model (Cavaliere & Fusco-Femiano 1976) is a poorer fit to the gas density profile at large r in both simulations and observations (Vikhlinin et al. 2006; Hallman et al. 2007; Nagai et al. 2007). (However, we follow Ota et al. 2004 and Zhang et al. 2005 and fit β-model profiles to the X-ray surface brightness profiles in our simulated observations.)

The gas density is also fitted to an exponential taper at r₂₀₀, extending to the radius at which it equals the mean baryonic density of the universe. For the gas profiles, we choose the core radius r_c to be roughly half the scale radius of the DM profile (Ricker & Sarazin 2001). We fix the gas mass fraction at r₂₀₀ to 0.12 in line with observations of real clusters (e.g., Vikhlinin et al. 2006, Sanderson et al. 2003), and the measurement of Ota et al. (2004).

From these parameters the central gas densities ρ_c can be determined. The equation of state for the gas is an ideal-gas equation of state with γ = 5/3 and mean atomic mass μ = 0.592. We determine the pressure, temperature, and internal energy profiles by assuming our functional forms for DM and gas density and numerically integrating the equation of hydrostatic equilibrium. The values used for the cluster parameters are given in Table 1.

Following Czoske et al. (2002) and Jee et al. (2007) we assume a mass ratio of 2:1 for the clusters. In our simulation the clusters are initialized so that their centers are separated by the sum of their respective r₂₀₀ values (approximately 3 Mpc), and they are given an initial relative velocity v_rel = 3000 km s⁻¹. This value is chosen because it is the inferred relative velocity of the two redshift components of Cl 0024+17 seen in the redshift histograms of Czoske et al. (2002). The results of Czoske et al. (2002) and Jee et al. (2007) suggest the clusters have already completed their first pericentric passage, and Czoske et al. (2002) show that in their collision simulation an initial velocity of ∼3000 km s⁻¹ is needed to reproduce the observed distribution of redshifts after the collision. This value for the velocity is approximately 2v_ff for the clusters, where v_ff is the free-fall velocity from infinity. Such an initially high relative velocity might be difficult to achieve in a ΛCDM universe (see, e.g., Hayashi & White 2006). However, our choice is motivated by our desire to match the relative inferred velocity of the cluster components of Cl 0024+17.

We refine the adaptive mesh using the second derivative of density, temperature, and pressure to capture the sharp features in the gas dynamics, and using the matter density to resolve the cores of the clusters. For our box size of 14.29 Mpc we achieve a finest resolution of Δx = 13.96 kpc.

To make a meaningful comparison with the X-ray observations of Cl 0024+17, we construct mock X-ray observations of the simulation output. To do this we use MARX, a suite of programs created by the MIT/CXC group that simulates the on-orbit performance of the Chandra X-ray Observatory.³ MARX provides detailed ray-tracing simulations of astrophysical sources as observed by Chandra. In this paper we use MARX to generate simulated surface brightness profiles and spectra.

For input into MARX we create projected flux maps from the FLASH simulation data. For each FLASH zone the MEKAL model (Mewe et al. 1995) is used to assign a spectral emissivity. For these simulations zero metallicity has been assumed. The spectra are projected along the appropriate line of sight, and flux maps are generated for photons in 198 separate energy bins from 0.1–10.0 keV, giving an energy resolution of ΔE = 50 eV, compared to the ΔE≈ 100 eV resolution of Chandra.

For our simulated Chandra observations we use the simulated ACIS-S detector and HRMA mirror system, as in the real Chandra observation of Cl 0024+17. Each observation is ∼40 ks of simulated exposure time, comparable to the Chandra observation of Cl 0024+17 (39 ks). The flux maps serve as a distribution function for the simulated photons and determine the direction cosines, energies, and time steps for the photons in MARX. For simplicity, effects such as foreground galactic absorption and point source contamination are not included. However, we do include a diffuse X-ray background component. Output is given as a FITS file that can be read and analyzed in the same way as real Chandra observations. To analyze our simulated data, we use CIAO 3.4 and XSPEC 12.3. Details of our procedure for generating these observations and a verification study are given in the Appendix.

Figures 1 and 2 show slices of density and temperature, respectively, through the z = 0 coordinate plane at different times in the simulation (the x-axis is the collision axis). From t = 0.0–0.9 Gyr the clusters approach each other and a shock front forms. At t = 1.0 Gyr the core of the smaller cluster collides with the core of the larger cluster, driving a shock wave forward into the larger core and displacing the larger cores gas from the collision axis into streams of cold gas. There is also a stream of dense, cold gas that is pulled directly between the two DM halos. Later on, the heads of these streams of gas, along with hotter gas from between the clusters, begin to fall into the potential of the larger halo and by t = 2.0 Gyr have fallen in completely. As for the gas core of the smaller cluster, it is penetrated shortly after the collision by a reverse shock, but this shock weakens as it traverses the core and so the latter remains cool, in agreement with similar investigations of high-velocity cluster mergers by Takizawa (2005) and Milosavljević et al. (2007). In this process some of the gas of the smaller cluster is also ram-pressure stripped, and a contact discontinuity forms that is similar to that observed in the 1E0657-56, and in the simulations of the latter (Springel & Farrar 2007; Mastropietro & Burkert 2008).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After penetrating the larger gas core the smaller core is forced to lag behind its DM core by ram pressure until it encounters regions of lower density (t > 1.2 Gyr). At this point the gas falls back toward the DM core but overshoots the potential and begins to slosh back and forth inside it. At late times the gas of the clusters is still settling into their respective DM potential wells, and the temperature structure of the gas is still complicated. Due to the high initial velocity of the collision, the clusters are in an unbound orbit and by the end of the simulation the separation between their respective DM halo centers is still increasing. This is in contrast to many previous works involving mergers (e.g., Ricker & Sarazin 2001; Poole et al. 2006), which assumed an initial velocity v ≈ v_ff.

For the X-ray observations we view the system along the line of centers. Since we do not know the current state of the actual collision, we set each individual observation as if it were being observed at a redshift z = 0.395. At this redshift for the assumed cosmology 1'' = 6.39 kpc. Figure 3 shows an example of a smoothed mock X-ray image. For each observed time we extract a radial surface brightness profile and a radial temperature profile.

Zoom In Zoom Out Reset image size

Cluster surface brightness profiles are extracted for the energy range 0.5–5.0 keV, within a circle centered on the surface brightness peak with a radius of 400'' (≈ 2500 kpc). Figure 4 shows the surface brightness profiles for times t = 2.0, 2.5, and 3.0 Gyr. The profiles are fitted with the β-model (Cavaliere & Fusco-Femiano 1976):

$$\begin{equation} S_X(r) = S_0\left[1+\left({r \over r_c}\right)^2\right]^{-3\beta +1/2}. \end{equation} \tag{ 5 }$$

We fit with a single β-model and a sum of two β-models. The χ²-statistic is used to fit the profiles. We follow Ota et al. (2004) and Jee et al. (2007) in fixing one of the β parameters to unity, as there exists a strong degeneracy between this parameter and its corresponding core radius r_c. We include the background constant as a free parameter in the single β-model fit, but freeze this value for the double β-model fit. The fitted profiles for the t = 3.0 Gyr case are shown in Figure 5. The values of the fitted parameters for epochs t = 2.0, 2.5, and 3.0 Gyr are given in Tables 2 and 3. Errors are quoted at the 90% confidence level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Fitted Cluster Parameters for a Double β-model

t (Gyr)	$r_{\rm c_1}$ (kpc)	S0,1 (counts s−1cm−2arcsec−2)	β1	$r_{\rm c_2}$ (kpc)	S0,2 (counts s−1cm−2arcsec−2)	β2	χ2/dof
0.0	89.24+3.16−0.32	(3.73+0.04−0.04) × 10−6	1.0 (F)	204.65+1.32−4.00	(7.03+0.07−0.07) × 10−7	0.981+0.003−0.014	260.01/193
2.0	39.91+34.41−17.45	(8.18+3.41−3.39) × 10−8	1.0 (F)	270.00+6.34−4.78	(1.11+0.02−0.02) × 10−7	0.970+0.010−0.025	245.83/193
2.5	79.28+9.04−5.36	(2.26+0.17−0.17) × 10−7	1.0 (F)	282.00+11.97−5.57	(7.70+0.21−0.21) × 10−8	0.970+0.021−0.021	226.75/193
3.0	82.12+8.12−3.94	(2.51+0.16−0.16) × 10−7	1.0 (F)	316.24+26.16−4.89	(4.94+0.16−0.16) × 10−8	0.980+0.015−0.049	205.86/193

Download table as:  ASCIITypeset image

The projected temperature profiles for t = 2.0, 2.5, and 3.0 Gyr are shown in Figure 6. The profiles were generated by extracting spectra from four annular regions in the 0.5–7.0 keV band centered on the peak of the surface brightness profile. Each is fitted with a MEKAL model, using the χ²-statistic with grouping to ensure at least 15 photons per bin. The radial extent of the temperature profiles is 125'' (≈ 800 kpc), which is slightly greater than the radial range for which the temperature profile was extracted in Ota et al. (2004). Single-temperature fits to the spectrum of the whole cluster in the 0.5–7.0 keV band for t = 2.0, 2.5, and 3.0 Gyr are given in Table 4.

The feature in the radial profile of Cl 0024+17 pointed out by Ota et al. (2004) and Jee et al. (2007) was that a double β-model fit is a better fit to the data, as a single β-model fit underestimates the central surface brightness. The surface brightness profile of Cl 0024+17 and the single and double β-model fits to the profile are shown in Figure 7. The difference between the fits is significant; there is a difference in the χ²-statistic of ∼60 for a difference of 2 degrees of freedom.

In our mock cluster observations, we find that at later times, it is also possible to make a distinction between the two surface brightness components of the two clusters from the difference in the fitted χ²-statistic (see Tables 2 and 3). At the t = 2.0 Gyr epoch, Δχ² ≈ 45 for a difference of only 2 degrees of freedom, and at later times is even larger (Δχ²≈ 200–260 for a difference of 2 degrees of freedom). We also note that we can use the Δχ²-test to distinguish between the two clusters in our initial conditions, which we also show in Tables 2 and 3.

In Cl 0024+17, the fitted spectral temperature of T = 4.47 keV seems low for a system undergoing a collision or merger. In our mock observations we note a similar phenomenon; the fitted temperatures in the mock observations are significantly lower than what might be expected given the temperatures observed in the simulation. In our simulated system, the hottest cluster gas is in the larger cluster and is at a temperature of ∼5–6 keV at later times. However, the measured temperatures in the mock X-ray observations barely reach ∼3 keV. A close look at the temperature distribution in the simulation itself (Figure 8) reveals that the highest temperature gas is confined to the central region of the larger cluster, and the denser gas of the smaller cluster is colder, around 3 keV. In projection along the line of sight, the lower-temperature, denser gas "washes out" the higher-temperature, less dense gas, and the complicated, two-component temperature structure of the clusters is lost because of the strong density dependence of the X-ray emission (∝ρ²). Figure 9 shows the profile of the mass-weighted temperature over the same range as Figure 6 for t = 3.0 Gyr. The lack of temperatures higher than ∼3 keV demonstrates that the denser, colder gas provides the dominant contribution to the X-ray emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The single-temperature fits are in general agreement with our projected temperature profiles. Attempting to fit the cluster temperature using a two-temperature model resulted in no significant decrease in the χ²-statistic from the single-temperature fits, so on the basis of this analysis we cannot distinguish the two different average temperatures of the two clusters.

At late times there is still some moderate variation in the projected temperature profile with time, indicative of the continuing evolution of the system, though the average temperature of the cluster is roughly constant with time. In Ota et al. (2004), the observed isothermal temperature profile of Cl 0024+17 (reproduced in Figure 10) was taken as evidence that the merger was not recent (i.e., it was more than a few Gyr ago). Our results show that if viewed along the line of sight ∼2 Gyr after the merger, the temperature structure of the system can appear more regular in projection than it actually is (e.g., the structure seen in Figure 2). If the results of Jee et al. (2007) regarding a ring-like DM structure indicate the clusters are being viewed along the line of sight shortly after the merger, our results show that the system can still appear relatively relaxed. This agrees with previous results of mock X-ray observations of simulations of galaxy cluster mergers (e.g., Ricker & Sarazin 2001, Poole et al. 2006).

Zoom In Zoom Out Reset image size

Finally, our observed spectral temperature from our simulations of T∼ 2.6–2.8 keV is somewhat lower than the T ∼ 4.47 keV for Cl 0024+17 measured by Ota et al. (2004) or the T ∼ 4.25 keV measured by Jee et al. (2007). If the collision scenario for Cl 0024+17 is correct, this indicates that the initial conditions for our simulation are not identical to the pre-merger conditions for the real cluster. A higher impact velocity or a change in the masses of the clusters could account for the temperature difference.

Though the temperature and surface brightness profiles are undergoing considerable evolution even at late times, it is still relevant to ask how accurate a mass estimate would be under the assumption of hydrostatic equilibrium. As previously noted, Ota et al. (2004) and Jee et al. (2007) made estimates of the mass of Cl 0024+17 using single and double β-model fits to the surface brightness profile together with an isothermal temperature profile. Although our temperature profile is not strictly isothermal, existing methods for temperature deprojection assume spherical symmetry, so we do not attempt such a deprojection here. If we assume a cluster temperature T_spec given by the spectral temperature fit to the entire cluster, we can arrive at a rough estimate for the mass of the system under the assumption of hydrostatic equilibrium.

For a β-model fit to the cluster gas, the projected cluster mass within a cylindrical volume for a given temperature, β, and core radius r_c can be estimated by (Ota et al. 1998; Jee et al. 2005)

$$\begin{equation} M_{X,\beta }(R) = 1.78 \times 10^{14}\beta \left({T \over {\rm keV}}\right)\left({R \over {\rm Mpc}}\right){R/r_{\rm c} \over \sqrt{1+(R/r_{\rm c})^2}}{\it M}_\odot \,, \end{equation} \tag{ 6 }$$

where R is the projected radius, T is the gas temperature, and β is the β-model index. We compute the estimated mass within the arc radius of the lensed galaxy in Cl 0024+17, r_arc = 35''/223.65 kpc. For the double β-model fits we assume the same temperature for both cluster components, and we add the mass contributions from the two fits together. Table 5 shows the estimated masses from the single and double β-model fits for times t = 2.0, 2.5, and 3.0 Gyr. Errors are quoted at the 90% confidence level. It is clear from the table that assuming a single-cluster model for the system will underestimate the projected mass by a factor ∼2–3, whereas assuming a double-cluster model will estimate the projected mass to within ∼10%. This is in agreement with the results of Jee et al. (2007).

In our simulated head-on collision, the gas in the clusters is significantly disrupted and by the end of the simulation is still settling into the potential wells created by the DM. If instead of colliding head-on the encounter were allowed to be slightly off-axis, this would result in less disruption of the cluster gas cores. Presumably the gas would relax more quickly and an even more accurate mass determination could be made under the assumption of hydrostatic equilibrium. In addition, in a slightly off-axis collision viewed along the line of sight the two cluster components would be more clearly distinguished in the X-ray emission. Two-component fits to both the X-ray surface brightness and the temperature of the different components would enable a more accurate mass determination. This is a large parameter space of simulations that will be explored in future investigations, but is beyond the scope of this paper.

These results have important implications for X-ray surveys. Neighboring galaxy clusters viewed in superposition, whether undergoing mergers or not, may bias temperature measurements based on spectral fitting. This can have an effect particularly on X-ray-based scaling relations. For example, previous simulated observations of merging clusters have demonstrated that for mergers occurring along the line of sight, the clusters appear to have higher temperatures for their mass (Ricker & Sarazin 2001; Poole et al. 2007). Our investigation shows that clusters viewed in superposition can make a cluster appear colder than it actually is if the colder cluster's gas is significantly denser. In addition, discrepancies between hydrostatic and lensing mass estimates may in some cases be attributed to the existence of multiple components viewed along the line of sight. as is apparent in the case of Cl 0024+17. Detailed analysis of X-ray observations would be required to discern such superpositions.

One possible way to correct for and identify such superposition effects is by using the optical redshifts of the cluster galaxies to identify separate cluster components, as Czoske et al. (2001, 2002) did in the case of Cl 0024+17. To reliably identify such separate cluster components requires accurate redshifts. Because spectroscopic redshifts are only feasible for the nearest or brightest cluster galaxies, cluster surveys must rely on photometric redshift estimators. Large optical surveys such as zCOSMOS (Lilly et al. 2007), Combo-17 (Wolf et al. 2003), and the upcoming Dark Energy Survey (The Dark Energy Survey Collaboration 2005) have photometric redshift errors σ_z ∼ (0.02 − 0.05)(1 + z).

If the only redshift determination for Cl 0024+17 were via photometry, would the separation between the two components still be discernible? To answer this question we constructed a mock "galaxy" catalog from our simulation. We chose DM particles from the two clusters as proxies for the cluster galaxies (120 from the smaller cluster, 240 from the larger) and assumed a value for σ_z. We constructed redshift histograms for the resulting galaxy distributions and fitted them to both a Gaussian distribution and a sum of two Gaussian distributions:

$$\begin{equation} f_{\rm single} = A_0e^{-(z-{\mu _0})^2/2{\sigma _0^2}} \end{equation} \tag{ 7 }$$

$$\begin{equation} f_{\rm double} = A_1e^{-(z-{\mu _1})^2/2{\sigma _1^2}} + A_2e^{-(z-{\mu _2})^2/2{\sigma _2}}. \end{equation} \tag{ 8 }$$

The fitted parameters for these distributions for varying σ_z, given a choice of mock galaxies from our simulation, are shown in Table 6. For low values of σ_z, the mean values of the two redshift distributions are statistically distinguishable, but for higher values they are not. The histograms shown in Figures 11 and 12 demonstrate that for σ_z∼ 0.005–0.01 the skewness of the distribution of galaxies is evident, but when the redshift error increases to σ_z ≳ 0.02, the resulting galaxy distribution is indistinguishable from that of a single cluster.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At a redshift of z ≈ 0.4 (Cl 0024+17), the expected photometric redshift errors from optical surveys would be too large at this time to distinguish between the two cluster components (σ_z∼ 0.03–0.07), and at a redshift of z ≈ 1 it would be worse (σ_z ∼ 0.04–0.1). Identification of line-of-sight mergers and collisions from redshift distributions likely will be restricted to clusters for which we have spectroscopic redshifts available.

Our X-ray fitting results might be taken to suggest that X-ray surface brightness profiles could be used in place of galaxy redshift distributions to detect clusters seen in superposition. To test this hypothesis, we placed our simulated clusters at the epoch t = 3.0 Gyr (∼2 Gyr after the collision, where we can clearly distinguish both cluster components) at two redshifts, the current redshift of z = 0.395 and a redshift of z = 1.0. The same analysis of the surface brightness profiles as above was performed for exposure times of 10 ks, 40 ks, and 240 ks; the former being a more typical exposure time for a given cluster in an X-ray survey. We focus here on the difference in the χ² between the single and double β-model fits. The resulting Δχ² for the fitted surface brightness models corresponding to the different exposure times for the different redshifts is shown in Table 7. At z = 0.395, it is possible to use the Δχ² between the two models to distinguish between the two clusters for all exposure times. At z = 1.0, for an exposure time of 10 ks, the Δχ² is ≈14 for a change of 2 degrees of freedom, and for 40 ks Δχ² ≈ 70. Therefore, we find that at exposure times more typical of X-ray surveys, even at high redshift we can distinguish between the two cluster components, though with more difficulty. This is potentially a promising avenue for identifying merging clusters and clusters in projection. However, the ability to identify merging clusters requires that the two clusters have surface brightness profiles that are distinct enough that the likelihood-ratio test can distinguish between them. For example, earlier on in the simulation, at the epoch t = 2.0 Gyr (≈ 1.0 Gyr after the collision, where we cannot distinguish between the two clusters as clearly), we find that the test can distinguish between the two components for our fiducial exposure time and redshift, but not at lower exposures and higher redshifts (Table 7). Thus, determining exactly the physical circumstances under which the likelihood ratio test is the most powerful way to identify merging clusters and clusters in projection requires further study.