Multi-phenomena Modeling of the New Bullet-like Cluster ZwCl 008.8+52 Using N-body/Hydrodynamical Simulations

We modeled ZwCl008 in 3D using the Eulerian N-body/hydrodynamic code FLASH (developed at the Center for Astrophysical Thermonuclear Flashes at the University of Chicago; Fryxell et al. 2000; Ricker 2008). FLASH is a publicly available AMR code, which can be run in parallel computer architectures. We assumed a binary merger for ZwCl008 and included dark matter and gas self-consistently, taking their gravity into account dynamically. We used our well-established method to carry out merging cluster simulations (Molnar et al. 2012, 2013a, 2013b; Molnar & Broadhurst 2015, 2017). For our simulations, we adopted a large box size (13.3 Mpc on a side) to capture the outgoing merger shocks and avoid losing mass during the time we ran our simulations. Our highest resolution, 12.7 kpc, was reached at the cluster centers, merger shocks, and in the turbulent regions behind the shocks. We chose a 3D Cartesian coordinate system, x, y, z, with the x, y plane containing the centers of the clusters and the initial (relative) velocity vector of the infalling cluster in the positive x direction. We included shock heating, the most important non-adiabatic process in merging clusters, and ignored other heating and cooling processes (e.g., supernova and AGN feedback, radiative cooling, etc.).

The initial models of the clusters were assumed to have spherical geometry with cutoffs of the dark matter and gas densities at the virial radius, R_vir. We assumed an NFW model (Navarro et al. 1997) for the dark matter distribution,

$$\begin{eqnarray}&&{\rho }_{\mathrm{DM}}(r)=\displaystyle \frac{{\rho }_{s}}{x(r){(1+x(r))}^{2}};\,r\leqslant {R}_{\mathrm{vir}},\end{eqnarray} \tag{ 1 }$$

where x(r) = r/r_s (r_s = R_vir/c_vir) and ρ_s are the scaling parameters for the radius and the amplitude of the density, respectively, and c_vir is the concentration parameter. We adopted a non-isothermal β model for the gas density,

$$\begin{eqnarray}&&{\rho }_{\mathrm{gas}}(r)=\displaystyle \frac{{\rho }_{0}}{{(1+{y}^{2})}^{3\beta /2}};\,r\leqslant {R}_{\mathrm{vir}},\end{eqnarray} \tag{ 2 }$$

where y = r/r_core is the scaling parameters for the radius r (in units of the core radius, r_core), and ρ₀ is the density at the center of the cluster. The exponent β determines the falloff of the density at large radii. We adopted β = 1, suggested by cosmological numerical simulations for the large-scale distribution of the intracluster gas in equilibrium (excluding the filaments; see Molnar et al. 2010).

We derived the gas temperature as a function of the radius, T(r), assuming hydrostatic equilibrium, adopting γ = 5/3 for the ideal gas equation of state. We used a gas fraction f_gas = 0.14 and represented baryons in galaxies together with collisionless dark matter particles, since, for our purposes, galaxies can also be considered collisionless.

It is less straightforward to model a stable dark matter density distribution than a gas distribution. The hydrostatic equilibrium assumption provides a stable distribution for the gas, but the dark matter, modeled as particles, has no pressure, interacting only gravitationally with itself and the gas, which means that they move on orbits in the potential of the cluster. We use the local Maxwellian approximation for the amplitude of the velocities of the dark matter particles: we randomly sample a Maxwellian distribution with a velocity dispersion as a function of radius, r, σ_v (r) derived from the Jeans equation assuming that the distribution of σ_v (r) is isotropic (Łokas & Mamon 2001). We assumed an isotropic distribution for the direction of the velocity vectors (for more details of the setup for our simulations see Molnar et al. 2012).

We have run a series of FLASH simulations varying the initial masses, concentration parameters, impact parameter, and infall velocity of our models. Our aim was to find a physical model for ZwCl008 with a reasonable agreement with observations.

