Reconstruction of Radio Relics and X-Ray Tails in an Off-axis Cluster Merger: Hydrodynamical Simulations of A115

To understand the merger history of A115 in a way that is consistent with the observed features in radio and X-ray, we perform idealized hydrodynamics simulations of off-axis cluster mergers using RAMSES (Teyssier 2002). The Euler equations for hydrodynamics are solved using the HLLC scheme (Toro et al. 1994) with a Courant number of 0.5 and Min-Mod limiter in the slope computation. We design our collision simulation in an isolated box with a volume of 15³ Mpc³ using outflow (zero-gradient) boundary conditions. The cells are resolved based on their density from the lowest level (cell size of 0.46 Mpc) to the highest level of refinement (cell size of 7.32 kpc). With higher-resolution simulations (∼3.6 kpc), we verify that this resolution gives converging features that we present in this paper. An additional refinement criterion based on the pressure gradient (ΔP/P > 1 ; Teyssier 2002) is applied so that the shock structures propagating at the cluster periphery are also maximally resolved. This allows us to model merger shocks outside the virial radius of both clusters with ∼20 million leaf cells at the maximum refinement level, out of the total 30 million leaf cells. To reduce complexity, neither radiative cooling nor baryonic feedback is implemented. We believe that the exact morphology of the X-ray surface brightness distribution depends on the implementation details and thus should be a subject of future studies. Nevertheless, Kang et al. (2007) and Hong et al. (2015) show that these two processes have only minor impacts on shock properties that we focus on in the current paper.

Two clusters, the more massive main cluster (hereafter represented with subscript "m") and the less massive subcluster (subscript "s"), are generated to represent A115S and A115N at the cluster redshift, respectively. Both clusters are spherically symmetric, comprised of dark matter and the ICM. We let dark matter follow a Navarro–Frenk–White (NFW; Navarro et al. 1996) profile and the ICM a beta profile (Cavaliere & Fusco-Femiano 1976). The two profiles are:

$$\begin{eqnarray}&&{\rho }_{\mathrm{DM}}(r)=\displaystyle \frac{{\rho }_{0,\mathrm{DM}}}{r/{r}_{{\rm{s}}}{\left(1+r/{r}_{{\rm{s}}}\right)}^{2}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{ICM}}(r)=\displaystyle \frac{{\rho }_{0,\mathrm{ICM}}}{{\left(1+{r}^{2}/{r}_{{\rm{c}}}^{2}\right)}^{\tfrac{3}{2}\beta }},\end{eqnarray} \tag{ 2 }$$

where ${\rho }_{0,\mathrm{DM}}$ and ${\rho }_{0,\mathrm{ICM}}$ are the normalization constants, r_s = r₂₀₀/c is the scale radius of the dark matter halo, and r_c is the ICM core radius. We use β = 2/3 (i.e., isothermal profile) in our simulations. We adopt the scale radius r_s and total mass M₂₀₀ of dark matter halo from the latest WL study of Kim et al. (2019). We then determine the ICM density normalization ${\rho }_{0,\mathrm{ICM}}$ so that it satisfies the baryon mass fraction ${f}_{\mathrm{bar},200}\sim 0.13$ at r₂₀₀. We choose r_c = 40 kpc (30 kpc) for the main (sub) cluster. Both ICM core radii are small fractions (10%–16%) of the dark matter halo scale radii r_s and are motivated by the presence of the cool cores. Using the Disk Initial Conditions Environment (DICE; Perret 2016) code, we generate 10⁷ dark matter particles for each cluster and assign gas densities to the individual AMR cells. The temperature for each cell is computed in such a way that an hydrostatic equilibrium is satisfied at the given potential.

Since both Equations (1) and (2) give a diverging total mass, we make the density drop exponentially beyond the cutoff radii. Following the suggestion of Donnert et al. (2017), we choose r_cut = 2.2 Mpc and 1.8 Mpc for the main clusters and subclusters, respectively, for both DM and ICM profiles. The exponential decrease of the gas density stops when it reaches the background density ${\rho }_{\mathrm{back}}=1.5\times {10}^{-30}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. We assign T_back = 10⁶ K to the background temperature, which is also the minimum ICM temperature of each cluster. This is needed to achieve the pressure equilibrium at the boundary between the cluster and background.

We tested the dynamical stability of each halo by simulating it in isolation for 4 Gyr. We use the change of the two radii that enclose 50% of the DM and the ICM, respectively, as a measure of the stability. We find that both radii decrease by up to 10% during the first ∼2 Gyr, but afterward they remain constant. Because our collision is designed to take place at ∼4 Gyr after the initial separation, our simulated clusters are in the stable state at the impact.

