Simulating the Bullet Cluster

Abstract

We present high resolution N-body/smoothed particle hydrodynamics (SPH) simulations of the interacting cluster 1E0657-56. The main and the subcluster are modelled using extended cuspy Λcold dark matter (ΛCDM) dark matter haloes and isothermal β-profiles for the collisional component. The hot gas is initially in hydrostatic equilibrium inside the global potential of the clusters. We investigate the X-ray morphology and derive the most likely impact parameters, mass ratios and initial relative velocities. We find that the observed displacement between the X-ray peaks and the associated mass distribution, the morphology of the bow shock, the surface brightness and projected temperature profiles across the shock discontinuity can be well reproduced by offset 6:1 encounters where the subcluster has initial velocity (in the rest frame of the main cluster) 2.3 times the virial velocity of the main cluster dark matter halo. A model with the same mass ratio and lower velocity (1.5 times the main cluster virial velocity) matches quite well most of the observations. However, it does not reproduce the relative surface brightness between the bullet and the main cluster. Dynamical friction strongly affects the kinematics of the subcluster so that the low-velocity bullet is actually bound to the main system at the end of the simulation. We find that a relatively high concentration (c= 6) of the main cluster dark matter halo is necessary in order to prevent the disruption of the associated X-ray peak. For a selected subsample of runs we perform a detailed three-dimensional analysis following the past, present and future evolution of the interacting systems. In particular, we investigate the kinematics of the gas and dark matter components as well as the changes in the density profiles and the motion of the system in the L_X–T diagram.

methods: N-body simulations, galaxies: clusters: individual: 1E0657-56, dark matter, X-rays: galaxies: clusters

1 INTRODUCTION

The ‘Bullet Cluster’ 1E0657-56 represents one of the most complex and unusual large-scale structures ever observed. Located at a redshift z= 0.296 it has the highest X-ray luminosity and temperature of all known clusters as a result of overheating due to a recent supersonic Mach (Markevitch et al. 2002; Markevitch 2006) central encounter of a subcluster (the bullet) with its main cluster. The 500-ks Chandra ACIS-I image of 1E0657-56 (fig. 1 in Markevitch 2006) shows two plasma concentrations with the bullet subcluster on the right-hand side of the image being deformed in a classical bow shock on the western side as a result of its motion through the hot gas of the main cluster. The analysis of the shock structure leads to the conclusion that the bullet is now moving away from the main cluster with a velocity of ∼4700 km s⁻¹ (Markevitch 2006). The line-of-sight velocity difference between the two systems is only 600 km s⁻¹ suggesting that the encounter is seen nearly in the plane of the sky (Barrena et al. 2002). As the core passage must have occurred ∼0.15 Gyr ago we have the unique opportunity to study this interaction in a very special short-lived stage, far away from thermal and dynamical equilibrium. As a result of the encounter, the collision-dominated hot plasma and the collisionless stellar and dark matter components have been separated. The galaxy components of both clusters are clearly offset from the associated X-ray emitting cluster gas (Liang et al. 2000; Barrena et al. 2002). In addition, weak and strong lensing maps (Clowe, Gonzalez & Markevitch 2004; Bradač et al. 2006; Clowe et al. 2006) show that the gravitational potential does not trace the distribution of the hot cluster gas that dominates the baryonic mass but follows approximately the galaxy distribution as expected for a collisionless dark matter component. The likelihood to find such a high-velocity cluster encounter in a Λcold dark matter (ΛCDM) cosmology has recently been investigated by Hayashi & White (2006) using the Millenium Run simulation. According to the newest estimates from X-ray and gravitational lensing results the Hayashi & White (2006) likelihood becomes ∼0.8 × 10⁻⁷ (Farrar & Rosen 2007) which means that 1E0657-56 represents an extremely rare system in a ΛCDM universe. Recent numerical works (Milosavljević et al. 2007; Springel & Farrar 2007, hereafter SF07) have demonstrated that the relative velocity of the dark matter components associated with the main and the subcluster is not necessarily coincident with the speed of the bullet as inferred from the shock analysis. In details, Milosavljević et al. (2007) using a two-dimensional Eulerian code well reproduced the observed increase in temperature across the shock front with a dark halo velocity ∼16 per cent lower than that of the shock, while SF07 found even a larger difference between the shock velocity (∼4500 km s⁻¹ in their best model) and the speed of the halo (only ∼2600 km s⁻¹). Moreover, according to Milosavljević et al. (2007), due to a drop in ram-pressure after the cores' interaction the gas component of the subcluster can eventually be larger than that of its dark matter counterpart.

The simulations of SF07 represent the most complete three-dimensional numerical modelling of the 1E0657-56 system so far. Nevertheless they focus preferentially on the speed of the bullet but fail in reproducing the observed displacement of the X-ray peaks which represent an important indicator of the nature of the interaction. In particular, they do not obtain any displacement in the X-ray distribution associated with the main cluster suggesting that the baryonic component is suffering too little ram-pressure. Moreover, the concentrations used for the main cluster (and obtained by modelling the lensing data) are much smaller than those suggested by ΛCDM (Macciò et al. 2007) for haloes of similar masses.

The aim of this paper is to investigate the evolution of the Bullet Cluster in details using high resolution smoothed particle hydrodynamics (SPH) simulations. We quantify the initial conditions that are required in order to better reproduce its observed state and predict its subsequent evolution.

Our model allows us to determine in details the spatial, thermal and dynamical state of the dark matter and hot gas distribution in the Bullet Cluster.

The paper is organized as follows. Section 2 describes the adopted cluster models and orbital parameters. In Section 3, we perform a projected analysis of our simulations comparing it with the latest X-ray and gravitational lensing results. In Section 4, we select some significant models and study in details the three-dimensional kinematics and morphology of the interacting systems and their past and future evolution with time, as well as the motion of the main cluster along the L_X–T diagram.

2 MODELS

Both the main and the subcluster are two components spherical systems modelled assuming a cuspy dark matter halo and a distribution of hot gas in hydrostatic equilibrium within the global potential of the cluster. The dark halo has a NFW (Navarro, Frank & White 1997) profile:  

1

where ρ_crit is the critical density of the universe at the time of the halo formation, r_s is a scale radius and δ_c the characteristic halo overdensity. The virial mass M_vir and radius r_vir are related by: M_vir=Δ_virρ_crit(4π/3)r³_vir, where the density contrast Δ_vir is set equal to 200. The concentration parameter c=r_vir/r_s is assumed to be dependent on the halo mass (Macciò et al. 2007). The velocity distribution at a given point in space is approximated by a Gaussian, whose velocity dispersion is given by the solution of the Jeans equation at this point (Hernquist 1993). The distribution of hot gas follows an isothermal β-model (Cavaliere & Fusco-Femiano 1976) of the form:  

2

We take the asymptotic slope parameter β= 2/3 (Jones & Forman 1984) and r_c= 1/2r_s (Ricker & Sarazin 2001). The adopted gas fraction ranges from a minimum value of 12 per cent, comparable to the gas mass fraction provided by X-ray observations of galaxy clusters (Vikhlinin et al. 2006; McCarthy, Bower & Balogh 2007), to 17 per cent, consistent with the recent WMAP results (Spergel et al. 2007).

Assuming a spherically symmetric model, the temperature profile is determined by the condition of hydrostatic equilibrium by the cumulative total mass distribution and the density profile of the gas (Mastropietro et al. 2005):  

3

where M(r) is the total mass within the radius r, m_p is the proton mass and μ the mean molecular weight. We assume μ= 0.6 for a gas of primordial composition, which appears to be a reasonable approximation since the mean temperature of 1E0657-56 is T∼ 14 keV according to Markevitch (2006) and cooling is dominated by bremsstrahlung and almost independent of the metallicity. G and k_B are the gravitational and Boltzmann constants.

Masses are assigned to the models according to the weak and strong lensing mass reconstruction of Bradač et al. (2006). In particular, we assume that the inferred mass enclosed within the field of the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) (Bradač et al. 2006) is comparable to the total projected mass of our simulated system (calculated when the two centres of the mass distribution are at a distance similar to the observed one) within the same area. Since the ACS field represents only the central fraction of the area covered by the entire system, this mass constraint is strongly influenced by the concentration of the dark matter haloes (Nusser 2007). With a cosmologically motivated choice of c= 6 (Macciò et al. 2007) for the main cluster initial halo model, we can reproduce the lensing mass reasonably well adopting a main cluster total mass (within the virial radius) of ∼8.34 × 10¹⁴M_⊙ (Table 1), almost a factor 1.8 smaller than the mass obtained by fitting lensing data with extremely low concentrated (c < 2) NFW haloes where the inner density profile is much flatter than the one suggested by ΛCDM simulations. One of the simulations presented in this paper (6:1c4) adopts a main halo with a lower (c= 4) concentration value. We will see that in this case the X-ray intensity peak associated with the main cluster is easily destroyed during the interaction.

Table 1

Initial conditions of the simulations.

Run	Mvir(main) (1014M⊙)	c (main)	Mvir(bullet) (1014M⊙)	c (bullet)	fg	di (kpc)	b (kpc)	vi (km s−1)	v (km s−1)	Model
6:1b0	7.13	6	1.14	8	0.17	5000	0	5000	4286	Adiabatic
6:1v3000b0	7.13	6	1.14	8	0.17	5000	0	3000	2571	Adiabatic
6:1	7.13	6	1.14	8	0.17	5000	150	5000	4286	Adiabatic
6:1v3000	7.13	6	1.14	8	0.17	5000	150	3000	2571	Adiabatic
6:1v2000	7.13	6	1.14	8	0.17	5000	150	2000	1714	Adiabatic
3:1	7.13	6	2.4	7	0.17	5000	150	5000	3750	adiabatic
8:1	7.13	6	0.91	8	0.17	5000	150	5000	4445	Adiabatic
6:1v3000big	14.2	6	2.4	7	0.17	5000	150	3000	2571	Adiabatic
6:1c4	7.13	4	1.14	8	0.17	5000	150	5000	4286	Adiabatic
6:1lfg	7.13	6	1.14	8	0.12	5000	150	5000	4286	Adiabatic
3:1lfg	7.13	6	2.4	7	0.12	5000	150	5000	3750	Adiabatic
6:1c	7.13	6	1.14	8	0.12	5000	150	5000	4286	Cooling
3:1clfg	7.13	6	2.4	7	0.12	5000	150	5000	3750	Cooling
10:1vb0c2	15	2	1.5	7.2	0.17	3370	0	1870	2057	Adiabatic
10:1vb0c2nfw	15	2	1.5	7.2	0.17	3370	0	1870	2057	Adiabatic