The long bright radio relic to the east most likely marks the location of the back shock (as in the Sausage cluster, as demonstrated by Molnar & Broadhurst 2017), due to the limb brightening of the spherical surface of the relic viewed in projection. The position of the bow shock associated with the infalling cluster is much less certain. The radio relics do not show a bright edge here, as in the case of the back shock, and this forward shock is not detected by the X-ray observations (Golovich et al. 2017). We estimate from our models that it should lie ahead of the small radio relic that we nevertheless associate with this shock. The reason for this displacement can be simply geometric in origin. We may approximate the surface of the bow shock with a hyperboloid of revolution soon after the first core passage (see the middle panel in Figure 1). If there are relics covering the shock surface well and the symmetry axis of the hyperboloid of revolution is not rotated far from the plane of the sky (say less than 40°), the relics will delineate the shock surface, showing a bright rim in projection because the radio emission, which is optically thin, adds up in projection. Indeed, notable examples of large anomalously bright shock fronts, in particular, the bow shock in the Sausage cluster (van Weeren et al. 2011a), are known, which modeling has shown to coincide with their observed projected shock front (Molnar & Broadhurst 2017). However, when there are only a few relics tracing the shock surface, as in the case of the bow shock in ZwCl008, the position of the shock is not well defined because we do not see a bright edge but only a few randomly scattered relics. When the rotation angle is reasonably small, no patches can fall outside of the edge of the hyperboloid of revolution in projection (as it surrounds a convex volume), thus we expect the projected shock surface to be outside of the relics, ahead of them, toward the direction of the assumed motion of the shock. This can be clearly seen in Figure 4 of van Weeren et al. (2011a), where, in projection, the simulated radio emission from relics is inside the projected shocks.

To find the best model, we used the following set of criteria: the position of the back shock as suggested by the radio relics (van Weeren et al. 2011b), the distance and positions of the mass centers of the subclusters derived from gravitational lensing (Golovich et al. 2017), the positions of the bright and secondary X-ray peaks associated with the infalling and main cluster, and the X-ray morphology (e.g., the cometary shape of the bright X-ray peak; the secondary peak, which is much fainter and forms a plateau; the two peaks, which are not well separated; Golovich et al. 2017).

A systematic parameter search is currently beyond reach of conventional high-speed computing resources based on CPUs. Therefore, first we analyzed simulations in parameters randomly distributed in a larger parameter space: initial cluster masses ${10}^{14}\,{M}_{\odot }$ ≤ M_vir1,2 ≤ 1.2 × ${10}^{15}\,{M}_{\odot }$, concentration parameters 3 ≤ c_vir1,2 ≤ 10, and infall velocities and impact parameters of 500 kpc ≤ V_in ≤ 3000 kpc and 0 kpc ≤ P ≤ 600 kpc. We choose simulations from our extended suite of simulations we carried out earlier, or performed new simulations. This process narrowed down the parameter space substantially. We also carried out a suite of simulations to provide a rough estimate of the errors on the infall velocity, impact parameter, and masses of the merging system before the collision.

Table 1 contains a list of the initial parameters of those simulations we discuss in this paper. This is a narrow subset of parameter space that contains our best solution with enough spread to illustrate the observable effects of moving away from it. The first column contains the IDs of our runs as "PijkVmn," where "Pijk" is the impact parameter in kpc, and "Vmn" is the infall velocity in units of 100 km s⁻¹. In columns 2–5, we show the virial masses, M_vir (in units of ${10}^{14}\,{M}_{\odot }$), and concentration parameters (c_vir) of the two subclusters. Columns 6 and 7 contain the impact parameters (P) and infall velocities (V_in) in km s⁻¹. We define the infall velocity as the relative velocity of the two clusters at the time of their first contact, i.e., when the distance between them is R_vir1 + R_vir2.

Table 1.  IDs and Input Parameters for the Different Models Used in Our Hydrodynamical Simulations