Finally, the two clusters are positioned on the x–y plane with a separation of Δy = 4.3Mpc. The initial velocities are chosen to produce a collision velocity of ∼2000 km s⁻¹ at the impact. We also apply an initial offset in the x direction to obtain a desired pericenter distance, b ∼ 500 kpc. The collision velocity ∼2000 km s⁻¹ is slightly greater than the freefall velocity ∼1800 km s⁻¹ estimated with the timing argument (Sarazin 2002). In the Appendix, we show that this collision velocity is bracketed by the range of the values observed in the cosmological simulation. The initial condition parameters are summarized in Tables 1 and 2, and we refer to this setup as the reference run. Both the initial velocity and pericenter distance for this reference are chosen so that the results reasonably reproduce the observed high Mach number (${{ \mathcal M }}_{\mathrm{DSA}}\sim 4.58$) with sustained cool cores.

As mentioned in Section 1, the current paper focuses on the reproduction of the morphology, orientation, and location of both the X-ray emission and radio relic of A115. These results depend on parameters such as the collision velocity, pericenter distance, and viewing angle; for instance, a lower initial velocity results in a more prominent X-ray tail because the X-ray gas suffers from a higher ram pressure. Therefore, in order to investigate their effects, we run three additional simulations by modifying the initial velocity and impact parameter from the reference run. These include two simulations with modified velocities in the y direction by ±200 km s⁻¹ and one simulation with a larger (Δb = +600 kpc) impact parameter, which in turn increases the pericenter distance by 300 kpc. Hereafter, we refer to these three setups as the V+, V−, and b+ runs, respectively. We also investigate one case when the merger axis is not exactly perpendicular to the line-of-sight direction. We are aware that different combinations of various initial velocities, impact parameters, and viewing angles can lead to non-negligible differences in the final results. However, it is infeasible to exhaustively search for the optimal combination in this computationally expensive AMR study, which has an effective dynamic range in volume of 2048³.

3.2.1. X-Ray Mock Observation