Note. For each run dark matter virial mass (M_vir) and concentration c of the main and the subcluster are indicated. f_g and d_i are the adopted gas mass fraction and initial separation of the cluster models. b represents the impact parameter while v_i and v are the initial velocity of the subcluster in the system of reference where the main cluster is at rest and in the centre of mass rest frame, respectively. The last column indicates whether radiative cooling is included.

We model encounters with mass ratios 3:1, 6:1 and 8:1 between the main and the subcluster in order to investigate the effects of tidal and ram-pressure stripping, which significantly reduce the mass associated with the subcluster and lead to values closer to the 10:1 ratio inferred from lensing observations (Clowe et al. 2004, 2006; Bradač et al. 2006). Initial conditions similar to those of SF07 are also investigated. In particular, 10:1vb0c2 assumes an encounter between massive haloes with mass ratio 10:1, zero impact parameter, initial relative velocity v= 1870 km s⁻¹ and a main halo with a concentration c= 2. The baryonic component is modelled according to equation (2). Run 10:1vb0c2nfw has the same initial conditions but – as in SF07 – the gas follows the same NFW profile as its dark matter counterpart.

The main cluster is initially at rest and the subcluster moves in the x direction with a velocity which ranges from ∼1900 to 5000 km s⁻¹. The initial conditions of the different runs are summarized in Table 1. The velocities of the subcluster relative to the centre of mass of the system are listed in the last but one column of the table.

All the simulations were carried out using GASOLINE, a parallel SPH tree-code with multistepping (Wadsley, Stadel & Quinn 2004). Most of the runs are adiabatic, with γ= 5/3. Radiative cooling for a primordial mixture of hydrogen and helium in collisional equilibrium is implemented in 6:1c and 3:1lfgc (star formation is not activated). The main cluster is modelled with 1.8 × 10⁶ particles, 10⁶ SPH and 8 × 10⁵ collisionless. The subcluster, with the exception of run 8:1 (where the number of gas particles in the bullet is 4 × 10⁵), has 9 × 10⁵ particles, 5 × 10⁵ collisional and the remainder dark matter particles. The gravitational spline softening is set equal to 5 kpc for the gaseous and dark component.

3 PROJECTED ANALYSIS

Each numerical work which aims to simulate the Bullet Cluster should be able to reproduce simultaneously the main features observed in X-ray maps (the bow shock, the relative surface brightness of the bullet and the main cluster), and the observed surface brightness and projected temperature profiles across the shock discontinuity. An additional constraint is provided by the observed displacement between X-ray and lensing maps, which is not negligible in both the main and the subcluster (∼110 and ∼270 kpc, respectively according to Clowe et al. 2006).

3.1 X-ray morphology: b= 0

A first indication about the validity of a model arises from the qualitative comparison of our simulated X-ray surface brightness maps with the X-ray 500-ks Chandra ACIS-I images provided by Markevitch (2006).

The impact parameter b is not strictly constrained by observations. Nevertheless a head-on merger, with b= 0, seems to be excluded by comparing deep X-ray observations and weak lensing maps. In particular, in the top panel of Fig. 1– which is fig. 1(b) of Clowe et al. (2006)– the brightest gas associated with the main cluster is not located along the line which connects the centres of the two total mass distributions. Moreover, the X-ray emission from the main cluster is asymmetric, with a peak in the north of the image and an extended tail of less bright material pointing south. These features are hardly associable with a zero impact parameter interaction, as shown in the middle and bottom panels of Fig. 1 which illustrate two simulated 6:1 head-on encounters where the subcluster moves from the left- to the right-hand side of the image (x-axis of the simulation) with decreasing initial velocities (the middle panel corresponds to run 6:1vb0 of Table 1, the bottom one to 6:1v3000b0). Images are projected along an axis perpendicular to the collision plane (the encounter is seen face-on) and the selected snapshot (the time interval between two subsequent outputs is 40 Myr) is the one which most closely matches the observed distance between the centres of the total mass distributions, associated with the two clusters – about 720 kpc from Bradač et al. (2006)– once the subcluster has passed through the core of the main system. Hereafter we will refer to the time corresponding to this snapshot as to the present time.

Figure 1

Upper panel: 500-ks Chandra image of the system with weak lensing k reconstruction shown in green (courtesy of D. Clowe and reproduced by permission of the AAS). Central and bottom panels: 0.8–4 keV surface brightness maps of runs 6:vb0 and 6:v3000b0. Logarithmic colour scaling is indicated by the key at the bottom of the figure with violet corresponding to 10³⁸erg s⁻¹kpc⁻² and white to 2 × 10⁴¹erg s⁻¹kpc⁻². White contours trace the total surface mass density of the system within 2.3 × 10³ and 2.3 × 10⁸M_⊙ kpc⁻². The box size is 1800 kpc.

In the lower two panels of Fig. 1, colours represent X-ray maps in the Chandra energy band (0.8–4 keV) generated using the Theoretical Image Processing System (TIPSY), which produces projected X-ray surface brightness maps with the appropriate variable SPH kernel applied individually to the flux represented by each particle. Assuming complete ionization and zero metallicity (metal lines are expected to provide a significant contribution to emissivity only at relatively low temperatures, less than 2 keV) the X-ray luminosity in a given energy band is defined as (Borgani et al. 2004)  

4

where Λ(T_i) is the cooling function in the specific band, T_i, ρ_i and m_i are temperature, density and mass associated with the ith hot (T_i > 10⁵K) gas particle, respectively, m_p is the proton mass and μ= 0.6 the mean molecular weight. The sum runs over all the N_gas particles within an oblong of base equal to the pixel size and major axis oriented along the line of sight. When the X-ray luminosity of the entire cluster is calculated N_gas is the number of hot particles within the virial radius r_vir. The cooling function is computed using a Raymond–Smith code (Raymond & Smith 1977) for a gas of primordial composition. The energy band (0.8–4 keV) is chosen in such a way to reproduce the 500-ks Chandra ACIS-I image of the Bullet Cluster (Markevitch 2006). The entire energy band, used to calculate bolometric X-ray luminosities in the following of the paper, goes from 5 eV to 5 × 10⁴ keV.

The lower two panels of Fig. 1 show that after an encounter with zero impact parameter the displacement of the main cluster's X-ray peak is aligned with the x-axis. A large relative velocity (central panel) induces a significant offset between the dark and baryonic component of the main cluster (see Table 2 for details) but it also leads to substantial disruption of the main cluster gaseous core. Moreover, the displacement of the bullet from its dark halo (278 kpc) is much larger than observed. Decreasing the relative velocity between the two clusters (bottom panel) two X-ray peaks are clearly visible but the displacement of the main cluster gas is now negligible due to the lower ram-pressure experienced by the main cluster core. SF07 (fig. 7 in their paper) provide further examples of head-on encounters with even lower mass ratios and relative velocities (v= 2600 km s⁻¹ in the centre of mass rest frame). Even assuming extremely low concentrations (c= 2) for the main halo, the authors never reproduce the displacement observed in the two systems. Increasing the concentration strongly increases the luminosity of the main cluster, which appears much brighter than the bullet, contrary to what is observed.

Table 2

Present time.

Run	Δ	offset (bullet)	offset (main)	vgas	vdark
6:1b0	753	278	188	3215	4715
6:1v3000b0	741	213	66	3131	3134
6:1	742	237	128	3609	4756
6:1v3000	729	185	172	2893	3137
6:1v2000	721	126	230	2849	2425
3:1	784	162	223	3908	4076
8:1	737	228	117	3647	4858
6:1v3000big	725	139	92	3927	3528
6:1c4	735	151	192	4168	4799
6:1lfg	718	200	189	3811	4804
3:1lfg	779	197	242	3746	4145
6:1c	736	234	127	3497	4806
3:1clfg	780	228	272	3595	4205
10:1vb0c2	757	22	119	3224	2718
10:1vb0c2nfw	756	26	35	3284	2507

Note. Δ is the projected (perpendicular to the plane of the encounter) distance between the peaks of the total mass distributions, associated with the two clusters. The third and fourth columns represent the projected offset between each X-ray peak and the associated mass density peak. v_gas and v_dark are the subcluster gas and dark matter velocity calculated in the centre of mass system of reference.

A bow shock is clearly visible on the right-hand side of each image. The shape of the shock front is only marginally dependent on the kinematics of the model while the distance between the edge of the bullet (the so called contact discontinuity) and the shock front becomes larger for decreasing bullet velocities. The contact discontinuity itself is much flatter in the case of 6:1b0 than in the low-velocity encounter 6:1v3000b0 and clearly not comparable to observations, which show a more narrow structure. In general, a more efficient ram-pressure during the phase of core–core interaction is associated with a larger opening angle of the contact discontinuity at the present time (Quilis & Moore 2001).

3.2 X-ray morphology: b > 0

The rest of the runs listed in Table 1 have an impact parameter b equal to 150 kpc, comparable to the core radius r_c of the main cluster gas distribution for most of the models. For a plasma distributed according to a β profile like the one adopted in this paper the assumption b= 150 kpc implies that the maximum external density crossed by the core of the subcluster is 1.5 times smaller than the one it would pass through for b= 0 kpc. For the same initial velocity, the subcluster also sees a weaker potential with respect to the case of a head-on encounter and experiences a smaller maximum orbital velocity. A choice of a much larger value of b would further decrease the mutual ram-pressure between the two systems and require much higher relative velocities in order to explain the observed offset between the gas and dark matter.

In Figs 2 and 3, we illustrate the projected X-ray surface brightness maps of some interesting models with impact parameter b= 150 kpc. The encounters are shown face-on. The box size and the surface mass density contours are the same as in the last two panels of Fig. 1. In order to underline the morphological details of the high emission regions the upper limit of the surface brightness scale varies in the different images. Individual values are indicated in the captions.