IDa	Mvir1 b	cvir1b	Mvir2 c	cvir2c	Pd	Vine
P400V18B	7.0	5	5.0	5	400	1800
P050V05G	5.7	4	1.2	4.6	50	500
P200V10G	5.7	4	1.2	4.6	200	1000
P400V16G	5.7	4	1.2	4.6	400	1600
P300V18	7.0	5	5.0	5	300	1800
P500V18	7.0	5	5.0	5	500	1800
P400V15	7.0	5	5.0	5	400	1500
P400V20	7.0	5	5.0	5	400	2000
P400C44	7.0	4	5.0	4	400	1800
P400C66	7.0	6	5.0	6	400	1800
P400C56	7.0	5	5.0	6	400	1800
P400C54	7.0	5	5.0	4	400	1800
P400C65	7.0	6	5.0	5	400	1800
P400C64	7.0	6	5.0	4	400	1800
P400C45	7.0	4	5.0	5	400	1800
P400C46	7.0	4	5.0	6	400	1800
P400C68	7.0	4	5.0	6	400	1800
P400V18M1	7.0	5	5.5	5	400	1800
P400V18M2	7.0	5	4.5	5	400	1800
P400V18M3	7.5	5	5.5	5	400	1800
P400V18M4	7.5	5	5.0	5	400	1800
P400V18M5	7.5	5	4.5	5	400	1800
P400V18M6	6.5	5	5.5	5	400	1800
P400V18M7	6.5	5	5.0	5	400	1800
P400V18M8	6.5	5	4.5	5	400	1800
P400V18C68	7.0	6	5.0	8	400	1800

Notes.

^aIDs of the runs indicate the impact parameters in kpc and the infalling velocities in units of 100 km s⁻¹. ^bVirial mass in units of ${10}^{14}\,{M}_{\odot }$ and concentration parameter for the main cluster (1). ^cVirial mass in units of ${10}^{14}\,{M}_{\odot }$ and concentration parameter for the infalling cluster (2). ^dImpact parameter in units of kpc. ^eInfall velocity of cluster 2 in km s⁻¹.

Download table as:  ASCIITypeset image

The second panel in Figure 1 shows a simulated X-ray color image of our best model at the epoch of t_b = 428 Myr after the first core passage for ZwCl008, run P400V18B, with infall velocity V_in = 1800 km s⁻¹, impact parameter P = 400 kpc, masses M_vir1 = 7 × ${10}^{14}\,{M}_{\odot }$ and M_vir2 = 5 × ${10}^{14}\,{M}_{\odot }$, and concentration parameters c_vir1 = 5 and c_vir2 = 5 (main and infalling cluster). The concentration parameters for the two colliding clusters we derived are close to c_virD,1 = 5.25 and c_virD,2 = 5.4, as suggested by the mass–concentration relation of Duffy et al. (2008) for virial overdensities. The yellow contours represent the projected dark matter distribution from our simulation. The white contours are based on radio observations; the white dashed circle and the horizontal bar represent R₅₀₀ and a physical length of 1 Mpc (from Golovich et al. 2017; as in the first panel). The viewing (rotation) angles were chosen the following way. First, we rotated the system by an angle φ = 29° ("roll angle") around the axis connecting the two dark matter centers (rotation around the y-axis) then we rotated this axis by an angle θ = − 70° out of the x–y plane. θ was chosen to provide a projected distance between the two dark matter centers, D = 940 kpc, to match the positions of the observed centers based on weak lensing mass reconstruction; the roll angle was chosen to find the best match with the observed X-ray morphology (Golovich et al. 2017). We chose the output (epoch) that could be rotated in such a way that the position of the back shock is near the eastern edge of the observed long radio relic to the east of the X-ray peak and the bow shock in the west is not inside the radio relic associated with it, as radio observations suggest (van Weeren et al. 2011b). The shocks in our simulated images were located based on projected pressure gradients. A detailed description of our method to generate mock X-ray and mass surface density images can be found in Molnar & Broadhurst (2015).

In Figure 2, we show images of models for ZwCl008, which do not satisfy our requirements for a good match with the data. In this figure, the color images represent the simulated X-ray surface brightness, and the color code for the contours is the same as in the second panel in Figure 1. The procedure is as follows. First, we aligned the position of the X-ray peak associated with the infalling cluster with that of the observed, then we chose the projection angles to match the distances and position angles between the dark matter centers to match those derived from weak lensing observations. In the third step, if it was possible, we selected cases where the position of the back shock is near the eastern edge of the long radio relic to the east of the X-ray peak.

Zoom In Zoom Out Reset image size