The X-ray emissivity of each gas cell is calculated based on the following description of bremsstrahlung radiation from thermal electrons in a fully ionized gas with a primordial composition (Marinacci et al. 2018):

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{{\rm{X}} \mbox{-} \mathrm{ray}} & = & 4.90\times {10}^{-28}\,\mathrm{erg}\ {{\rm{s}}}^{-1}\\ & & \times {n}_{{\rm{H}}}^{2}{T}^{0.5}\left[\exp \left(-\displaystyle \frac{{T}_{1}}{{k}_{{\rm{B}}}T}\right)-\exp \left(-\displaystyle \frac{{T}_{2}}{{k}_{{\rm{B}}}T}\right)\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where n_H and T are the hydrogen number density (in units of cm⁻³) and temperature of each cell (in units of K), respectively. Note that we do not include metal lines, which can non-negligibly alter the results in the low-temperature region, but are not likely to significantly modify the results in the high-temperature region that we focus on in the current study. We choose T₁ = 0.5 keV and T₂ = 7 keV in order to match the energy range used in our Chandra image creation. The X-ray surface brightness map is obtained by integrating the emissivity contributed from every cell along the line-of-sight (z-axis) direction.

The Chandra temperature map shown in the middle panel of Figure 1 (Hallman et al. 2018) is estimated from the local X-ray spectral fitting. To simulate this Chandra spectroscopic temperature map, we average the temperature of each cell using the weight suggested by Mazzotta et al. (2004):

$$\begin{eqnarray}&&{T}_{\mathrm{spec}}=\displaystyle \frac{{\sum }_{\mathrm{LOS}}{n}_{i}^{2}{T}_{i}^{0.25}{{dx}}_{i}^{3}}{{\sum }_{\mathrm{LOS}}{n}_{i}^{2}{T}_{i}^{-0.75}{{dx}}_{i}^{3}},\end{eqnarray} \tag{ 4 }$$

where n_i, T_i, and dx_i represent the gas density, temperature, and cell size of the ith cell along the line-of-sight (z-axis) direction.

3.2.2. Shock Detection Algorithm

We adopt the shock criteria suggested by Ryu et al. (2003) and follow the implementation for the AMR structure described by Skillman et al. (2008). We provide a brief summary below.

A cell is flagged as a shock center candidate if its gas velocity converges (∇ · v < 0). We walk through the candidate cells in the ±∇T directions and find the cell that first violates either of the two conditions: ΔT × ΔS > 0 or $\mathrm{log}\left|{\rm{\Delta }}T\right|\gt 0.11$, where ΔT and ΔS are the temperature and entropy differences between the two adjacent cells. The cell that first violates either of the two conditions is classified as a post-shock (pre-shock) cell if the search direction is ∇T > 0 (<0). Prior to identifying both pre-shock and post-shock cells, if the divergence, ${\rm{\nabla }}\cdot {\boldsymbol{v}}$, of the visited cell is found to be lower than the shock center candidate that initiated this walk, the shock center candidate is discarded. Once both post- and pre-shock cells are identified, the shock center candidate is promoted to a shock center with a Mach number (${ \mathcal M }$) derived from the following relation:

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

where T₁ and T₂ are the temperatures of the pre- and the post-shock cells, respectively. Thus, each shock structure is a column of cells comprised of one shock center, one post- and one pre-shock cell, and multiple shock cells; see Figure 1 of Skillman et al. (2008) for the schematic description. Note that we must perform this shock identification with the smallest AMR cells in order to avoid facing multiple adjacent cells. However, we verify that no cells at a coarser level satisfy our shock conditions described above. This is attributed to our refinement criteria ΔP/P > 1 , which ensures that the shock cells are always maximally resolved.

3.2.3. Generation of Mock Radio Observation

The radio emissivity is estimated using the prescription of Hong et al. (2015). For each shock structure, we first compute the kinetic energy flux f_k, the main power source of the radio emission, as follows:

$$\begin{eqnarray}&&{f}_{{\rm{k}}}=\displaystyle \frac{1}{2}\rho {\left({ \mathcal M }{c}_{{\rm{s}}}\right)}^{3},\end{eqnarray} \tag{ 6 }$$

where c_s is the sound speed of the pre-shock cell given by c_s = (γP/ρ)^1/2 with γ = 5/3. We then calculate the efficiency of acceleration as an increasing function of the Mach number ${ \mathcal M }$. Kang & Ryu (2013) show that the dissipation efficiency $\eta ({ \mathcal M })$ becomes negligible for ${ \mathcal M }\lt 3$, approaches ∼0.8% at ${ \mathcal M }=3$, and saturates to ∼20% at ${ \mathcal M }\gt 10$. This determines the momentum distribution of cosmic-ray electrons (CRe) for each shock structure:

$$\begin{eqnarray}&&{f}_{\mathrm{CRe}}(p)\propto {K}_{{\rm{e}}/{\rm{p}}}\eta ({ \mathcal M }){f}_{{\rm{k}}}{p}^{-q},\end{eqnarray} \tag{ 7 }$$

where q = 3σ/(σ − 1) is the spectral index, σ is the compression rate given by $\sigma =(\gamma +1){{ \mathcal M }}^{2}/[(\gamma -1){{ \mathcal M }}^{2}+2]$, p is the momentum of electron, and K_e/p is the CR electron-to-proton ratio, which is set to 0.01 (Hong et al. 2015). Following Hong et al. (2015), we also apply the exponential cutoff due to cooling in the momentum distribution. Finally, the radio emissivity (j_ν) is obtained from

$$\begin{eqnarray}&&{j}_{\nu }\propto {\int }_{{p}_{\min }}^{\infty }{BF}\left(\displaystyle \frac{\nu }{{\nu }_{{\rm{c}}}(p)}\right){f}_{\mathrm{CRe}}{d}^{3}p,\end{eqnarray} \tag{ 8 }$$

where B is the magnetic field strength of the post-shock cell, ν_c is the characteristic frequency, and $F(x)\equiv x{\int }_{x}^{\infty }d\eta {K}_{5/3}(\eta )$ with K_5/3(η) being a modified Bessel function. In our study, we assume a constant magnetic field of B = 1 μG, which is derived by Govoni et al. (2001) from the equipartition assumption at the A115 radio relic. The minimum momentum p_min contributing to the radio emission is set to ${p}_{\min }\,=\,0.01{m}_{{\rm{p}}}c$, which assumes that DSA acts on $p\gtrsim 3{p}_{\mathrm{th},{\rm{p}}}$, where ${p}_{\mathrm{th},{\rm{p}}}$ is the most probable momentum of thermal protons. Accordingly, radio emissivity is calculated along the downstream regions of the identified shock structures (i.e., from the shock center to the post-shock cell). The projected radio flux map obtained in this way is smoothed to match the beam size of the VLA-C observation (15'' × 15''; Botteon et al. 2016). We refer readers to Hong et al. (2015) for more details on the individual steps in estimating the radio emission.

4.1.1. X-Ray observables

Figure 2 shows the time evolution of our simulated off-axis cluster merger presenting both projected X-ray surface brightness and temperature distributions. Overall, the result is consistent with general expectations from previous studies (e.g., Ricker & Sarazin 2001). The X-ray luminosity of the system rises as they approach, peaks at the closest passage (t = 0 Gyr), and decreases afterward. The temperature in the intermediate region between the two clusters increases as the two clusters approach each other. The maximum temperature reaches up to > 10 keV when the two clusters are near their impact (t = −0.1 Gyr). Because of the large pericenter distance (∼500 kpc), the two clusters maintain their cool cores after the impact. We note that temperature folds exist near the density cutoff regions, although we employ a large cutoff radius (yellow arrow in Figure 2). Thus, interpretation of the temperature features beyond the cutoff radius should use caution. Apart from these temperature folds, both the X-ray luminosity (${L}_{0.1-2.4\mathrm{keV}}\sim 1.3\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$) and the projected temperature of the intermediate region after the closest passage (t = [0.1,0.4] Gyr) agree with the observed level ${L}_{0.1-2.4\mathrm{keV}}=1.5\times {10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (Botteon et al. 2016) and T = 11.03 ± 1.7 keV (region D in Table 2 of Hallman et al. (2018)). This agreement in luminosity must be interpreted with caution, however, because it depends on how exactly one models AGN feedback and radiative cooling, neither of which are included in the current study.

Zoom In Zoom Out Reset image size

We note that the "S"-shaped feature of the central high-temperature region seen at or near the closest encounter (t = ±0.1 Gyr) is not clear in the observed temperature map (Figure 1). Assuming that the S/N value of the observed temperature map is sufficient, we suspect that the difference may be attributed to the intrinsic inhomogeneity of A115, viewing angle, and/or the lack of numerical instability, all of which can blur the boundary of the "S"-shape feature.

For both clusters, cometary tails can be identified from tapering surface brightness (${{\rm{\Sigma }}}_{\mathrm{Xray}}\gtrsim {10}^{-4}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$) as well as low-temperature regions (T ≲ 2 keV; Figure 2). These features start to develop when the two clusters are near the core passage (t = −0.1 Gyr) and become more distinct at t = 0.1 Gyr. The direction of the tails traced by the surface brightness are antiparallel to the initial velocities. However, at t ∼ 0.4 Gyr after the core passage, the tails are somewhat shortened, with their directions nearly perpendicular to the initial velocities. This X-ray tail rotation is known to take place as the gas that was stripped by ram pressure earlier begins to fall back to the cluster potential with a nonzero angular momentum. This results in a "slingshot" feature in the X-ray images (Sheardown et al. 2019) that shows a morphology different from the ram pressure stripping tail observed at t = 0.1 Gyr. Remember that, at this epoch (t ∼ 0.4 Gyr), the two cluster are still moving away from each other (they reach their apocenters much later, at t ∼ 2.0 Gyr) and the orientation of the t ∼ 0.4 Gyr velocity vectors are within  20° from that of the collision axis. The time evolution of the subcluster X-ray tail orientation in comparison with the velocity vector is further discussed with shock structures in Section 4.1.2.

Dissociation between dark matter halo and gas is often reported in cluster mergers (e.g., Clowe et al. 2006; Dawson et al. 2012). The gas-mass offset is consistent with our expectation in our cold dark matter (CDM) paradigm, because the plasma is subject to Coulomb forces whereas the dark matter is believed to be collisionless. In off-axis collisions, the offset is expected to be smaller and decrease with impact parameter. Both mass peaks of A115 from our WL analysis are well-aligned with the X-ray peaks, as shown in Figure 1. From the current hydrodynamical simulation, we confirm that no significant offsets between the dark matter and X-ray peaks are present during the period t = [−0.1,0.4] Gyr.

4.1.2. Merger Shocks and X-Ray Tail Alignments

The shock structures found in the simulation results are displayed in Figure 3 using the integrated kinetic energy flux (F_k, top panel) and the average Mach number in the central Δz = 600 kpc slice (bottom panel). The integrated kinetic energy flux, F_k, is obtained by integrating f_k × A_sh along the line-of-sight direction, where A_sh is the shock surface area.

Zoom In Zoom Out Reset image size

As shown by Ha et al. (2018), the shock layers (see bottom panel of Figure 3) can be categorized into three groups: (1) the merger shocks propagating along the merger axis in front of the X-ray cores, (2) the equatorial shock with high Mach numbers (${ \mathcal M }\gt 10$) expanding in a direction perpendicular to the merger axis, and (3) the spherical accretion shock at the outskirt of the cluster.⁷

Among the three shock groups, the integrated kinetic energy flux, F_k is dominated by the merger shocks and reaches the highest values along the front edge (see top panel of Figure 3). The kinetic energy flux of the northern merger shock extends more than 1 Mpc and diminishes precipitously at the eastern density cutoff. As the kinetic energy flux f_k powers the radio emission, we can expect that the radio emission will also be dominated by the merger shocks and become brightest along the edge of the subcluster merger shock. Scrutiny of the eastern side of the northern merger shock reveals a sudden increase in the Mach number toward the eastern edge. This feature is attributed to our assignment of a low temperature to the cells in the density cutoff region. Considering the small temperature gradient of the β-profile, we think that this artificial temperature drop at the density cutoff may lead to overestimation of the Mach number (from ∼4 to ∼7), which can in turn non-negligibly boost the acceleration efficiency (Kang & Ryu 2013). Also, this artifact may result in the premature truncation of the predicted radio relic and overestimation of the radio flux on the eastern side (see Section 4.2).

Understanding how the orientations of the X-ray tails and merger shocks evolve at different stages of an off-axis merger with respect to the collision axis and the halo velocity vector is one of the key goals of the current study. From visual inspection, one can notice that, unlike the orientation of the X-ray tail of the subcluster, the orientation of the subcluster merger shock does not vary greatly during the t = [−0.1,0.4] Gyr period. In order to discuss the orientation evolution quantitatively, we use the second moments of the features to evaluate the orientation angle. Employing a proper window function is needed in this calculation because the orientation depends on the scale (in particular for the X-ray tail). For the integrated kinetic energy flux, we consider the data where the values are within the top 0.1%, which defines an ellipse whose major axis is ∼1 Mpc. The threshold is chosen to be $\gt 2\times {10}^{-4}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ for the X-ray surface brightness. The resulting major axis of the ellipse is ∼300 kpc. We determine the cluster motion from the average velocity vector of the dark halo. In Figure 4, we display the orientation angle as a function of time since the impact.

Zoom In Zoom Out Reset image size

As qualitatively described before, the direction of the X-ray tail (red arrow in the top panel and red solid line in the bottom panel of Figure 4) is nearly aligned (∼30°) with the collision axis at the time of the impact (epoch B). However, its deviation from the collision axis quickly increases in a short period of time and reaches > 90° at epoch C (t = 0.3 Gyr). Because of the large collision velocity, the subcluster velocity vector (see the black solid line in the bottom panel of Figure 4) remains well-aligned with the collision axis until t = 0.6 Gyr within ∼20°.

As expected, the orientation of the merger shock is approximately perpendicular to the collision axis from the onset of the merger (epoch B) and changes only by a small amount throughout the merger duration t = [−0.2,0.6] Gyr investigated here (see the blue solid line in the bottom panel of Figure 4).⁸ Between epochs B and C, the shock orientation also undergoes a rotation, which however is much smaller than that of the X-ray tail.

Again, this quantitative analysis confirms that the orientation of the X-ray tail is not a reliable indicator of the subcluster motion in an off-axis merger, as already demonstrated by Sheardown et al. (2019). On the other hand, merger shocks, being relatively insensitive to the stage of the merger, provide more reliable information on the direction of the collision axis.

The difference in the rotation rate between the X-ray tail and the merger shock gives an epoch where the orientations of the two structures align (t ∼ 0.1 Gyr; see the middle panel of Figure 4). Because we witness a similar tail-shock alignment in the current A115 observation, this issue deserves a further test. One can imagine a number of parameters affecting the exact alignment epoch. However, here we mainly focus on its dependence on the choice in our window function. Lowering the threshold of the X-ray tail selection gives a larger area for the evaluation of the orientation. As one can expect from the X-ray map, it decreases the angle between the tail orientation and the collision axis. Specifically, when we lower the threshold by a factor of two from $2\times {10}^{-4}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ to $1\times {10}^{-4}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$, the decrease is ∼30° when averaged over different epochs. However, the rotation "rate" remains similar, and thus the total amount of the orientation change Δθ is also similarly large (∼146° ± 87°) even with this new window function choice.⁹ We adopt this window function–dependent change of ∼30° as a systematic uncertainty in our determination of the orientation angle. For comparison with the observation, we applied the same procedures to our X-ray and radio data, except that we masked out the point sources. Our analysis shows that the observed alignment occurs 0.1–0.4 Gyr after the collision.

We find that the projected radio flux based on the procedure summarized in Section 3.2.3 is lower than the noise level $1\sigma =70\,\mu \mathrm{Jy}\,{\mathrm{beam}}^{-1}$ of the VLA-C observation (Botteon et al. 2016). Therefore, in our presentation (Figure 5) we arbitrarily scale up the radio flux by a factor of 10 (see Section 4.2.1 for discussion). Readers are reminded that the radio flux contours in Figure 5 and hereafter are displayed with solid (dashed) lines for the 3 and 6σ (1 and 2σ) levels.

Zoom In Zoom Out Reset image size

Figure 5 shows the northern and southern radio relics along the edge of the simulated merger shocks seen in Figure 3. The orientation of the northern radio relic roughly agrees with that of the subcluster merger shock, however, with some notable differences. At the 0.1 < t < 1.0 Gyr epoch, during which the peak of the radio flux is > 3σ, the orientation of the relic remains nearly constant (Δθ ≲ 3°).¹⁰ The merger shock undergoes a larger change (Δθ ≲ 20°) in orientation. This is because the part of the subcluster merger shock ahead of the subcluster has a low Mach number (∼2, see Figure 3), which results in a low radio emissivity. The nearly constant radio relic orientation angle, ∼62° with respect to the collision axis, is similar to the merger shock orientation at t ∼ 0 Gyr. Thus, our simulation result illustrates that the orientation of the observed radio relic might serve as a reliable indicator of the merger axis, insensitive to the stage of the merger.

Our simulated relic is somewhat shorter than the observed radio relic. One might suspect that this may arise from the eastern side truncation mentioned in Section 3.2.2. However, one should also remember that the Mach number in this region is likely to be overestimated because of the same truncation effect (thus the relic on this side is estimated to be brighter). Correcting for this Mach number overestimation would result in a factor of ∼4 reduction in radio flux, providing a smoothly decreasing radio flux toward the eastern boundary before the truncation.

4.2.1. Enhancement of Radio Relics through Stronger Magnetic Fields and Fossil Electrons

As mentioned before, the reconciliation in radio flux between our simulation and the observation requires a factor of 10 boost in the simulated radio flux. This supports the well-known claim that the acceleration of electrons from the thermal pool (via DSA) is too inefficient to explain the observed luminosity of typical radio relics (e.g., Botteon et al. 2020). One suggestion to resolve the issue is the so-called "reacceleration model," which hypothesizes that fossil relativistic electrons from radio galaxies are reaccelerated by the DSA mechanism. The theory has been gaining support with the discovery of a growing number of clusters hinting at such radio galaxy–relic connections (e.g., van Weeren et al. 2017).

Perhaps A115 can also serve as such an example, because there are a few confirmed radio sources blended with the radio relic (Harwood et al. 2015). In this paper, we test the enhancement of the radio relic flux by assuming a factor of 3 stronger magnetic field (B = 3 μG) with the availability of significant fossil electrons. We implement the impact of fossil electrons through a higher dissipation efficiency ${\eta }_{\mathrm{eff}}({ \mathcal M })$. As shown in Equation (6), the stronger magnetic field increases the overall radio emissivity, whereas the modified dissipation efficiency assuming the pre-existing CR (i.e., P_CR = 0.05P_th) enables the weak shocks (${ \mathcal M }\sim 2$) to accelerate particles (Kang & Ryu 2013).

We display the resulting radio flux map at t = 0.1 Gyr in Figure 6. With the aforementioned recipe, this new radio flux matches the observed level without the arbitrary 10× upscaling employed in Figure 5. Moreover, the northern relic is stretched further to the west than the previous case, because the enhanced dissipation enables weak shocks in the western region to emit radio. This increase in length by a factor of two also better matches the observation. In the eastern shock region, the simulated radio signal remains high (>6σ) even after we correct for the overestimated Mach number, which implies that the northern radio relic could potentially extend further beyond the truncation and reach the observed length of ∼2 Mpc, if there were no boundary effect.

Zoom In Zoom Out Reset image size

Figure 6 shows that the southern main cluster merger shock also becomes detectable. The absence of the corresponding signal in the current radio observation may be attributed to the lack of radio point sources in the southern region, if indeed fossil electrons are mainly responsible for the enhancement of the radio relic flux.

We adopt a uniform magnetic field strength (1 μG) to estimate our radio emissivity. Considering the rather turbulent nature of the ICM magnetic fields (e.g., Roh et al. 2019), we think that it is probably unrealistic to assume such homogeneity over the ∼2 Mpc scale. Moreover, if this radio relic was formed by reacceleration of low-energy cosmic-ray electrons in a dead radio jet, which has been mixed with entrained ICM, the fossil plasma may possess remnant magnetic fields as well. Detailed modeling of these features is beyond the scope of the current paper.

In addition to the reference run, we test three additional cases, where the collision velocity is lower (V− run), the collision velocity is higher (V+ run), and the pericenter distance is larger (b+ run). Figure 7 presents the resulting X-ray density and temperature maps overlaid with radio contours enhanced through the same recipe explained in Section 4.2.1. Because these three simulations are now at different stages for the same simulation time t, we compare the epochs when the distance between the clusters is ∼900 kpc after the core passage.

Zoom In Zoom Out Reset image size

On large scales, the three simulations in Figure 7 share similar distributions in density and temperature to those of the reference run, featuring cool cores and hot intermediate regions. On small scales, a few notable differences are present. First, the X-ray tails become more prominent when the collision velocity is lower. We attribute this to the smaller cluster separation at the impact, which results in the stronger drag force. Second, the length of the radio relic also increases for the lower velocity because more space is accessible to the merger shock.

The b+ run is selected at t = 0 Gyr since the smallest separation for this case is already comparable to the observed value ∼900 kpc. The simulated X-ray map shows a distinct ram pressure stripping tail, which approximately corresponds to the merger scenario proposed by Hallman et al. (2018). The triangular shape of the cometary tail and the hot temperature in the intermediate region resembles the observed features. Nevertheless, the anticipated merger shock, even after we consider reacceleration of cosmic rays, is too weak to produce the observed radio flux. Furthermore, the expected peak location of the radio emission is different from the observed position.

The radio relics in Figure 7, except for the b+ run, all show an apparent alignment with the X-ray tails. For the three simulation runs with a distinct radio relic, we compare the following four physical properties with the observations: the distance between the subcluster merger shock and X-ray core D_shock, the temperature in between the clusters, the distance between the two X-ray cores D_core, and the alignment between the subcluster X-ray tail and the subcluster merger shock Δθ in Figure 8.

Zoom In Zoom Out Reset image size

The distance to the merger shock D_shock is computed using the separation between the kinetic energy flux (f_k)-weighted center of the merger shock and the X-ray core for the simulations. For the observation, we use the radio flux-weighted center of the radio relic and the X-ray center of A115N. The radio point sources are masked out when computing the flux-weighted center of the observed radio relic.

The angle Δθ is measured in the same way as in the creation of Figure 4 for both simulations and observations. The systematic uncertainty due to the window function choice is considered to be +30°, as before.

The central temperature in the simulation is computed using Equation (4), with the cells within the 400 kpc × 200 kpc box positioned between the clusters and with the shorter side aligned parallel to the axis connecting the two clusters. Within the box, we averaged the projected temperature to derive the central temperature of the simulation. To determine the temperature from the observation, we first extract the X-ray spectrum from the equivalent region (∼2' × 1') using ciao-4.11; we merged the X-ray spectra from the observations 3233, 13458, 13459, 15578, and 15581, and binned the counts in such a way that each energy bin has a signal-to-noise ratio of 5. We fit an absorbed MEKAL model (Mewe et al. 1985; Kaastra & Mewe 1993; Liedahl et al. 1995) to the binned spectrum using the 0.5–7 keV energy range and determine the temperature to 7.46 ± 0.57 keV. We verify that a highly consistent result is obtained when we use the Astrophysical Plasma Emission Code (APEC; Smith et al. 2001) plasma model instead of the MEKAL model.

The shock distance (measured from the subcluster X-ray core) D_shock slowly decreases until it reaches the minimum at t ∼ 0.1–0.3 Gyr (Figure 8). Although the exact values differ among the simulation setups, they all are similar to the observed value (dashed). After this epoch, D_shock gradually increases with a moderate slope, which indicates the launch of the merger shock (Sheardown et al. 2018).

As one could expect, the temperatures in between the clusters show a gradual decrease after they peak at the impact. The peak temperatures are similar to the observed value. Nevertheless, one should not use this temperature comparison to constrain the merger scenario/epoch, because a few keV difference in temperature can be easily resolved by "fine-tuning" simulation setups, which are beyond the scope of the current paper.

The distances between the cores D_core match the observed value at t ∼ 0.3 Gyr for the reference and V+ runs, and at t ∼ 0.4 Gyr for the V− run.

Considering the systematic uncertainty due to the window function, we can estimate that the alignments between the X-ray tail and merger shock happen during the t = [0,0.4] Gyr period.

In summary, our simulations show that the three observed properties of A115 can be reproduced at t ∼  0.3 Gyr after the core impact.

Our simulations described in Sections 4.1 and 5.1 show that many observed features such as the morphology, location, and orientation of the X-ray tail and radio relic can be reproduced by an off-axis merger ∼0.3 Gyr after the core passage, which favors a "slingshot" mechanism rather than a ram pressure stripping for the origin of the observed X-ray morphology. These simulations assume a plane-of-the-sky merger, which is motivated by the small line-of-sight velocity difference between the member galaxies of A115N and A115S (Golovich et al. 2019). However, detailed analysis of the velocity structure near the cluster cores with the cluster galaxies and gas hints at the possibility that the A115 merger may have a significant LOS velocity component (Barrena et al. 2007; Kim et al. 2019). For example, Kim et al. (2019) show that the relative velocity between the two BCGs of A115N and A115S is large (853 ± 6 km s⁻¹). They also demonstrate the robustness of this result by including in their analysis the member galaxies within 250 kpc radius from each BCG, although the significance is somewhat lower (838 ± 549 km s⁻¹). Using X-ray spectra, Liu et al. (2016) report an extremely large velocity difference (∼4600 ± 1100 km s⁻¹) between A115N and A115S. Although we think that more studies are needed to settle the A115 LOS velocity issue, here we discuss how the results change when we assume a nonzero viewing angle.

If we hypothesize that the relative LOS velocity between the BCGs  (∼850 km s⁻¹) is representative of the LOS velocity difference between A115N and A115S, the observed projected distance ∼900 kpc between the X-ray cores and the relative LOS velocity can be obtained at t = 0.4 Gyr by rotating the merger axis of the reference run by ∼30°. Figure 9 shows the resulting X-ray density map and the radio relics.

Zoom In Zoom Out Reset image size

When rotated by ∼30°, the X-ray morphology shown in Figure 9 features an enhanced bridge between the two clusters. Also, the X-ray tails look somewhat more compact than the unrotated version, more similar to the observation (Figure 1).

The peak radio flux is reduced by ∼80%, although the overall flux and morphology still resemble the radio data. The properties described in Figure 8 do not change very much because of the small viewing angle; the angle Δθ is ∼50° and the distance to the shock center D_shock from the X-ray core is ∼500 kpc. In this rotated projection, the orientation of the radio relics is still nearly perpendicular to the cluster motion, whereas the projected X-ray tail becomes further rotated because the epoch selected here (to match the projected distance between the cores) corresponds to a later merger phase. The rotated temperature map (not included in the current paper) shows a somewhat reduced temperature contrast compared to Figure 2. However, because of the small viewing angle considered here, the overall difference is insignificant.

From this simple experiment, we can reason that the introduction of a small viewing angle (∼30°) still provides a consistent picture of the merger scenario of A115. However, we believe that a viewing angle significantly larger than this (e.g., Liu et al. 2016) cannot reconcile with the observation. A larger viewing angle requires a larger three-dimensional separation between A115N and A115S. Therefore, we do not expect that any noticeable X-ray tails with the current observed orientations would be produced at this large separation (i.e., at the very late phase). Of course, the formation of radio relics matching the current radio data would be even more challenging.

Our simulations show that the observed morphology, orientation, and location of both the X-ray emission and radio relic of A115 can occur during the t = 0.1–0.3 Gyr period after an off-axis collision, before the two subclusters reach their apocenters. In particular, this scenario can explain the parallel alignment of the radio relic and the northern X-ray tail within the conventional paradigm, in which the relic traces the location of the merger shock. Although the details change when different initial velocities and impact parameters are selected, this alignment is always reproduced in our experiments with a zero or modest viewing angle. We present the schematic diagram showing our merging scenario in Figure 10.

Zoom In Zoom Out Reset image size

Golovich et al. (2019) suggested that the two subclusters in A115 might be in their returning phase. This scenario is motivated by the direction of the observed X-ray tails, which are interpreted as ram pressure features. The main drawback of this scenario is its inability to explain the proximity of the relic to the northern X-ray core. As our simulations show, the onset of the second infall happens at t = 1.8 Gyr, by which time the merger shock is located well outside the simulation box.

The Hallman et al. (2018) scenario is similar to that of Golovich et al. (2019), in that they attribute the X-ray tails to ram pressure stripping. The difference is that they explain that "the radio relic results from the Fermi acceleration of a relic population of cosmic rays ejected from the local radio galaxies and advected behind the motion of A115N." Although the authors are implicit regarding the merger phase, the proximity of the relic requires the relic to be young. Hence, although we agree that the local galaxies might be associated with the radio relic, it is difficult to produce such a fully developed merger shock at this early phase.

We cannot completely rule out the possibility of the ram pressure stripping tail scenario because this scenario can better reproduce the triangular morphology of A115N X-ray tail (see the b+ run in Figure 7). However, the difficulty in this scenario is that our estimation of the radio emission via merger shock is incompatible with the observed radio emission, in both flux level and the location. This challenge (if one favors a ram pressure scenario) can be overcome if what we have called the radio relic so far might, in fact, have a different origin. For example, as mentioned in Botteon et al. (2016), energized trails of the old AGN tails (i.e., radio phoenix; de Gasperin et al. 2017) may mimic a relic-like morphology as observed. Future deeper radio observations can test this possibility.

Despite the ease with which we can explain the key observed features, our favored merger scenario has some limitations regarding detailed structures. Our cluster merger is modeled as a collision between the two spherically symmetric halos, which inevitably produces a smooth merger shock. On the other hand, the merger shocks in both cosmological simulations and observations show inhomogeneous Mach numbers (e.g., Hong et al. 2015) and filamentary substructures (e.g., Vazza et al. 2016; Rajpurohit et al. 2018). Although the main X-ray and shock properties are determined by the cluster collision itself, by and large (Kang et al. 2007), real-world inhomogeneities in halo and background properties can modify the simulated features. For instance, turbulence in the X-ray tail can be generated (e.g., Donnert et al. 2017). In our simulation, the turbulent instability along the shear of the X-ray tail (i.e., Kelvin–Helmholtz instability) is only weakly developed. As other studies with comparable or larger resolution also show similar features to our results (e.g., Brzycki & ZuHone 2019), we speculate that perhaps the limited resolution might have prevented a full development of the instability. Nevertheless, since our analysis focuses on the X-ray tail near the center, where the instability plays a negligible role (see Figure 2 in Sheardown et al. 2019), we expect our conclusions to remain valid even with the presence of instabilities.

The X-ray tail in our favored merger scenario is somewhat curved, whereas this feature is not clear in the Chandra observation. As mentioned earlier, this morphological discrepancy can be mitigated with a nonzero viewing angle (Figure 9). Obviously, an introduction of contrived inhomogeneous environment might also produce X-ray features more similar to the observation.

The model used in our estimation of the radio flux is rather limited, since we adopt a simple DSA prescription for the proton acceleration at parallel shocks (Kang & Ryu 2013). The shock criticality and obliquity are not considered, because we do not follow the evolution of magnetic fields in our hydrodynamic simulations. Recently, Kang et al. (2019) suggested through particle-in-cell simulations that electrons may be injected to the DSA process only at supercritical (${{ \mathcal M }}_{{\rm{s}}}\gt 2.3$), quasi-perpendicular shocks in high plasma beta (=P_g/P_B) ICM plasma. Future studies using MHD simulations will be useful to quantifying the electron acceleration at ICM shocks driven by off-axis cluster mergers.