Figure 2

0.8–4 keV surface brightness maps of runs (from the top left-to-bottom right-hand panels) 3:1, 8:1, 6:1, 6:1v3000, 6:1v2000 and 6:1c4. Logarithmic colour scaling is indicated by the key at the bottom of the figure, with violet corresponding to 10³⁸erg s⁻¹kpc⁻² and white to 1.8 × 10⁴¹erg s⁻¹kpc⁻² in runs 8:1 and 6:1c4, to 2.34 × 10⁴¹erg s⁻¹kpc⁻² in run 1:3 and to 2 × 10⁴¹erg s⁻¹kpc⁻² in the remaining cases. Projected isodensity contours of the total mass distribution are shown. Limits are 2.3 × 10³ and 2.3 × 10⁹M_⊙ kpc⁻². Each box size is 1800 kpc.

Figure 3

Same as in Fig. 2 for the cooling run 6:1c (the X-ray upper limit is 1.8 × 10⁴¹erg s⁻¹kpc⁻²) and run 6:1v3000big (3.9 × 10⁴¹erg s⁻¹kpc⁻²).

Among this subsample of runs, 3:1 produces the largest displacement of the X-ray peak associated with the main cluster, but the X-ray map differs from the observations. In particular a large strip of strongly emitting gas still connects the two X-ray peaks while the morphology of the main cluster peak is much more elongated than observed. Increasing the mass ratio between the two interacting systems (runs 6:1 and 8:1) the displacement in the bullet becomes larger than that in the main cluster. Run 6:1 is characterized by an initial subcluster velocity of 5000 km s⁻¹ (as well as runs 3:1 and 8:1) in the system of reference where the main cluster is at rest, that corresponds to a present-time velocity of ∼4300 km s⁻¹ in the centre of mass rest frame. The same model is simulated assuming lower relative velocities (6:1v3000 and 6:1v2000). With decreasing velocities the offset in the subcluster becomes smaller and the X-ray emission from the bullet less pronounced with respect to the bright X-ray emitting region at the centre of the main cluster. At the same time the shape of the contact discontinuity changes, getting progressively more narrow while the distance between the contact discontinuity and the shock front increases as will be shown more quantitatively in the next section. The displacement associated with the main-cluster is determined by the distance of closest approach between the centres of the two clusters, which becomes smaller – assuming the same initial impact parameter b– with decreasing bullet velocities. Indeed, the separation between the main cluster X-ray emission peak and its dark matter counterpart is maximum in the case of the low-velocity run 6:1v2000.

A non-negligible fraction of the X-ray emission visible at the present time near the centre of the main-cluster is actually associated with hot gas stripped from the external regions of the bullet. Fig. 4 refers to run 6:1v3000. It shows the individual projected distribution of gas originating from the main (top panel) and subcluster (bottom panel) and lying at the present time within 1 Mpc from the centre of the system. Comparing these images with the middle right-hand panel of Fig. 2 it appears evident that the bright, elongated X-ray feature crossing the second innermost isodensity contour is associated with the displaced gaseous centre of the main cluster while the surrounding, more diffuse region hosts a significant amount of subcluster gas. Indeed, as illustrated in the upper panel of Fig. 4, the motion of the bullet across the inner regions of the main system creates a low density ‘tunnel’ in the main cluster gas distribution. At the same time the subcluster loses a large amount of baryonic material during the phase of core–core interaction. This material, which fills the tunnel, falls back into the gravitational centre of the main cluster and resides now more than 500 kpc distant from the X-ray bullet. As we will show in the next section, the amount of gas deposited by the subcluster in the central regions of the main system increases with decreasing relative velocities. This trend explains the relative increase in luminosity of the diffuse strongly emitting component if compared to the peak associated with the main cluster core gas (always in the upper right-hand region with respect to the centre of the mass density distribution), when we compare run 6:1 with the low-velocity run 6:1v2000 where it becomes the primary peak of X-ray emission.

Figure 4

Run 6:1v3000. Gas originating from the main (top panel) and the subcluster (bottom panel) is projected individually along the z-axis perpendicular to the plane of the encounter. Violet corresponds to a surface density of 2.3 × 10⁵M_⊙ kpc⁻² and white to 2.3 × 10⁸M_⊙ kpc⁻². Projected isodensity contours of the total mass distribution are drawn on top of the image. The box size is 1 Mpc.

As shown in the bottom right-hand panel of Fig. 2, a main halo with low concentration (c= 4) does not survive a 6:1 subcluster encounter with velocity v= 5000 km s⁻¹ and its X-ray peak is destroyed.

As previously noticed by SF07 the choice of a lower gas fraction (f_g= 0.12 in three of the last four runs of Table 1) does not affect significantly the X-ray map morphology although the displacement of the two luminosity peaks with respect to their dark matter counterparts changes slightly.

Finally, we tested the consequences caused by including radiative cooling and choosing a larger main halo model (Fig. 3). Cooling makes the contact discontinuity narrower and the amount of diffuse X-ray gas around the peak associated with the main cluster smaller. The offset in the main and subcluster remains unaltered. Run 6:1v3000big is characterized by a main cluster total mass of ∼1.64 × 10¹⁵M_⊙, which is closer to the value adopted by SF07 and predicted by fitting the large field weak lensing data with extremely low concentrated (c∼ 2) NFW haloes. Although the initial relative velocity is only 3000 km s⁻¹, due to the large mass of the host halo the present-time velocity of the bullet in the centre of mass rest-frame is much higher than in the corresponding 6:1v3000 run. Consequently the offset of the X-ray peak associated with the main cluster is less than 100 kpc and the amount of gas lost by the bullet in the core of the main system is closer to that observed in 6:1 than in 6:1v3000.

Summarizing, run 6:1v3000, with mass ratio 6:1, initial relative velocity v= 3000 km s⁻¹ and present-time subcluster velocity v∼ 3100 km s⁻¹ in the centre of mass rest frame, best reproduces the main features observed in X-ray maps, in particular the peculiar morphology of the X-ray emission associated with the main cluster, the relative surface brightness between the main and the subcluster, the shape of the shock front and of the contact discontinuity. In the rest of the paper we will indicate this model as our favourite model. Although the low-velocity run 6:1v2000 leads to a X-ray displacement closer to that observed by Clowe et al. (2006) both in the main and in the subcluster, this model seems to be excluded on the basis of a pure morphological comparison with the observational data. Indeed, the bullet seems to be much less bright than the centre of the main cluster and the morphology of the main cluster X-ray peak itself does not match observations.

3.3 X-ray morphology: comparison with SF07

In runs 10:1vb0c2 and 10:1vb0c2nfw we adopted initial conditions similar to those of SF07. The two runs differ in the gas density profile: in 10:1vb0c2 the baryonic component follows a β-profile, while in 10:1vb0c2nfw it has the same NFW profile as the dark matter, except for a small core in the centre (100 kpc in the main cluster and 50 kpc in the bullet).

Their X-ray maps at the present time are shown in Fig. 5. In 10:1vb0c2 (left-hand panel) – despite the initial very low concentration of the main halo – the peak associated with the main cluster gas survives to the core–core encounter due to the relatively low initial velocity and high mass ratio of the interaction. It gets displaced by ∼120 kpc with respect to the centre of the mass density distribution while the offset in the bullet is almost negligible. Note that the offset in the main cluster in 10:1vb0c2 is larger than in 6:1v3000b0, although the latter run is characterized by stronger ram-pressure values. Indeed, the low concentration of the main halo makes it simpler to displace the baryonic component from the centre of the dark matter potential. The choice of a different initial gas density distribution influences the intensity and position of the main cluster X-ray peak after the central phases of the interaction. Gas in the right-hand panel of Fig. 5 is initially set in hydrostatic equilibrium with NFW profiles. We note that in this case the main cluster does not show any significant displacement, while its luminosity is much lower than that of the bullet. The displacement in the subcluster appears negligible as well. The present-time configuration is similar to that found by SF07 for an encounter with concentration of the main cluster c= 2. The shape of the contact discontinuity appears rounder than in SF07 (compare with their fig. 2). This could be due to a resolution effect (SF07 have a gas mass resolution 1.5 higher than in our simulations) or to the choice of the core size in the gas density profile of the bullet. Indeed the linear dimension of the contact discontinuity is slightly more than twice the gas core size.

Figure 5

Same as in Fig. 2 for runs 10:1vb0c2 (left-hand panel) and 10:1vb0c2nfw (right-hand panel) (the X-ray upper limit is 2 × 10⁴¹erg s⁻¹kpc⁻²).

3.4 Projected mass