In each panel in the first row, the positions of the X-ray peak associated with the infalling cluster and the back shock match the observations. The first and second panels in the first row show runs with only the impact parameters changed to P = 300 and 500 kpc (runs P300V18 and P500V18) from the best model (run P400V18B with P = 400 kpc). The X-ray morphology of P300V18 (first panel) is elongated along the line connecting the two dark matter centers, but it has only one X-ray peak, no secondary peak. Run P500V18 displays two X-ray peaks with a large separation and very small offset from the dark matter centers, as well as the brightest X-ray peak associated with the main cluster. The third and fourth panels show simulations changing only the relative (infall) velocity of the infalling cluster to V_in = 1500 and 2000 km s⁻¹ (runs P400V15 and P400V20) to bracket our best model with V_in = 1800 km s⁻¹. The X-ray morphology in the third panel (V_in = 1500) is similar to that of the first panel, the brightest X-ray peak is associated with the main cluster, not to the infalling cluster as observed, and the fourth panel (V_in = 2000 km s⁻¹) show two bright X-ray peaks.

The first and second panels in the second row display simulations changing only the concentration parameters to c_vir1 = 4 and c_vir2 = 4 (run P400C44), and the masses of the colliding clusters to M_vir1 = 6.5 × ${10}^{14}\,{M}_{\odot }$ and M_vir2 =4.5 × ${10}^{14}\,{M}_{\odot }$ (run P400M8). The X-ray emission in these two panels show two bright X-ray peaks, again differing significantly with respect to the observed morphology. In the first panel, the X-ray peaks are well separated; in the second panel the main X-ray peak is associated with the main cluster and not with the infalling cluster. The third and fourth panels in the second row display our best model (run P40V018) but before and after the best epoch, at t₁ = 412 Myr and t₂ = 444 Myr after the first core passage. The best epoch is t_obs = 428 Myr. The third panel in the second row shows only one X-ray peak and the back shock is closer to the peak than observed, the fourth panel shows an epoch when the back shock is too far from the X-ray peak, thus its position does not match that of the observed. See Table 1 for a list of initial parameters for the runs. Figure 2 demonstrates that changing the virial masses of the subclusters by ΔM = 5  × ${10}^{13}\,{M}_{\odot }$, the concentration parameters by Δc_vir = 1, the impact parameter by ΔP = 100 kpc, and the infall velocity by ΔV_in = 300 km s⁻¹ changes the X-ray morphology substantially, and the colliding systems do not resemble ZwCl008, giving a rough estimate of the errors in our derived parameters.

Sampling a larger parameter space initially, we found that values of the initial cluster masses, concentration parameters, infall velocities and impact parameters very different from those of our best model provide very different shock positions and X-ray morphology. The positions of the mass centers are less constraining. It is unlikely that, using the parameters we varied, we missed the "best model." This can be proved only with a densely sampled parameter space, which is beyond the scope of our paper.

We conclude from our suite of N-body/hydrodynamical simulations that ZwCl008 is viewed at about t_obs = 428 Myr after first core passage, and that the infalling cluster has a mass M_vir2 ∼ 5 $\times {10}^{14}\,{M}_{\odot }$ moving to the west with an impact parameter P ∼ 400 kpc and infall velocity of V_in ∼ 1800 km s⁻¹, disrupting the gas of the main cluster with mass M_vir1 ∼ 7 $\times {10}^{14}\,{M}_{\odot }$, so that it lies to the west of the main cluster at the observed epoch. Our best model also clearly confirms that the disrupted gas of the main cluster is offset from its dark matter center, as the observations seem to indicate. Our simulations clearly demonstrate that the merging cluster, ZwCl008, is in the outgoing phase just after the first core passage, before the first turnover. The gas and dark matter associated with the two components are moving outward (the infalling cluster moving to the east, the main cluster to the west; see fist and second panels in Figure 1).