Fig. 6 compares the integrated mass profiles of the main and the subcluster for different mass models with the gravitational lensing results from Bradač et al. (2006) and previous simulations by SF07. The profiles inferred from observations are obtained by measuring the enclosed mass in cylinders centred on the southern cD of the main cluster and on the BCG of the subcluster. Here we calculated the present-time projected mass within cylindrical bins centred on the centre of mass of the two clusters. The massive cluster run 6:1v3000big (green lines) overestimates the cumulative profile of the main cluster, while an initial main cluster mass of ∼8.3 × 10¹⁴M_⊙ (6:1v3000, in red) fits the observational points better than the large-mass–low-concentration model of SF07 in the central regions and underestimates the projected mass by less than 20 per cent between 250 and 450 kpc from the centre. Blue lines refer to run 3:1lfg where the main halo has the same dark matter mass as 6:1v3000 but lower baryonic fraction. It is interesting to notice that using almost the same initial subcluster mass (the subcluster model in run 6:1v3000big and 3:1lfg only differs in the baryonic fraction) we can nearly reproduce the observed projected mass in the bullet only in the case of a very massive host halo (green dashed line, 6:1v3000big run) which overpredicts the main cluster profile (green solid line). On the other hand a smaller main halo – which better fits the observational data (blue solid line) – associated with a similar subcluster mass (3:1lfg) underpredicts the lensing results (blue dashed line). Our favourite run 6:1v3000 gives even lower projected subcluster mass estimates (red dashed line) and lies close to the model proposed by SF07 (their subcluster mass is similar to the one used in our 6:1 models). These results seem to indicate that the mass of the subcluster is a significant part of the main cluster mass, with initial mass ratio between the main and the subcluster much smaller than 3:1. Nusser (2007) successfully reproduces the projected mass associated with the subcluster only assuming very massive haloes and small mass ratios. In particular, the highest mass ratio allowed by their models is 2.7:1, with a main cluster and bullet mass of 3.2 × 10¹⁵ and 1.2 × 10¹⁵M_⊙, respectively. Such high masses are however clearly incompatible with the values inferred by galaxy kinematics. In particular, the virial mass derived by Barrena et al. (2002) for the main cluster is 1.24 × 10¹⁵ (comparable to the value obtained by Girardi & Mezzetti 2001). They also measure the main cluster total R-band luminosity and find L_R= 10¹²L_⊙. Using L_R to estimate the cluster mass (Miller et al. 2005) we obtain M₂₀₀∼ 10¹⁵M_⊙. The bullet mass determination of Barrena et al. (2002) is less reliable since it is based on only seven galaxies and on the assumption of equilibrium while the subcluster seems to be tidally perturbed beyond ∼200 kpc from its centre. They find M₂₀₀= 1.2 × 10¹³M_⊙. From the total R-band luminosity (L_R= 0.2 × 10¹²L_⊙) we get a value almost five times larger (M₂₀₀∼ 5 × 10¹³M_⊙). The bullet mass proposed by Nusser (2007) is even one order of magnitude larger than that obtained from weak lensing analysis (1.5 × 10¹⁴M_⊙) of a new larger lensing field which covers most of the area occupied by the system (Clowe et al., in preparation). With such low mass ratios we cannot reproduce the morphology of the X-ray maps (upper left-hand panel of Fig. 2) since the luminosity peak associated with the main cluster core is almost completely destroyed by the interaction.

Figure 6

Cumulative projected mass profiles. Points represent the observational results (Bradač et al. 2006) for the main (dots) and the subcluster (stars). Black curves refer to the model of SF07 (solid and dashed black line indicates the main and subcluster mass, respectively). Coloured lines refer to the models in this work, always with solid curves indicating the main cluster and dashed ones the bullet. In details the red lines are associated with our favourite model 6:1v3000, blue to 3:1lfg and green to 6:1v3000big.

Actually any attempt to derive the mass of the main and subcluster simply by projecting the mass of two isolated NFW haloes and fitting the observed lensing data – like the analysis of Nusser (2007)– is going to provide wrong results. Indeed, as we will show in Section 4, the dark matter profiles of the main and subcluster change significantly during the central phases of the interaction, with a density increase in the inner regions and a decrement beyond 0.4 r_vir. This effect, which characterizes both the main and the subcluster, is not related to stripping processes but can be explained by the fact that the two haloes are strongly perturbed and cannot be idealized as the sum of two isolated equilibrium models. Fitting the projected cumulative mass around the centres of the two interacting clusters with NFW models we would overestimate either the virial mass or the concentration of the original clusters.

3.5 Inclination effects

According to Barrena et al. (2002), the line-of-sight velocity difference between the main and the subcluster is relatively small, about 600 km s⁻¹, which implies that the encounter must be occurring nearly in the plane of the sky. In Fig. 7, we show the X-ray surface brightness map of run 6:1v3000 when the orbital plane is viewed inclined by 45° left (top panel) and down (bottom panel). In the upper image the bullet (which is closer to the observer) appears to be much brighter than the main cluster while the contact discontinuity and the shock front are much rounder than in the case where the encounter is seen perpendicular to the orbital plane. Such a large inclination has to be excluded by a comparison with observations. Moreover, the line-of-sight velocity difference between the two clusters (we assume the line-of-sight kinematics of the dark matter to be coincident with that of the stellar component observed by Barrena et al. 2002) is much larger than the observed one (∼2700 km s⁻¹) due to the fact that the main component (∼3200 km s⁻¹ in the centre of mass rest frame) of the subcluster velocity is oriented in the direction of the positive x-axis in the original face-on encounter. A velocity difference comparable to the observed one is reached for a right-to-left-hand inclination of ∼10°. On the opposite, a top-down inclination of the plane of the encounter with respect to the line of sight (bottom panel) does not produce drastic changes in the morphology of the bow shock or of the edge of the bullet. Indeed the main and the subcluster still lie close to the plane of the sky, but the relative surface brightness between the dense and the diffuse strongly emitting components at the centre of the main system changes. In particular the peak of X-ray emission due to the presence of hot gas deposited by the subcluster within the core of the main cluster is now comparable to the main cluster X-ray peak. The difference in line-of-sight velocity in this case is only ∼100 km s⁻¹ and would probably not be distinguishable with respect to a pure face-on encounter.

Figure 7

X-ray surface brightness maps for different inclinations of run 6:1v3000. In the top panel the orbital plane is on the left-hand side rotated by 45° with respect to the projection plane while in the panel on the bottom the encounter is seen with a top-down inclination of 45°. The colour scaling is the same as in the corresponding face-on map (Fig. 2, middle right-hand panel). The box size is 1.8 Mpc.

3.6 Brightness and temperature profiles across the shock discontinuity

Fig. 8 represents the X-ray brightness profile measured in a narrow slit (with thickness 20 kpc) parallel to the x-axis across the shock front. The contact discontinuity (the ‘edge’ of the bullet) is located at x= 0 with x increasing towards the pre-shock region. The general trend does not depend on the model and is comparable to the profile suggested by observations (Markevitch 2006): all the simulations show an inner bump associated with the bullet and an outer one (between x= 20 and 100 kpc) caused by the shock. The abrupt jump at x∼ 100 kpc is the shock front while the pre-shock region is well fitted by a two-dimensional β-profile S_X(x) =S_X0[1 + (x/x_c)²]^−3β+1/2. In the case of run 6:1 x_c= 150 kpc, β= 2/3 and S_X0= 8 × 10⁴⁰erg s⁻¹kpc⁻². A different orientation (45°, bottom-up inclination of the plane of the encounter) of the line of sight does not affect substantially the surface brightness profile. If cooling is activated the surface brightness of the pre and post-shock region is smaller, but the thickness of the shock-front and the jump in surface brightness are similar.

Figure 8

Runs 6:1 and 6:1c: X-ray surface brightness profiles across the shock discontinuity. The bullet is located at x= 0 kpc and the shock front at ∼100 kpc.

In order to calculate projected temperatures we need to define a weighting function. The emission weighted temperature T_ew was originally introduced to provide a better comparison between simulations and observations with respect to a simple mass weighted temperature definition and has been commonly used in the analysis of simulations (Borgani et al. 2004). It assumes a weighting function proportional to the emissivity of each hot gas particle and is defined as  

5

Recently Mazzotta et al. (2004) have demonstrated that for clusters with a complex thermal structure the emission weighted temperature always overestimates the spectroscopic temperature obtained from X-ray observations due to fact that the source is not a single temperature plasma. For clusters with a temperature T > 3 keV the discrepancy is ∼20−30 per cent. This difference becomes particularly large in the presence of strong temperature gradients, such as shocks that appear to be much weaker in observations than what is predicted by emission-weighted temperature maps (Mathiesen & Evrard 2001; Gardini et al. 2004; Rasia et al. 2005).

Mazzotta et al. (2004) proposed a new definition of temperature, the spectroscopic-like temperature T_sl:  

6

When applied to clusters hotter than 2–3 keV, this equation, with α= 0.75, gives a good approximation (within few per cent) of the spectroscopic temperature obtained from data analysis of Chandra.

The Bullet Cluster 1E0657-56 is characterized by the highest luminosity and temperature and the strongest bow shock of all known clusters (Markevitch 2006) and seems indeed to be the ideal candidate for adopting the spectroscopic-like temperature definition. Therefore in the remainder of this paper T_sl will be analysed and T_ew only indicated for a comparison. Fig. 9 illustrates the projected spectroscopic-like temperature profiles across the shock for the different runs of Table 1. x= 0 corresponds to the position of the contact discontinuity while the pre-shock gas is located at x > 0. All the values refer to the present time. The upper left-hand panel shows the temperature jumps associated with different relative velocities of the two clusters. Decreasing the initial relative velocity from 5000 to 2000 km s⁻¹ (in the system of reference where the main cluster is at rest) reduces the temperature peak by ∼7 keV while the peak itself becomes broader since the thickness of the shock region, between the shock front and the contact discontinuity, increases by almost a factor of 2 due to the lower pressure exercised by the pre-shock gas. Both the 6:1v3000 and 6:1v2000 models seem to fit quite nicely the observed height (∼27–30 keV) and thickness (150–200 kpc) of the shock front (Markevitch 2006) while the 6:1 run produces a peak which is too narrow (∼100 kpc) and pronounced (∼35 keV). The upper right-hand panel refers to different mass ratios. If the encounter is characterized by the same initial relative velocity and gas fraction the strongest shock is associated with the most massive subcluster. A 3:1 adiabatic encounter produces a ∼45 keV temperature peak with thickness ∼150 kpc while the 8:1 temperature profile is not substantially different from that of the 6:1 run. In the same plot we also show the shock created in the massive run 6:1v3000big. Clearly the maximum temperature is much higher than the observed one. Including radiative cooling (bottom left-hand panel) has the effect of reducing the peak in temperature but does not influence the thickness of the shock region. The adoption of a simple cooling model without non-gravitational heating is actually questionable. Indeed, although the estimate of the temperature jump is reasonable, the entire temperature profile drops by 5 keV due to the fact that once cooling is activated the main cluster gas component becomes thermally unstable in the early phases of the interaction and overcools in the central regions. If cooling is important, models with smaller mass ratios and higher relative velocities (like the 3:1lfg run in plot) have still to be taken in account and cannot be excluded as the high-temperature peaks could actually cool significantly. The choice of a different line of sight (bottom right-hand plot) does not affect significantly the temperature profile across the shock. Even decreasing the baryonic fraction in the clusters and assuming a much less concentrated main halo, the height and thickness of the temperature peak do not change.

Figure 9

Spectroscopic-like temperature profiles measured in a narrow slit (20 kpc) across the shock. The bullet is located at x= 0 kpc with x increasing towards the pre-shock region. Upper left-hand panel: 6:1 runs with different relative velocities are compared. Upper right-hand panel: different mass ratios. Bottom left-hand panel: comparison among adiabatic and cooling 6:1 and 3:1 encounters. Bottom right-hand panel: effects of inclination, lower gas fraction and lower concentration of the main halo.

In Fig. 10, we compare the spectroscopic-like temperature profile across the shock region with the emission weighted one for our favourite run 6:1v3000. While the projected temperature profile calculated according to the two definitions is similar (T_ew is only slightly higher than T_sl) in the pre and post-shock regions, the emission weighted temperature T_ew in the region 0 ≤x≤ 150 kpc is ∼20 per cent higher than T_sl. T_shock represents the actual temperature along the x-axis through the shock. Indeed the blue curve in Fig. 10 gives the exact temperature jump across the shock (the ‘true temperature’ of the shock) which is characterized by an even higher peak with respect to the projected ones. This deprojected temperature profile is actually the one shown by SF07 in their Fig. 9 where they compare their model with observations which on the other hand refer to projected quantities (Markevitch 2006). As the calculated local temperatures of SF07 fit the projected observed temperature profile very well we conclude that their projected temperature profiles are actually inconsistent with the observations by Markevitch (2006).

Figure 10

Run 1:6v3000: spectroscopic-like (T_sl), emission weighted (T_ew) and true (T_shock) temperature profiles across the shock.

3.7 Evolution around the present time

In Figs 11 and 12, we show for our selected adiabatic run 6:1v3000 the evolution with time of the 0.8–4 keV X-ray surface brightness and spectroscopic-like temperature during the central phases of the interaction. All the quantities are projected along an axis perpendicular to the plane of the collision. On top of the colour maps we draw in black the projected isodensity contours of the total mass distribution, which is dominated by the dark matter component. Time increases from the left-to-right-hand and from the top to the bottom. The sequence of six panels covers an interval in time of 400 Myr, between ∼1.3 and 1.7 Gyr from the beginning of the simulation. The snapshot corresponding to the present time is the one on the bottom left-hand panel. The bullet approaches the main cluster from the left-hand side, with an initial velocity of 3000 km s⁻¹ in the main cluster rest frame. The shock front has an arc-like shape (Ricker & Sarazin 2001) and becomes progressively more asymmetric as the bullet moves closer to the core of the main cluster. When the centres of the two clusters are less than 250 kpc apart, ram-pressure becomes effective in producing a displacement (visible in the projected density maps as well as in surface brightness) between the gaseous core of the bullet and the peak of its associated mass distribution. The offset in the main cluster is evident only when the bow shock passes through its core. In the central panels of Fig. 11 the X-ray luminosity saturates in order to distinguish features in the maps at later times. Nevertheless, it is visible how the main cluster core gets compressed and displaced from the centre of the potential towards the top right-hand side of the image and appears in the X-ray maps (last two panels) as an elongated structure, characterized by high surface brightness.

Figure 11

Time evolution of run 6:1v3000 (from the top left-to-bottom right-hand panels: 1.28, 1.36, 1.44, 1.52, 1.6, 1.68 Gyr from the beginning of the simulation). Colours represent X-ray (0.8–4 keV) surface brightness maps. Logarithmic colour scaling is indicated by the key to the bottom of the figure, with violet corresponding to 10³⁸erg s⁻¹kpc⁻² and white to 2.5 × 10⁴¹erg s⁻¹kpc⁻². Projected isodensity contours of the total mass distribution are drawn on top of temperature maps.

Figure 12

Time evolution of run 6:1v3000 (from the top left-to-bottom right-hand panel: 1.28, 1.36, 1.44, 1.52, 1.6, 1.68 Gyr from the beginning of the simulation). Colours represent spectroscopic-like X-ray temperature maps. Logarithmic colour scaling is indicated by the key to the bottom of the figure, with violet corresponding to 0.9 keV and white to 86 keV. Projected isodensity contours of the total mass distribution are drawn on top of temperature maps.

Temperature maps better describe the evolution of the shock region, which gets compressed and hotter during the core passage while in a later phase it cools down and becomes thicker due to lower pressure of the pre-shock gas. The bullet itself expands as it leaves the central regions of the main cluster. As observed by Markevitch (2006) and previously noticed by SF07, despite its strong X-ray emission (it is the brightest feature in the post core–core interaction X-ray maps) the bullet remains relatively cold. Even if radiative cooling is not activated, the core of the subcluster is heated only to a maximum temperature of 10⁸ K, while the shock front is much hotter (5 × 10⁸ K). The projected temperature associated with the bullet is higher when it passes through the core of the main system (central panels of Fig. 12) due to the line-of-sight overlap with the hot gas from the main cluster and the shock heated material stripped from the bullet itself which surrounds it. As soon as the subcluster moves out into the cool external regions of the main cluster and loses part of the hot envelope of stripped gas, its projected temperature decreases to values comparable to the observed ones. Another peculiar feature in the temperature maps is the high-temperature region next to the innermost total density contour of the main cluster and visible in the middle right- and bottom left-hand panels of Fig. 12. This area could be associated with the south-eastern high-temperature region observed by Markevitch et al. (2002) (regions f, i and l in their fig. 2) in the main cluster X-ray map and assumed to be coincident with the main merger site. In our simulations the high-temperature region is filled with hot gas stripped from the external regions of the subcluster and deposited within the core of the main system (compare with Fig. 4). This high-temperature material, combined with the main cluster gas which lies in the same projected region, produces the diffuse X-ray high emission feature visible at the present time in the main cluster below the primary peak (bottom left-hand panel of Fig. 11 or middle right-hand panel of Fig. 2 for a better colour contrast). The eastern side of the high-temperature region is less bright in X-ray emission and lies beyond the luminosity peak associated with the main cluster. As the subcluster moves to larger radii, the high-temperature region expands and cools. Shortly (∼70 Myr) after the present time this gas is no longer clearly distinguishable in X-ray maps.

4 THREE-DIMENSIONAL ANALYSIS

In this section we will focus on a sample of runs from Table 1– namely 6:1, 6:1v3000, 6:1v2000 and 3:1clfg – and investigate in detail their three-dimensional characteristics. Three of these simulations are adiabatic models characterized by different initial relative velocities and will permit us to study the effects of the subcluster speed on the present and future state of the encounter. The remaining run 3:1clfg is the one which better reproduces – together with 6:1v3000 and 6:1v2000 – the observed jump in temperature across the shock front.

4.1 The shock discontinuity

Fig. 13 illustrates the behaviour of hydrodynamical quantities across the shock discontinuity for a snapshot which corresponds to the present time. The horizontal axis is centred on the gaseous bullet and oriented perpendicular to the bow shock nose. The bow shock location is indicated by a vertical line while the edge of the bullet corresponds to the peak in density. The subcluster is moving outward from the main cluster core towards positive values of x. All the physical quantities are calculated and mass averaged on a one-dimensional grid where each grid element has a volume of 15³kpc³⋅v is the component of the velocity perpendicular to the shock front and is calculated with respect to the system centre of mass. The Mach number is determined from the temperature jump – which shows a better defined discontinuity compared with the density jump – using the Rankine–Hugoniot conditions. We find values in good agreement with observations (Markevitch 2006) for 6:1v3000 and the cooling run 3:1clfg (both with ⁠). 6:1v2000 is characterized by a slightly lower value of while the high-velocity run 6:1 has a stronger shock, with ⁠.

Figure 13

Gas temperature, density and one-dimensional velocity profiles across the shock discontinuity. The bullet is located at x= 0 while the dashed vertical lines indicates the position of the bow shock.

As seen in the previous section, the pre-shock temperature (∼11–12 keV for the adiabatic runs) is slightly higher than the one obtained by projection along the line of sight (Fig. 10). Assuming T= 11 keV and we predict a pre-shock sound speed c_S= 1700 km s⁻¹ and a shock velocity ⁠. This value reduces to the observed shock velocity v_S∼ 4700 km s⁻¹ if we use the projected average pre-shock temperature (T∼ 9 keV) adopted by Markevitch (2006).

As previously noticed by other authors (SF07) the velocity jump shown in the last panel of each plot is much smaller than the theoretically inferred shock velocity. Actually the pre-shock gas is not at rest but shows a negative velocity along the x-axis which however can only be partially explained by the fact that the centre of mass of the system is moving in the positive direction of the x-axis following the bullet. Indeed the upstream velocity maintains a negative sign even with respect to the rest frame of the parent cluster, indicating a pre-shock infall towards the bullet. This effect is explained in SF07 by studying the dynamical evolution of the system's global potential, which becomes deeper after the core–core interaction and induces infall of material from the region ahead of the shock. Nevertheless the infall velocity which characterizes our models is significantly smaller than the subcluster velocity in contrast to the values found by SF07 and ranges only between 500 and 900 km s⁻¹ in the different runs.

4.2 Dark matter, gas and shock velocities