These results are in broad qualitative agreement with the results of Golovich et al. (2017). However, our self-consistent simulations break the degeneracy between the outgoing phase after the first core passage and the infalling phase after the first turnover found by Golovich et al. and clearly prefer an epoch close to post-collision and exclude later ones, since by that time the shock fronts will have long left the system. Our best model has a somewhat larger total mass, 1.2 × 10¹⁵ M_⊙, than the one derived by Golovich et al. using gravitational lensing, ${M}_{\mathrm{virG},\mathrm{total}}={9.8}_{-2.6}^{+4.4}\,$ $\times {10}^{14}\,{M}_{\odot }$ (note we converted their M_{200 G,total} to virial mass based on the scaling of Hu & Kravtsov 2003); however, these masses agree within 1σ. Golovich et al. (2017) also derived total mass estimates based on velocity dispersion measurements and X-ray temperature and luminosity scaling relations, which were smaller but within 1σ of the estimate based on lensing (M_200,total = 6, 6.5, and $4.66\times \,{10}^{14}\,{M}_{\odot }$ respectively, versus ${\text{}}{M}_{G200,\mathrm{total}}\,=\,{8.0}_{-2.1}^{+3.6}\,$ $\,\times {10}^{14}\,{M}_{\odot }$; see their Table 5). In their dynamical model, Golovich et al. adopted the masses they derived based on gravitational lensing as the more reliable value. However, our best model has a pericenter velocity of V_perictr = 3515 km s⁻¹, which is much larger than the V_perictr,G = 1800 km s⁻¹ derived by Golovich et al. (2017). This is partly due to the fact that our best model has a larger total mass than the one used by Golovich et al.

We performed a simulation, run P05V05G, with initial conditions suggested by Golovich et al. (2017). See Table 1 for the initial parameters of this run. We used V_in = 500 km s⁻¹ to get the velocity at the pericenter within 1σ of the velocity quoted in Golovich et al. (${V}_{\mathrm{perictr},G}={1800}_{-300}^{+400}\,$km s⁻¹). Note that in the case of a small impact parameter, the infalling cluster reaches the deep potential well of the main cluster, acquiring higher velocities. We show the result in the third panel in Figure 1. There are two X-ray peaks, but the main X-ray peak is associated with the gas of the main cluster. The infalling cluster has a very faint X-ray peak, contrary to observations. The back shock has already left the cluster gas. We chose a projection, which provides a rough match to the positions of the two X-ray and dark matter peaks. Most of the gas of the infalling cluster has been stripped off as a result of the low mass of the infalling cluster and the relatively high velocity at pericenter (∼2000 km s⁻¹; note that in the case of lower pericenter velocity, we found that the infalling cluster does not get as far from the main cluster as the gravitational lensing observations suggest). Golovich et al. (2017) assumed zero impact parameter, which would make it even easier to strip the gas from the infalling cluster. We chose a non-zero, but small, impact parameter since the observed X-ray morphology is not symmetric. We carried out simulations using the Golovich models but with different infall velocities and impact parameters (runs P200V10G and P400V16G) to check if using these models provide a better match with observations, but we found that they also do not resemble the observed features of ZwCl008.

X-ray observations measure a gas temperature weighted with emission in the LOS convolved with the response of the detectors, which we approximate with the spectroscopic temperature (Mazzotta et al. 2004). In Figure 3, we show the spectroscopic temperature, T_spec, of our best model (see, e.g., Molnar & Broadhurst 2017 for our method to derive T_spec from simulations) with contours of the X-ray emission overlaid. The X-ray peak is close to the center in a cooler region where the gas density is higher; the shocks are not visible in the observations (see the left panel of Figure 1). We obtain 6 keV for the X-ray emission dominating the high-density region (see Figure 3), which is somewhat higher than the observed global temperature, 4.9 ± 0.13 keV (Golovich et al. 2017).

Zoom In Zoom Out Reset image size

We derive the Mach numbers for the merging shocks based on our best model using the temperature jump at the shocks, as it is done using temperature based on X-ray observations. The Rankine–Hugoniot jump conditions provide the connection between the temperature jump, T₂/T₁, and the Mach number,

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{2}}{{T}_{1}}=\displaystyle \frac{5{{ \mathcal M }}^{2}+14{{ \mathcal M }}^{2}-3}{16{{ \mathcal M }}^{2}},\end{eqnarray} \tag{ 3 }$$