Figs 14 and 15 illustrate the characteristic velocities of the subcluster in the orbital plane for the three runs 6:1v3000, 6:1 and 6:1v2000. All the velocities are calculated in the centre of mass rest frame and the time corresponding to the present position is indicated by a vertical line. The velocity of the dark matter component peaks at the moment of closest approach between the two cores and then decreases faster than for a ballistic orbit as a result of dynamical friction. The escape velocity at a given subcluster position is calculated assuming a spherical unperturbed host potential and is indicated by a black dotted curve. All the runs have initially unbound subclusters. Due to the effects of dynamical friction after the phase of core–core interaction the 6:1v2000 subcluster is actually bound to the main system, while subclusters with initial velocities v= 3000 km s⁻¹ and 5000 km s⁻¹ (in the rest frame of the host system) have velocities slightly or much larger than the escape velocity from the main cluster. The gaseous bullet initially follows its dark matter counterpart but before the point of closest approach it is slowed down by ram-pressure. At the same time the morphology of the subcluster gas distribution changes. The contact discontinuity assumes an arc-like shape (partially reducing the effect of ram-pressure) and the bow shock forms. It is interesting to notice that after the point of closest approach the evolution of the relative velocity between gas and dark matter in the subcluster strongly depends on the intensity of the ram-pressure force. In particular for relatively low ram-pressure values (6:1v2000 and 6:1v3000) the gaseous bullet is accelerated towards its dark matter counterpart as soon as it leaves the core of the host cluster where it experienced the largest external densities and ram-pressure. As a result the relative velocity between the gaseous and dark component of the subcluster (whose two components in the orbital plane are represented for 6:1v3000 by the blue solid and dashed curves in Fig. 14) is larger than zero. At the present time the two velocities look comparable in the case of 6:1v3000 while in 6:1v2000, where the ram-pressure acting on the bullet is lower, the acceleration starts earlier and the gaseous bullet is already ∼500 km s⁻¹ faster than its dark counterpart. For larger impact velocities – as in the case of 6:1 – ram-pressure is effective in slowing down the gaseous bullet even at large distances from the centre of the main-cluster. The velocity of the gaseous bullet is therefore smaller than that of its dark counterpart until the subcluster is well outside of the virial radius of the main system.

Figure 14

Run 6:1v3000. Characteristic bullet velocities plotted as a function of time. The black and red solid curves represent the velocity of the dark matter and gaseous component of the subcluster in the plane of sky, while the black dotted line indicates the escape velocity from the main cluster. The blue curves show the relative velocity of the two components of the subcluster in the two directions perpendicular to the line of sight. The red solid and dashed curves show the velocity of the edge of the bullet in the centre of mass and pre-shock gas rest frame, respectively. Finally, the green curves represent the shock velocity obtained by differentiating the shock position. The vertical line indicates the present time. More details are given in the text.

Figure 15

Run 6:1 (top panel) and 6:1v2000 (bottom panel). Same as in Fig. 14.

The green solid curve in each plot represents the velocity of the front shock obtained by differentiating the positions of the shock front at increasing times. The shock velocity rapidly increases after the core–core interaction and at the present time it is always larger than the velocity of the gaseous bullet. In order to calculate the bullet and shock velocities in the system of reference of the pre-shock gas (red and green dashed curve, respectively) the infall velocity of the upstream gas is calculated at different times before and after the present one.

The shock velocity v_S at the present position is 4100, 4500 and 5100 km s⁻¹ in the cases of 6:1v2000, 6:1v3000 and 1:6, respectively, consistently with the shock velocity v_S=Mc_S inferred from the Rankine–Hugoniot jump conditions and with the value provided by observations. Only 6:1v2000 shows a two-dimensional shock velocity well below the observational uncertainties. As found by SF07, after the point of closest interaction the shock velocity is always larger than the velocity of the subcluster mass centroid, but the amount of the difference strongly depends on the model. In particular, in the case of the two low-velocity runs 6:1v2000 and 6:1v3000 v_S is ∼65 and ∼40 per cent larger than the velocity of the dark matter component, while in 6:1 the difference is almost negligible (only 6 per cent).

4.3 Comparison with SF07

Fig. 16 illustrates the three-dimensional velocities of runs 10:1vb0c2 and 10:1vb0c2nfw. The time interval on the horizontal axis is the same as in fig. 4 of SF07. The strength of ram-pressure during the phase of core–core interaction is sensible to the inner gas density profile and the gaseous bullet slows down more efficiently in 10:1vb0c2nfw (bottom panel). However, at the present time the gas velocity is comparable in the two runs, with the gaseous bullet moving slower than the shock front: in 10:1vb0c2 the velocity difference is ∼330 km s⁻¹, in 10:1vb0c2nfw only ∼100 km s⁻¹. The latter result looks comparable with the velocity difference measured by SF07 at time t= 1.36 Gyr, which corresponds to our present position. The infall velocity of the pre-shock gas is 660 km s⁻¹ in 10:1vb0c2 and about 900 km s⁻¹ in the corresponding NFW run. The present-time shock velocity in the system of reference of the pre-shock gas is ∼4270 km s⁻¹ and ∼4370 km s⁻¹ in 10:1vb0c2 and 10:1vb0c2nfw, respectively. These values are both below the best measurement provided by Markevitch (2006) but lay within the observational uncertainties. In particular v_s in 10:1vb0c2nfw is consistent with the result obtained by SF07.

Figure 16

Run 10:1vb0c2 (top panel) and 10:1vb0c2nfw (bottom panel). Same as in Fig. 14.

4.4 Evolution in the L–T diagram

The cluster collision produces a drastic increase in luminosity and temperature. Figs 17 and 18 show the bolometric X-ray luminosity L_Xbol and average spectroscopic-like temperature T_sl as functions of time. Both quantities are calculated for the entire simulated box and scaled to their initial values. The first phase of the interaction – which involves only the external low density regions of the two clusters – is characterized by an identical slow increment of luminosity and temperature in all the adiabatic runs. The jump is associated with the phase of core–core interaction and peaks right after the time of closest approach. The present time is indicated with an empty square and sits on the downturning curve. The high-velocity run 6:1 is associated with the largest increase in temperature (T_sl/T_sl(0) ∼ 3.5) and with the smallest jump in luminosity (L_Xbol/L_Xbol(0) ∼ 4). For decreasing subcluster velocities the peak in temperature becomes slightly less pronounced while the luminosity jump rises by a factor of 1.5. The loss of baryonic material from the subcluster within the high density core of the main system is indeed larger for low-velocity encounters and leads to higher luminosities even if the increment in temperature is smaller with respect to the high-velocity runs. Excluding the bound run 6:1v2000, the amount of gas stripped from the bullet and lying within the virial radius of the main cluster at t= 4 Gyr is 1.6 × 10¹³M_⊙ in 1:6 and 1.8 × 10¹³M_⊙ in 6:1v3000, which in both cases corresponds to almost ∼60 per cent of the initial baryonic content of the subcluster. The difference becomes more pronounced if we consider only the core (r < r_s) of the main cluster, where the mass of subcluster gas is 1.4 × 10¹²M_⊙ in the 6:1v3000 run and one order of magnitude less in the case of 1:6. At 4 Gyr after the beginning of the simulation the total luminosity of the system is similar to the initial one in the case of the bound system 6:1v2000 where the centre of mass of the bullet does not move out to distances beyond the virial radius of the main system. The luminosity drops to 50 per cent or even less of the initial luminosity for runs 6:1v3000 and 6:1, respectively. This decrease in luminosity is motivated by the fact that at the final stages of the simulations a significant fraction of the gas is unbound and very extended. Due to its low density it does not contribute to the luminosity of the system despite the high temperature. In particular, for the same intracluster distance, the high-velocity encounter 6:1 is associated with the highest fraction of unbound material, as will be shown in the next subsection. On the other hand, the final temperature of the system is higher than the one associated with the initially isolated clusters and converges to a value T_sl/T_sl(0) ∼ 1.5 almost independently of the subcluster velocity. The cooling run 3:1clfg shows a somehow different behaviour: both luminosity and temperature profiles have an initial decrement due to the cooling of the central regions of the two approaching clusters. Already during the early phases of the interaction the cooling run moves out of the equilibrium. The peaks in luminosity and temperature are much smaller than the corresponding adiabatic ones (not illustrated). The final luminosity approaches zero.

Figure 17

Bolometric X-ray luminosity L_Xbol as function of time for four selected runs. Luminosity is scaled to its initial value and calculated for the entire simulated volume.

Figure 18

Average spectroscopic-like temperature T_sl as function of time for four selected runs. Temperature is scaled to its initial value and averaged for the entire simulated volume.

Markevitch (2006) found that the Bullet Cluster lies exactly on the L_X–T relation for nearby clusters (Markevitch 1998) but its temperature is much higher than the one expected according to weak lensing mass estimates. Fig. 19 illustrates the drift of the simulated systems along the L_X–T diagram. The position of 1E0657-560 is indicated by a star. L_Xbol and T_sl are calculated in a cylindrical region centred on the centre of mass of the main cluster and with radius equal to its initial virial radius. The luminosity of each model is normalized in such a way that the initial main cluster lies on the L_X–T relation for local clusters despite of the different initial gas fractions. The starting time of the simulations is indicated with a black solid square. All the adiabatic runs present a similar evolution and move roughly parallel to the L_X–T relation, as previously noticed by Rowley, Thomas & Kay (2004) for major mergers in cosmological simulations. During the early stages of the encounter the cluster moves along a curve which is flatter than the observed L_X–T relation: the compression of the low density gas at the outskirts of the cluster produces an increase in temperature which is only marginally accompanied by a luminosity growth. The time when the core of the subcluster enters the virial radius of the main system represents an inversion point in the diagram: despite the formation of a bow-shock the temperature decreases due to the expansion of the main-cluster gas and the presence of low-temperature baryons belonging to the subcluster within the virial radius of the main system. At the same time the luminosity rises as a result of the shock and the cluster moves perpendicularly towards the L_X–T relation. This phase is actually quite short (on average ∼0.12 Gyr) but it characterizes all the adiabatic runs. Later on the cluster moves almost parallel to the L_X–T relation toward larger values of T_sl, with the peak in temperature being reached at the point of closest approach. Most of the runs show a small delay (∼40 Myr) between the time characterized by the highest temperature and the time with highest luminosity, with the curve making a knot in the diagram. The branch of the curve associated with the post core–core interaction is parallel to the increasing one but shifted to smaller luminosities: during the strongest phase of the interaction some of the hot material is lost beyond the virial radius of the main cluster and indeed the largest shift in luminosity is observable in the high-velocity run 6:1. Both luminosity and temperature now decrease until they reach a second inversion point in the curve (in the case of 6:1 it is only a change in slope), associated with the egress of the bullet from the virial radius of the main system. The cooling run 3:1clfg is characterized by a first decrease in temperature which corresponds to the initial phase of thermal instability and central cooling of the main cluster. Later on it moves in the L_X–T diagram similarly to the adiabatic runs although the peaks in luminosity and temperature are much less prominent. Nevertheless, in a pure cooling model the main cluster does not return to a state of thermal equilibrium at the end of the interaction since nothing prevents the central regions from cooling and the system moves toward extremely low values of luminosity and temperature.

Figure 19

Evolution of the interacting system along the L_X–T diagram. Here L_Xbol and T_sl are the bolometric X-ray luminosity and the spectroscopic-like temperature within the virial radius of the main system. The dotted line represents the L_X–T relation by Markevitch (1998) for local clusters. The initial and present time of each simulation, as well as the observed position of the Bullet Cluster on the L_X–T relation (Markevitch 2006) are indicated.

4.5 Density profiles

Fig. 20 illustrates the final structure of the remnants. In the case of the two highest velocities the final time is chosen in such a way that the distance between the centres of the two systems is about 1000 kpc larger than the sum of the virial radii at Time= 0. This occurs at Time= 2.2 and 4 Gyr in the case of 6:1 and 6:1v3000, respectively. The low-velocity run 6:1v2000 is analyzed at the time corresponding to the first apocentre when the core of the subcluster is close to the virial radius of the main cluster. The collisionless component of the main cluster (top left-hand panel) in all cases is not substantially affected by the interaction while the subcluster dark matter halo appears to be strongly perturbed and retains the original spherical symmetry only within its scale radius. The subcluster central density profile (top right-hand panel) is vertically shifted downward without a significant change of slope while for r > 0.2r_vir the loss of material becomes more significant and in the case of low-velocity encounters (v= 2000 and 3000 km s⁻¹) with mass ratio 6:1 the profile shows a large jump of more than one order of magnitude between 0.3r_vir and 0.5r_vir. The amount of dark matter stripped from the bullet in the unbound runs 6:1v3000 and 6:1 is 67 per cent M_vir and 34 per cent M_vir, respectively. Beyond the scale radius the isodensity contours appear elongated and show a large plateau (see also Fig. 21) associated with tidally stripped material.

Figure 20

Dark matter (upper panels) and gas (bottom panels) density profiles of the main (left-hand panels) and subcluster (right-hand panels). Initial values and profiles at the final time (see text) are shown. Blue curves in the bottom left-hand image refer to gas stripped from the subcluster and lying in the potential of the main system. Radius is scaled to the virial radius r_vir of the dark matter distribution.

Figure 21

X-ray surface brightness maps showing the late phases of the evolution of run 6:1v3000. Two successive times are considered, when the centres of the projected total mass distribution (white contours) are ∼2500 kpc (top panel, t= 2.4 Gyr) and ∼5500 kpc (bottom panel, t= 4 Gyr) apart. Logarithmic colour scaling is indicated by the key to the bottom of the figure, with violet corresponding to 0.9 keV and white to 86 keV.

The evolution of gas density profiles is represented by the two panels on the bottom row of Fig. 20. In the 6:1 runs the interaction affects the central slope of the main cluster which becomes shallower while the 3:1 encounter produces the largest deviation from the initial values, with the final density profile shifted down by ∼25 per cent. As mentioned earlier, a not negligible part of the gas within the virial radius of the main cluster originally belonged to the subcluster and was subsequently stripped by ram-pressure during the central phases of the interaction. The density profiles of the stripped gas are drawn in blue for the different runs. In the case of 6:1v2000 only the gas outside the virial radius of the subcluster is considered. Both 6:1v3000 and 6:1v2000 are characterized by a large fraction of subcluster gas lost to the core of the main system, with a flat (6:1v2000) or even positive (6:1v3000) central slope, while the high-velocity run 6:1, despite the larger ram-pressure values, has less time to deposit gas in the central regions and shows a clear cut-off for r < 0.1r_vir. At larger radii, baryonic material is still accreting on to the remnant. Indeed, part of the main cluster gas has been pushed out by the bow shock and is now falling back into the cluster potential together with a fraction of the material lost by the subcluster. Curves in the bottom right-hand panel of Fig. 20 represent the subcluster gas density profiles. The low-velocity subclusters 6:1v3000 and 6:1v2000 retain less than one tenth of their initial gas within 0.1r_vir while the density profiles in the external regions drop by a factor ≥5. The baryonic fraction of the unbound subcluster in our favourite configuration 6:1v3000 is f_g= 6 per cent, comparable with the values provided by (McCarthy et al. 2007) for clusters of similar temperature (the unbound bullet has T_spec∼ 2 keV). In general, the encounter flattens the gas density profile of the bullet core and this effect is evident in the cooling simulation as well. The bullet at Time= 2.2 Gyr has lost almost all the gas within r_vir. Half of this baryonic material is now accelerating towards its dark matter counterpart and will be accreted by the subcluster halo at later times.

Fig. 20 compares dark matter and gas density profiles at times when the main and subcluster are close to a state of virial equilibrium. Despite the changes observable in the gas density profiles and the tidal stripping affecting the bullet beyond its core radius, the central slope of the collisionless component seems to be not strongly perturbed by the interaction (Kazantzidis et al. 2004; Kazantzidis, Zentner & Kravtsov 2006). The situation changes comparing the density profiles of the initial systems with those of the main and subcluster at the present time, when the bullet is located at almost one third of the host virial radius. As illustrated in the two upper panels of Fig. 22, the density profile of the interacting dark matter haloes increases in the inner regions and shows a decrement beyond 0.4r_vir. This behaviour is similar for different models and is therefore independent of the orbital details and of the mass ratios, and is observed for both the main and the subcluster. The two-body relaxation time-scale beyond the softening radius is longer than 10⁴ Gyr (Boylan-Kolchin & Ma 2004; Arad & Johansson 2005), thus implying that dark matter density profiles are not affected by numerical relaxation. As noticed in Section 3, fitting the projected masses around the centres of main and subcluster at the present time we would get more concentrated or more massive systems with respect to the real ones. The last two panels of Fig. 22 represent gas density profiles at the present time. The loss of gas in the central regions of the subcluster is expected as a result of ram-pressure stripping. However, interestingly the inner 200 kpc of the main cluster are also completely devoid of gas originally belonging to the main cluster and replenished with baryonic material stripped from the subcluster. At later times the bullet moves toward the outskirts of the main system and the main cluster gas has time to collapse again into the centre of the potential (bottom left-hand panel of Fig. 20).

Figure 22

Dark matter (upper panels) and gas (bottom panels) density profiles of the main (left-hand panels) and subcluster (right-hand panels). Initial values and profiles at the present time are shown. Radius is scaled to the virial radius r_vir of the dark matter distribution.

4.6 Likelihood of the models

According to Hayashi & White (2006) the probability to find a cosmological configuration where the most massive subhalo has a velocity larger than v_sub is well fitted by the function:  

7

where v₂₀₀ is the virial velocity of the main cluster. Adopting for v_sub the present-time velocity of the subcluster dark matter component (in Table 1) we find a probability of 5 × 10⁻¹⁴, 4.3 × 10⁻⁴ and 0.036 for runs 1:6, 1:6v3000 and 1:6v2000 respectively.

5 CONCLUSIONS

We used high resolution N-body/SPH simulations to perform an extensive parameter study of the ‘Bullet Cluster’ system 1E0657-56. Collisions of two NFW haloes with hot isothermal gas components in hydrostatic equilibrium were studied, adopting initial relative velocities of 5000, 3000 and 2000 km s⁻¹ (in the system of reference of the host system), which is 3.8, 2.3 and 1.5 times the main cluster's virial velocity, respectively. We varied masses, mass ratios, impact parameter, baryonic fraction and concentrations. We also ran simulations with orbital and structural parameters similar to those of SF07. Most of the runs are adiabatic. Radiative cooling is included in two cases, adopting a standard cooling function for primordial gas. We analysed the projected properties of the system at a time where the distance between the centres of the mass distribution associated with the main and subcluster is comparable to the value provided by lensing observations. We have referred to this time as to the present time in our simulations. For a selected subsample of runs we then performed a detailed three-dimensional analysis following their past, present and future evolution.

We have derived constraints on the structure of the progenitors and orbit of the interaction. In what follows we list bounds and limits on the different free parameters considered in this work together with their most likely values.

Impact parameter. Most of the main features in the observed X-ray maps are not well reproduced by encounters with b= 0. Indeed, depending on the relative velocity of the bullet, a perfectly head-on encounter either destroys the X-ray peak associated with the main cluster or does not produce a significant displacement between gas and dark matter in the main system. Moreover, in a pure head-on collision, the gaseous core of the main cluster is displaced along the line which connects the centres of the two total mass distributions, contrary to observations. We find that an impact parameter corresponding to the core radius of the main cluster gas distribution (b= 150 kpc in this paper) provides enough ram-pressure to produce a displacement comparable to observations and introduces asymmetries in the main cluster emissivity map, similar to those detected in X-ray.

Mass ratio. An initial cluster mass ratio as large as 8:1 produces an offset in the bullet much larger than that in the main cluster, contrary to what observed. A similar phenomenon is observed in SF07 – with larger absolute cluster masses and even larger mass ratios – when the concentration of the main halo is c= 5 (bottom right-hand panel of their fig. 7). The peak in projected temperature associated with the shock does not differ from the case of a 6:1 encounter with the same infall velocity.

Adiabatic encounters with mass ratio 3:1 do not match the observed X-ray morphology and the size of the projected temperature jump across the shock. Introducing cooling in the simulations, the temperature peak is cooled to a value comparable with the observed one, but the X-ray map is still not consistent with observations.

We find that the morphology of X-ray maps, the height of the projected temperature jump and the displacement of the gas peaks are best simulated by 6:1 encounters. A significant fraction of the X-ray emission next to the centre of the main cluster is associated with gas stripped from the external regions of the bullet. In these runs the main cluster X-ray peak presents two distinct components: a compact strongly emitting region associated with the displaced core of the host system and a more diffuse component spatially coincident with gas stripped from the subcluster during the central phases of the interaction. The relative luminosity of this secondary component increases with decreasing bullet velocities.

The limit of these mass constraints is that we cannot reproduce the cumulative projected mass associated with the subcluster as inferred by lensing observations (Bradač et al. 2006). The subcluster mass profile in runs 6:1 is close to the one shown by SF07, who use a similar subcluster mass. The only way to fit the observational data using NFW profiles would be to assume a mass ratio as small as 2.7:1 (Nusser 2007), which according to our results would not reproduce X-ray data. This remains an open question. Gravitational lensing analysis of the simulations presented in this paper will be part of a follow-up paper (Mastropietro, Macciò & Burkert, in preparation).