where T₁ and T₂ are the pre- and post-shock temperatures (e.g., Akamatsu et al. 2015; Molnar & Broadhurst 2017; for a review, see Markevitch & Vikhlinin 2007). Using the temperature jump from our simulations at the shocks in Equation (3), we obtain Mach numbers for the merging shocks directly from the physical (not the observed) temperature ratios (the bow shock and back shock to the west and east): ${{ \mathcal M }}_{w,\mathrm{simu}}=5.5$ and ${{ \mathcal M }}_{e,\mathrm{simu}}=6.5$.

Recently, Mach numbers significantly larger than X-ray measurements of some merging clusters were derived from radio observations (e.g., the Sausage cluster: ${{ \mathcal M }}_{\mathrm{radio}}=4.58\pm 1.09$, ${{ \mathcal M }}_{{\rm{X}} \mbox{-} \mathrm{ray}}={2.7}_{0.4}^{+0.7}$, van Weeren et al. 2010, Stroe et al. 2013, and Akamatsu et al. 2015; RX J0603.3+4214, the "Toothbrush" cluster: ${{ \mathcal M }}_{\mathrm{radio}}=2.8$, ${{ \mathcal M }}_{{\rm{X}} \mbox{-} \mathrm{ray}}=1.2$, van Weeren et al. 2016; A2255: ${{ \mathcal M }}_{\mathrm{radio}}=2.7$, ${{ \mathcal M }}_{{\rm{X}} \mbox{-} \mathrm{ray}}=1.4$, Akamatsu et al. 2017; A3667: ${{ \mathcal M }}_{\mathrm{radio}}=4.6$, ${{ \mathcal M }}_{{\rm{X}} \mbox{-} \mathrm{ray}}=2.09$, Sarazin et al. 2014).

Using our N-body/hydrodynamical simulations, we demonstrated that gas temperatures around shocks in merging clusters derived from X-ray observations may be much lower than physical temperatures depending on the viewing angle of the collision (Molnar & Broadhurst 2017). In Figure 4, we show the physical temperature and T_spec (solid and dashed lines) along a cut through the center of the X-ray peak based on our best model for ZwCl008. Using Equation (3) with the values of T_spec at the shocks, we estimate the Mach numbers that would be derived from X-ray observations: ${{ \mathcal M }}_{w,{\rm{X}} \mbox{-} \mathrm{ray}}=2.8$ and ${{ \mathcal M }}_{e,{\rm{X}} \mbox{-} \mathrm{ray}}=2.2$. Interestingly, these Mach numbers for the shocks in ZwCl008 we derived assuming X-ray observations based on our N-body/hydrodynamical simulations using spectroscopic temperatures are close to those derived from radio observations assuming a linear diffusive shock acceleration (DSA) mechanism, ${ \mathcal M }=2.2$–2.4 (van Weeren et al. 2011b; Kierdorf et al. 2017).

We have found a similar discrepancy modeling the northern shock in the Sausage cluster using N-body/hydrodynamical simulations (Molnar & Broadhurst 2017). The Mach numbers we derived for this shock based on our simulations using physical and spectroscopic temperatures were ${{ \mathcal M }}_{\mathrm{phys}}=6.5$ and ${{ \mathcal M }}_{{\rm{X}} \mbox{-} \mathrm{ray}}=4.5$. Radio and X-ray observations of this shock suggested ${{ \mathcal M }}_{\mathrm{radio}}=4.58\pm 1.09$ (assuming a linear DSA mechanism) and ${{ \mathcal M }}_{{\rm{X}} \mbox{-} \mathrm{ray}}={2.7}_{0.4}^{+0.7}$ (van Weeren et al. 2010; Stroe et al. 2013; Akamatsu et al. 2015). The ${{ \mathcal M }}_{\mathrm{phys}}$ we derived for this shock is within 1.5σ of the Mach number derived from these radio observations.

Our hydrodynamical simulations of merging clusters (ZwCl008, the Sausage cluster) suggest that Mach numbers of merging shocks based on X-ray observations can be significantly underestimated due to projection effects. Other explanations for this Mach number discrepancy are based on using more complicated particle acceleration models than the linear DSA. Analyzing the radio observations of the Toothbrush cluster, van Weeren et al. (2016) suggested that either the shock has a complex morphology in the LOS or the radio emission is due to a reaccelerated population of electrons. Van Weeren et al. conclude that more theoretical work is necessary to understand the acceleration and reacceleration processes in merger shocks.