Baryonic fraction. The choice of a different gas fraction f_g does not seem to affect significantly the results. We tested f_g= 17 per cent, consistent with WMAP results, and f_g= 12 per cent, closer to the value provided by X-ray observations.

Concentration of the main halo. Assuming NFW profiles for the two interacting systems, lensing observations provide a very low concentration (c < 2) for the main cluster dark matter halo. Such a low value of c is clearly incompatible with the predictions of ΛCDM for haloes with a similar mass (Macciò et al. 2007). We find that NFW models with a cosmologically motivated concentration c= 6 can reproduce quite well the morphology of the X-ray maps and the observed offset between the X-ray peaks and the centres of the projected total mass distribution. A low concentrated (c= 4) main cluster does not survive the collision and its X-ray emissivity peak is destroyed. Runs 10:1vb0c2 and 10:1vb0c2nfw, with initial conditions close to those adopted in SF07, provide further examples of interactions with low concentrated main haloes. Here the mass ratio is very large (10:1) and the peak associated with the main cluster gas survives the core–core encounter due to the relatively lower ram-pressure generated during the interaction. Nevertheless – as observed also by SF07 – the displacement in the bullet is negligible.

We have shown that the apparent concentration of the two interacting haloes, if measured at the present time, when the bullet is still located well within the virial radius of the main cluster, is higher than that measured in a situation of virial equilibrium. Indeed the three-dimensional dark matter density of the main and the subcluster increases in the inner regions and shows a decrement beyond 0.4r_vir, as a consequence of the interaction. This effect goes in the direction of increasing the conflict between lensing data and simulations. The new large field analysis will likely provide us with some new indications (D. Clowe, private communication).

Subcluster velocity. For our favourite mass ratio 6:1 we have tested different subcluster velocities. The projected temperature jump across the shock discontinuity gives important indications about the relative velocity of the interaction. Indeed, both, the height and the thickness of the peak change, with the peak becoming broader for decreasing bullet velocities due to the lower pressure exercised by the pre-shock gas after the central phases of the encounter. Both runs 6:1v3000 and 6:1v2000 match the observed temperature jump quite well, while for larger velocities (6:1) the peak becomes too narrow and pronounced compared to observations. An initial velocity as large as v= 5000 km s⁻¹ (∼4300 km s⁻¹ in the centre of mass system of reference) has to be excluded also due to the shape of the contact discontinuity, which appears to be much broader than observed. The shock Mach number in this run is ⁠, higher than the one provided by observations but still within the observational uncertainties.

Run 6:1v3000, with initial velocity v= 3000 km s⁻¹ (∼2570 km s⁻¹ in the centre of mass system of reference) and present-time dark matter velocity v∼ 3100 km s⁻¹ (again in the centre of mass rest frame), reproduces most of the main X-ray features: besides the projected temperature jump, the peculiar morphology of the X-ray emission associated with the main cluster, the relative surface brightness between the main and the subcluster, the shape of the bow shock and of the contact discontinuity. The Mach number is ⁠, which is the best value provided by Markevitch (2006).

Decreasing the relative velocity (6:1v2000 has a present-time dark matter velocity which is almost 300 km s⁻¹ smaller than that used by Springel & Farrar 2007) the bullet becomes much less bright with respect to the centre of the host system and the morphology of the main cluster X-ray peak is not reproduced.

Likelihood of our favourite model. ln run 6:1v3000 the bullet is initially unbound, with a large infall velocity of 2.3v₂₀₀ in the system of reference of the main cluster. At the end of the interaction its velocity has decreased to a value comparable to the escape velocity but the subcluster remains unbound and leaves its host halo with positive energy. According to Benson (2005) the assumption that the subcluster is initially unbound is not unlikely in a ΛCDM universe. However, due to the effects of dynamical friction the probability that a satellite remains unbound after the interaction with its host halo is quite small, only about 0.3 per cent of all orbits in Benson (2005). This probability is not far from the likelihood of finding such a high-velocity encounter in the Millenium Run simulation as calculated by Hayashi & White (2006) (0.43 per cent for the choice of masses adopted in this paper). Reducing the impact velocity (run 6:1v2000) would increase this probability to 3.6 per cent, but would also reduce the bullet's brightness beyond what is observed. Concerning the morphology of the main cluster X-ray peak, we have seen in Section 3 that it can be affected by a top-down inclination of the plane of the encounter with respect to the line-of-sight velocity. The Bullet Cluster seems to be a very peculiar case in a ΛCDM universe, but it could also represent an interesting issue for cosmological structure formation in a cold dark matter scenario, especially since another interacting cluster (Abel 576) showing similar velocity gradients has been observed (Dupke et al. 2007). Another possibility is that a more complex thermodynamical model might lead to a different X-ray morphology, with more luminous bullets for lower impact velocities. We have seen that a pure cooling model neglecting energetic feedback does not give a realistic description of the interacting system since overcooling in the central regions follows a thermal instability in the early phases of the interaction.

During the interaction the subcluster loses nearly 67 per cent of its initial dark matter mass. The final baryonic fraction of the unbound substructure is f_g∼ 6 per cent, which is quite low but comparable with the values provided by X-ray observations (McCarthy et al. 2007) for clusters of similar temperatures (about 2 keV).

Further observations.

• Projected temperature maps reveal some interesting features. As already found by SF07, the bullet remains relatively cold despite its dominant X-ray emission. For a short time immediately after the central phases of the interaction a high-temperature region appears next to the centre of the main cluster mass distribution. This feature – still visible at the present time – is not spatially coincident with the main cluster gaseous core and appears associated with hot gas stripped from the subcluster within the core of the host system. Its location partially corresponds to the X-ray diffuse peak, while its eastern component could be related to the high-temperature region observed by Markevitch et al. (2002) southeast of the main-cluster peak.

• Due to the drop in ram-pressure after the central phases of the interaction and gravitational acceleration by its dark halo counterpart, at the present time the gaseous bullet is not necessarily slower than the dark matter halo (Milosavljević et al. 2007). We explored this question in greater details and found that the behaviour of the subcluster gas after the point of closest approach strongly depends on the initial velocity of the bullet. For relatively low velocities (and ram-pressure values, run 6:1v3000 and 6:1v2000) the gaseous bullet is accelerated towards its dark matter counterpart as it leaves the core of the main cluster. At the present time the gaseous bullet moves as fast as its dark matter halo in 6:1v3000 and even 500 km s⁻¹ faster in 6:1v2000. For larger encounter velocities (6:1) ram-pressure is more effective in slowing down the bullet even beyond the core radius of host cluster and the gaseous bullet is always slower than dark matter.

• It has been already noticed by Milosavljević et al. (2007) and SF07 that the subcluster velocities do not coincide with the shock velocity v_s as measured by observers. In the case of the three runs (6:1v3000, 6:1v2000 and 6:1clfg) which best reproduce the projected temperature jump, the Mach number and shock velocity determined using the Rankine–Hugoniot conditions across the shock discontinuity show a good agreement with the values provided by Markevitch (2006). However, if the deprojected temperature of the pre-shock gas is used, v_s is higher than observed (5100 km s⁻¹ in the case of 6:1v3000). The two-dimensional shock velocity calculated by tracking the shock position as function of time has values between ∼4000 and 5000 km s⁻¹ in our 6:1 runs. We find that, although after the point of closest approach v_s is always larger than the velocity of the subcluster dark matter halo, the difference depends on the model and becomes less significant for higher velocities. In particular in the high-velocity run 6:1 the present time shock velocity is only 6 per cent higher than the dark bullet.

• The collision produces a drastic increase in luminosity and temperature. The highest velocity impacts are associated with the largest increase in temperature (and shocks) and the smallest peaks in luminosity due to the fact that in high-velocity encounters a significant fraction of the baryons stripped from the bullet are deposited at large radii within the main cluster. We followed the evolution of the main cluster in the L_X–T diagram during the different phases of the interaction and find that it moves roughly parallel to the L_X–T relation for nearby clusters. The maximum temperature is reached at the point of closest approach and the peak in luminosity immediately afterwards. The branch of the curve associated with the post core–core interaction is still parallel to the L_X–T relation but shifted towards smaller luminosities with respect to the early increasing branch, due to the loss of gas beyond the virial radius of the main cluster.

• After the encounter, as soon as the bullet is close to or beyond the virial radius of its host system, the dark matter density profile of the main cluster does not deviate any more significantly from the original one, while the gas profile becomes shallower in the central regions. The situation changes drastically if we compare the density profiles of the initial systems with those at the present time, when the bullet is still located well within the virial radius of the main cluster. In this case the interacting systems are not in virial equilibrium and the dark matter densities of both the main and subcluster increase in the inner regions and show a decrement beyond 0.4 r_vir. Comparing these present time profiles with NFW haloes would give wrong estimates of the halo parameters and suggest higher concentrations. At the present time the centre of the host system is completely devoid of main cluster gas – as a result of the displacement of the main cluster gaseous core by the bullet – while the baryons stripped from the bullet have replenished the central regions of the host. This is only a temporary situation, however, since a few 10⁹ years later the main cluster gas falls back and again becomes the dominant component in the central regions, although a significant fraction of gas stripped from the bullet is still present within the host system. The final subcluster dark matter density profiles seem significantly affected by the interaction beyond their scale radius, where the isodensity contours are elongated and show a plateau.

We would like to thank R. Piffaretti, S. Borgani, A. Biviano, T. Naab, A. Maccio' and M. Girardi for useful discussions. We acknowledge T. Quinn for support with the TIPSY X-ray package, D. Clowe for providing the upper panel of Fig. 1 and S. Ameglio for the Raymond-Smith code. The work was partly supported by the DFG Sonderforschungsbereich 375 ‘Astro-Teilchenphysik’. The numerical simulations were performed on a local SGI-Altix 3700 Bx2 which was partly funded by the cluster of excellence ‘Origin and Structure of the Universe’.