Stroe et al. (2014), analyzing radio observations of the northern relic in the Sausage cluster, used different spectral aging models to address the discrepancy between Mach numbers derived from radio and X-ray observations of the shock. For the northern shock, they derived ${{ \mathcal M }}_{\mathrm{radio}}={2.90}_{-0.13}^{+0.10}$, which is consistent with the Mach number derived from X-ray observations. However, as they noted, their solution has an aging problem: the aging of the accelerated high-energy electrons derived from the shock speed is not compatible with the aging derived from the radio spectrum data. Stroe et al. concluded that the models they investigated have difficulties in explaining all of the observations. Recently, Hoang et al. (2017) carried out deep radio observations of the Sausage cluster. Assuming linear DSA, they derived ${{ \mathcal M }}_{\mathrm{radio}}={2.7}_{-0.3}^{+0.6}$ for the northern shock. They concluded that high acceleration efficiency is necessary to explain the radio observations, which is problematic for the DSA mechanism.

However, using time-dependent DSA simulations, Kang et al. (2012, 2017) showed that both larger and smaller Mach numbers derived from radio and X-ray observations can be explained. Kang et al. (2012) carried out time-dependent DSA simulations and found that models with and without preexisting cosmic-ray electrons can explain the radio observations of the northern shock in the Sausage cluster providing Mach numbers ${ \mathcal M }\approx 2$ and ${ \mathcal M }\approx 4.5$, respectively. Recently, Kang et al. (2017) performed time-dependent DSA simulations including radiative cooling and post-shock turbulent acceleration to explain the different Mach numbers derived for the shock in the Toothbrush cluster based on radio and X-ray observations. Kang et al. derived synchrotron emission from the accelerated electrons based on a model for low Mach numbers, ${{ \mathcal M }}_{\mathrm{radio}}\leqslant 2$, which includes a preexisting population of cosmic-ray electrons assuming a flat energy spectrum, and a model for larger Mach numbers, ${{ \mathcal M }}_{\mathrm{radio}}\approx 3$, which assumes either shock-generated suprathermal electrons or preexisting cosmic-ray electrons. They found that both models can explain reasonably well the observed integrated radio spectrum of the northern relic in the Toothbrush cluster. However, they prefer their model for the higher Mach number, because explaining the lower Mach number with a uniform population of preexisting cosmic-ray electrons with a high cutoff energy (≥40 GeV) on a large scale of 400 kpc would be difficult, due to the short cooling time of those electrons.

The properties of merging shocks were studied in detail by making use of N-body/hydrodynamical simulations (e.g., Ha et al. 2018; Molnar & Broadhurst 2017). In this section, we study the evolution of merging shocks as a function of time around the first core passage before the first turnover using our best solution for ZwCl008 as an example.

Zoom In Zoom Out Reset image size

The evolution of merging shocks is illustrated in Figure 5. The panels in this figure show a color image of the X-ray emission with contours of the dark matter distribution (yellow) and the merging shocks (green) superimposed. We used a projection along the z-axis, i.e., we assumed that z coincides with the LOS, and the collision takes place in the plane of the sky (x, y plane). The infalling cluster passes the main cluster center from below, moving east to west (second panel in the first row; run P40V18). The time is measured in Myr relative to the first core passage (t₀ = 0 Myr; left to right, first row: t = −475, −317, 0, and 238 Myr; second row: t = 396, 459, 555, and 713 Myr).

Zoom In Zoom Out Reset image size

The first two panels in the first row show the epochs before the first core passage, when only the outer regions of the gas in the two clusters take part in the collision. The shocks, generated in the same direction as the respective dark matter (the shock on the east moves to the west, the shock on the west moves to the east), move ahead of the two clusters. The third panel represents the epoch of the first core passage. The two dark matter peaks overlap at the center; they are slightly ahead of the X-ray peaks. The structure of the merging shocks change; multiple shocks are generated due to the collision of the denser gas in the merging clusters. The fourth panel shows the phase right after the first core passage, when there is only one X-ray peak, and the bow shock in the gas west of the main cluster is moving to the west ahead of the infalling cluster. The back shock in the east propagates to the east in the gas of the infalling cluster. Panels in the second row show phases after the first core passage, when there are two X-ray peaks. The first two panels show two shocks ahead of the cluster centers moving outward. The third panel shows the phase when the back shock has already left the gas of the infalling cluster. The bow shock is still moving outward (toward west). The fourth panel shows a later epoch, when both shocks have left the gas of the merging system, and the gas is falling back to the cores of the dark matter of the two components.

We study these outgoing merging shocks in detail because of their importance in determining the dynamical state of the merging cluster and the phase of the collision. The phase of the collision is important for particle acceleration models, because the physical properties of the shock depend on it, and these properties have a large impact on particle acceleration (e.g., Stroe et al. 2014; Fujita et al. 2016; Kang & Ryu 2016), as well as testing cosmological models (e.g., Molnar 2015).

The reason why the outgoing shocks constrain well the phase of the collision right after the first core passage in the outgoing phase is that they move fast and leave the gas as they propagate with high speed in the low-density gas at the outskirts of the system. Because of these relatively high outer shock speeds, these shocks can (re)accelerate particles and generate luminous radio relics for a limited time. As a consequence, shocks and bright relics can only be expected to be detected in a pair of merging clusters relatively soon after the first core passage, before the first turnaround, after which the relics become fainter without a shock to provide additional energy.

In general, as we see here, if the observed shocks were to lie in between the two X-ray peaks associated with relaxed cluster centers, we would conclude that the merging system is being viewed before the first core passage, as shown in the first panel of Figure 5, such as in the case of, e.g., Abell 1750 (Molnar et al. 2013a). A disturbed X-ray morphology, displaying one or two X-ray peaks within one or two shocks, instead suggests a system caught after first core passage (e.g., Mastropietro & Burkert 2008; Molnar & Broadhurst 2015).

Figure 5 also illustrates how fast the outgoing shocks propagate in a merging cluster similar to ZwCl008. After ∼500 Myr, both the bow shock and the back shock propagating in the gas of the main and the infalling cluster have left the system (panel 4 in Figure 5). The velocity of the bow shock relative to the ambient gas (the shock on the west moving west) is 4500 km s⁻¹; the back shock propagates faster (to the east) with 4900 km s⁻¹. At these speeds, the merging shocks (bow and back shocks) would have left the gas 470 and 380 Myr after the first core passage, but our numerical simulations suggest that they leave after 618 and 522 Myr. The turnover occurs 1.5 Gyr after the first core passage; thus, well before the turnover, both merger shocks leave the system. This is a general feature of merging shocks in most colliding clusters with moderate mass ratios and infall velocities $\gtrsim 1000$ km s⁻¹.

ZwCl008 has been suggested to be an older, less massive version of the Bullet Cluster (Golovich et al. 2017), where the less massive infalling cluster has pushed the gas of the main cluster out of equilibrium, and as a result, the X-ray morphology, which traces the gas of the main cluster, is irregular (X-ray feature in the east in the first panel of Figure 1). However, our hydrodynamical simulations suggest that the cometary-like X-ray peak in ZwCl008 near the center of mass of the smaller western cluster (see the first panel of Figure 1) marks the gas density maximum of the infalling cluster, similar to El Gordo (Molnar & Broadhurst 2015) and unlike the Bullet Cluster, in which the wedge-shaped bright X-ray feature marks a contact discontinuity. The shock in the Bullet Cluster is a faint X-ray feature ahead of the merging shock to the west (e.g., Markevitch & Vikhlinin 2007). Based on our best model, ZwCl008 is only 428 Myr past the first core passage and in a much earlier phase than the turnover, 1.5 Gyr, like the Bullet Cluster.

Our results suggest instead that, because of a larger impact parameter for the collision in ZwCl008 and the lower infall velocity relative to the Bullet Cluster, the ram pressure cannot fully hold back the gas relative to the dark matter in the infalling cluster, as opposed to the case of the Bullet Cluster. Instead, the infalling cluster in ZwCl008 passes relatively unhindered through the gas of the main cluster since it does not penetrate through the dense central gas of the main cluster.