Runaway Merger Shocks in Galaxy Cluster Outskirts and Radio Relics

TL;DR

This paper investigates how runaway shocks in galaxy cluster outskirts can sustain or amplify their strength due to steep density gradients, potentially powering radio relics with confined radio emission regions.

Contribution

It introduces the concept of runaway shocks maintaining or increasing strength in steep density profiles, explaining the formation of narrow radio relics in cluster outskirts.

Findings

Runaway shocks can persist or strengthen in steep density gradients beyond R500.

Such shocks produce narrow regions of enhanced radio emission due to adiabatic compression.

Steep gas density profiles facilitate the survival of moderately strong shocks at large radii.

Abstract

Moderately strong shocks arise naturally when two subclusters merge. For instance, when a smaller subcluster falls into the gravitational potential of a more massive cluster, a bow shock is formed and moves together with the subcluster. After pericenter passage, however, the subcluster is decelerated by the gravity of the main cluster, while the shock continues moving away from the cluster center. These shocks are considered as promising candidates for powering radio relics found in many clusters. The aim of this paper is to explore the fate of such shocks when they travel to the cluster outskirts, far from the place where the shocks were initiated. In a uniform medium, such a "runaway" shock should weaken with distance. However, as shocks move to large radii in galaxy clusters, the shock is moving down a steep density gradient that helps the shock to maintain its strength over a large distance. Observations and numerical simulations show that, beyond $R_{500}$, gas density profiles are as steep as, or steeper than, $\sim r^{-3}$, suggesting that there exists a "Habitable zone" for moderately strong shocks in cluster outskirts where the shock strength can be maintained or even amplified. A characteristic feature of runaway shocks is that the strong compression, relative to the initial state, is confined to a narrow region just behind the shock. Therefore, if such a shock runs over a region with a pre-existing population of relativistic particles, then the boost in radio emissivity, due to pure adiabatic compression, will also be confined to a narrow radial shell.

Tables

Table 1. Table 1 : Parameters of 1D blast-wave simulations.

IDs	$\omega$666The parameter used in Eq. (5) to set the slope of the initial gas density profile.	$\xi$777The factor $\xi$ by which the gas pressure in the innermost cell is increased to initiate the blast wave.	${ℳ}_s^{100}$	${ℳ}_s^{950}$888${ℳ}_s^{100}$ and ${ℳ}_s^{950}$ are the Mach numbers of the blast wave when it reaches $R_s=100$ and $950\mathrm{kpc}$, respectively.
W1P5	$1$	${10}^5$	1.28	1.03
W1P6	$1$	${10}^6$	2.17	1.08
W1P7	$1$	${10}^7$	5.88	1.26
W2P5	$2$	${10}^5$	1.51	1.09
W2P6	$2$	${10}^6$	3.15	1.30
W2P7	$2$	${10}^7$	9.18	2.32
W3P5	$3$	${10}^5$	1.89	1.30
W3P6	$3$	${10}^6$	4.77	2.43
W3P7	$3$	${10}^7$	14.6	7.11
W4P5	$4$	${10}^5$	2.59	2.21
W4P6	$4$	${10}^6$	7.45	6.81
W4P7	$4$	${10}^7$	22.9	21.7

Equations

$${\eta }_{\mathrm{s}}\equiv \frac{\mathrm{d}\ln u_{\mathrm{g}}}{\mathrm{d}\ln R_{\mathrm{s}}}=\frac{\omega }{2}-1,$$

$$u_{\mathrm{g}}\propto \frac{1}{R_{\mathrm{s}}\sqrt{\ln (R_{\mathrm{s}}/\alpha )}},$$

$${\eta }_{\mathrm{s}}=\frac{\omega }{4}-1.$$

$${\eta }_{\mathrm{s}}=\frac{\omega }{2}-\frac{3}{2},$$

$$\rho_{\rm gas,0}(r)=\rho_{\rm c}\Big{(}\frac{r}{r_{\rm core}}+1\Big{)}^{-\omega},$$

$${\int }_0^r{\rho }_{\mathrm{gas}}(r^{\prime })r^{\prime 2}\mathrm{d}{r}^{\prime }={\int }_0^{r^{\ast }}{\rho }_{\mathrm{gas},0}(r^{\prime })r^{\prime 2}\mathrm{d}{r}^{\prime }.$$

Full text

††pubyear: 2019††pagerange: Runaway Merger Shocks in Galaxy Cluster Outskirts and Radio Relics–References

1 Introduction

Shocks naturally form during the merger process of galaxy groups and clusters and frequently lead to remarkable observational signatures, e.g. sharp X-ray surface brightness edges, gas temperature discontinuities, and radio relics (see Markevitch & Vikhlinin 2007; Feretti et al. 2012; Bykov et al. 2015; van Weeren et al. 2019 for reviews). In X-ray observations, these shocks are usually found well within the cluster virial radius, and have a moderate Mach number ($\mathcal{M}_{\rm s}\sim{\rm few}$), in contrast to the cluster virial shocks ($\mathcal{M}_{\rm s}\sim 10-10^{3}$; see e.g. Borgani & Kravtsov 2011).

These merger shocks have also been extensively explored through numerical methods (e.g. Ricker & Sarazin 2001; Vazza et al. 2009; Paul et al. 2011; Ha et al. 2018). In general, the evolution of merger shocks in galaxy clusters could be summarized as a two-phase process (see Fig. 1 for an example):

• •

In the initial “driven phase”, a bow shock is formed ahead of the gaseous core of the infalling subcluster111The boundary separating the gaseous atmosphere of the main cluster from that of the subcluster is a contact discontinuity, a.k.a. cold front.. Zhang et al. (2019) have made a systematic study of this phase, and have shown that the velocities of the cold front, shock, and dark matter (DM) halo of the subcluster could be quite different, especially after the passage of the subcluster through the main cluster core. This is caused by (1) the acceleration/deceleration of the subcluster while moving in the gravitational potential well of the main cluster; and (2) the density/pressure gradients of the intracluster medium (ICM) of the main cluster.

• •

In the second phase, the shock velocity deviates strongly from the velocity of the subcluster (in amplitude and/or direction), for example, the velocities differ when the subcluster is decelerated by the gravity of the main cluster after core passage, and finally falls backwards after crossing the apocenter (see the right panel in Fig. 1). The merger shock consequently transits into a “detached phase” – the bow shock detaches from the body which drives it (i.e. the gas core of the infalling subcluster) and moves farther towards the main-cluster outskirts. In this phase, the subcluster has only limited impact on the runaway shock properties, whose propagation is mainly determined by the shock Mach number and the density/temperature profiles of the main cluster.

In this paper, we concentrate on the fate of merger shocks in the “detached phase”. Unlike the case of a shock moving in a homogeneous medium, the steep gas density gradient of the main cluster atmosphere helps the runaway shock maintain its strength over a long distance (even possibly be amplified, if the gradient is sufficiently steep).

The current/recent generation of X-ray observatories (e.g. Chandra, XMM-Newton, Suzaku) have routinely measured the gas density profiles of galaxy clusters within $\sim R_{500}$ (the radius within which the mass density is 500 times the critical density at the cluster redshift) and sometimes beyond (see e.g. Vikhlinin et al. 2006; Bautz et al. 2009; George et al. 2009; Morandi et al. 2015; Ghirardini et al. 2019; Okabe et al. 2014). Observations of the Sunyaev-Zel’dovich (SZ) effect (e.g., Planck, SPT) extend the radial range to a few virial radii ($R_{\rm vir}\simeq 1.7R_{500}$). In general, the gas density profiles of galaxy clusters are fairly flat in the central region ($\sim r^{-1}$), and become progressively steeper with radius (e.g. $\sim r^{-2.5}$ around $R_{500}$). This implies that, beyond $\sim R_{500}$, a “Habitable zone” of runaway merger shocks exists, where shocks are able to maintain their strength and are, therefore, “long-lived”.

The aim of this work is to understand the behaviors of the runaway merger shocks in galaxy clusters (particularly in the “Habitable zones”). The radio relics discovered in cluster peripheries are possibly associated with these “long-lived” runaway shocks (see e.g. Vazza et al., 2012; Hong et al., 2014; Kang, 2015; Ha et al., 2018), which, on the other hand, provide a unique opportunity to identify and study these shocks whose X-ray signature has become too weak to be detected. In this study, instead of full cosmological runs, we use idealized/controlled simulations that elucidate conditions required for “survival” of the shocks. A one-dimensional (1D) blast wave, propagating through a power-law density distribution, is the simplest model that captures the essential properties of runaway shocks and probes a wide range of shock strengths and gas density slopes. Various models, combining shocks with the diffusive shock acceleration and/or adiabatic compression of pre-existing relativistic electrons have already been considered (see e.g. Enßlin et al. 1998; Roettiger et al. 1999; Enßlin & Gopal-Krishna 2001; Markevitch et al. 2005; Kang & Ryu 2015; see also Brunetti, & Jones 2014; van Weeren et al. 2019). Here we focus on a combination of a pure adiabatic compression (no re-acceleration) with a spherical runaway shock. We show that such shocks have a characteristic pattern of adiabatic compression/expansion, namely, the region of strong compression is confined to a relatively narrow shell on the downstream side. Donnert et al. (2016) have already considered a set of heuristic models where a parameter characterizing the expansion of the gas on the downstream side was introduced. We argue that this parameter is in fact linked to the Mach number of the shock and we model the evolution of the gas compression factor self-consistently.

This paper is organized as follows. In Section 2, we present qualitative analytic arguments relevant to shock propagation in a non-uniform medium. In Section 3, we show results of numerical simulations, which specifically aim to explore runaway shocks in the context of galaxy clusters and their synchrotron radiation. In Section 4, we discuss the existence of the “Habitable zone” of runaway shocks in the cluster outskirts. In Section 5, we summarize our findings.

2 Shock propagation in a non-uniform medium

As mentioned above, once the bow shock runs away from its driving subcluster (Fig. 2), the subsequent shock evolution is determined by its Mach number and the gas distribution in the atmosphere of the main cluster. During this phase, the motion of the subcluster has essentially no impact on the properties of the shock. Such runaway shocks are reminiscent of the spherical blast waves propagating in a non-uniform medium222The discussion in this section could also be applied to shocks/sound waves driven by outbursts of supermassive black holes in cluster center, where gas density profiles are shallower.. It is therefore useful to consider, how a spherical shock is attenuated/amplified when it travels in a stratified atmosphere. For simplicity, we assume a radial gas density profile of the medium that follows a power law, i.e. $\rho_{\rm gas,0}(r)\propto r^{-\omega}\ (\omega\geq 0)$333The subscript [math] in gas density $\rho_{\rm gas,0}$ indicates the unperturbed state of the atmosphere..

In general, the strength of the propagating spherical shock varies with distance due to three effects: (1) the non-uniform gas density profile, (2) the geometric factor ($\propto r^{2}$), and (3) the shock dissipation process. For linear waves, only the first two effects are relevant. The energy conservation law (i.e. $\rho_{\rm gas}(R_{\rm s})R_{\rm s}^{2}u_{\rm g}^{2}={\rm constant}$) shows that, the velocity amplitude $u_{\rm g}$ of a linear wave changes with radius $R_{\rm s}$ as a power law $u_{\rm g}\propto R_{\rm s}^{\eta_{\rm s}}$, with

[TABLE]

where $u_{\rm g}$ and $c_{\rm s}$ are the gas velocity in the rest frame of the upstream medium and the sound speed, respectively; $R_{\rm s}$ is the position of the wave (see green line in Fig. 2) or the shock (see discussions below).

Once the wave becomes nonlinear, a discontinuity (weak shock) is formed and dissipation contributes to the decline of the amplitude. Landau (1945) has shown that for $\omega=0$, the velocity amplitude of a very weak shock, far from its place of origin, declines with distance $R_{\rm s}$ as

[TABLE]

where $\alpha$ is a certain constant. Extending Landau’s argument for $\omega>0$, we find that the asymptotic behavior of a weak shock also follows a power law, whose exponent (see blue line in Fig. 2) is given as

[TABLE]

It is worth noting that, unlike linear waves, the characteristic length scale of weak shocks changes with distance as $R_{\rm s}^{\omega/4}$ (in the limit of very large distance). We thus limit our discussion here to $\omega\leq 4$, since for steeper profiles the behavior is more complicated (the shock length scale would become longer than the density scale height of the atmosphere $H_{\rho}\equiv|{\rm d}\ln\rho_{\rm gas,0}/{\rm d}r|^{-1}=r/\omega$ in a finite time).

It is not surprising that strong shocks behave differently from weak shocks. The former are described by a self-similar (SS) solution if the shock itself is sufficiently strong that the pre-shock pressure can be ignored. In strong shocks, the gas velocity behind the shock is $u_{\rm g}=3V_{\rm s}/4$, where $V_{\rm s}$ is the shock velocity (the adiabatic index $\Gamma=5/3$ is used in this work, unless stated otherwise). When $\omega<3$, a first-type SS solution exists (Sedov, 1959; Book, 1994), which gives

[TABLE]

see red solid line in Fig. 2. For $\omega<3$, the shock velocity decreases with distance. When $\omega$ approaches $3$, the shock velocity is expected to remain approximately constant as a function of radius. When $\omega\geq 3$, the situation is different (Eq. 4 is invalid in this limit), because the gas mass and energy of the atmosphere diverge at the origin ($r=0$). For this case, Waxman & Shvarts (1993) found a second-type SS solution444More accurately speaking, the second-type SS solution exists when $\omega_{\rm g}<\omega<\omega_{\rm c}$, where $\omega_{\rm g}$ and $\omega_{\rm c}$ are functions of the gas adiabatic index $\Gamma$. For $\Gamma=5/3$, Waxman & Shvarts (1993) showed that $\omega_{\rm g}\simeq 3.26$ and $\omega_{\rm c}\simeq 7.69$. When $3\leq\omega\leq\omega_{\rm g}$, the change rate $\eta_{\rm s}=0$ gives an asymptotic solution (Kushnir & Waxman, 2010)., where only the outer part of the flow follows a SS solution and $\eta_{\rm s}$ is found from the numerical solution (e.g., $u_{\rm g}\propto R_{\rm s}^{0.21}$ when $\omega=4$; see red dashed line in Fig. 2).

Fig. 2 shows the expected behavior of $\eta_{\rm s}$ as a function of $\omega$ for the three different cases discussed above (namely, a linear wave, a very weak shock, and a strong shock). It is clear that (i) all three curves monotonically increase with $\omega$ and (ii) the values of $\eta_{\rm s}$ for strong and weak shocks are smaller than those for linear waves.

3 Numerical results

The arguments presented in the previous section suggest that a shock, traveling radially, will be attenuated or amplified depending on the strength of the shock and on the slope of the density profile. However, application of these arguments to merger shocks in clusters faces two difficulties. First, merger shocks usually have moderate Mach numbers (e.g. $\mathcal{M}_{\rm s}\sim 2-4$). Neither the results for a very weak shock nor for a strong shock can be directly adopted. Second, weak and strong shocks both require a long period of time (or long travel distance) to evolve towards their asymptotic or SS solutions. This condition, however, might not be satisfied in galaxy clusters. For strong shocks, this issue is particularly serious when $\omega$ is around $3$ (see e.g. Kushnir & Waxman 2010). For these reasons, we perform numerical simulations555All simulations presented in this section are done by the mesh-based code OpenFOAM (Open Source Field Operation and Manipulation, www.openfoam.org). in this section to specifically model the motion of runaway shocks in galaxy clusters.

In this study, we applied a 1D blast-wave model, where the shock is initiated by an instantaneous energy release in a small region at the cluster center (see Section 3.1). In all our simulation runs, a static gravitational field with spherical symmetry is assumed to model the cluster environment. Since the radial temperature variations in galaxy clusters (a factor of $\sim 2-3$ within the virial radius) are much smaller than the corresponding gas density variations ($2-3$ orders of magnitude, see e.g. Ghirardini et al. 2019), the ICM is assumed to be initially isothermal ($T_{\rm gas,0}=1{\rm\,keV}$) and in hydrostatic equilibrium. The corresponding gas density profile follows

[TABLE]

where $\rho_{\rm c}=10^{7}{\,M_{\odot}}{\rm\,kpc}^{-3}$ and $r_{\rm core}=50{\rm\,kpc}$ are the central density and the core radius of the ICM, respectively. When $r\gg r_{\rm core}$, the logarithmic slope of the gas density profile ($\equiv{\rm d}\ln\rho_{\rm gas,0}(r)/{\rm d}\ln r$) approaches $-\omega$. We note that, the problem discussed here is almost scale-free (radiative cooling is ignored in this study). Our conclusions are only sensitive to the shock Mach number and the gas density slope, but not to the specific choice of parameters (e.g. $T_{\rm gas,0},\ \rho_{\rm c}$) used to describe the cluster atmosphere.

3.1 1D blast-wave model

Our simulation of a 1D blast-wave assumes a spherically symmetric cluster. Four groups of simulations with different initial gas density profiles ($\omega=1,\ 2,\ 3$, and $4$ in Eq. 5, respectively) are performed to study the role of the density profile slope (of the main cluster) on shock propagation. In each group of runs, blast waves are initiated by a sudden increase of pressure in the innermost cell ($r<1{\rm\,kpc}$) by a factor of $\xi\ (=10^{5},\,10^{6},\,10^{7})$ relative to the equilibrium value. Shocks with moderate Mach numbers $\mathcal{M}_{\rm s}\simeq 1-10$ are formed outside $r_{\rm core}$. The Mach numbers of the excited blast waves when they reach $R_{\rm s}=100$ and $950{\rm\,kpc}$ are summarized in Table 1.

Fig. 3 shows the evolution of the blast wave Mach numbers in runs with $\xi=10^{6}$ (solid lines). These results follow the trend discussed in Section 2 (see Fig. 2), namely, the shallower the gas density profile is, the more rapidly the shock fades with distance. We have further analyzed the rate of the wave amplitude change with radius ($\eta_{\rm s}$ introduced in Section 2). The values shown in Fig. 4 (symbols) broadly agree with theoretical expectations for the weak and strong shocks (see previous section), which are shown with blue and red lines. In our simulations, all shocks have faster attenuation rates during the initial period than at later times. This is mainly caused by the presence of a core ($r_{\rm core}=50{\rm\,kpc}$) in the adopted gas density profile (see Eq. 5). For example, at $r=100{\rm\,kpc}$, the logarithmic slope of the gas density profile is ${\rm d}\ln\rho_{\rm gas,0}(r)/{{\rm d}\ln r}=-2\omega/3$, i.e. flatter than the asymptotic value ($-\omega$) at large radii. This leads to the faster attenuation of the shock. When the shock leaves the core ($R_{\rm s}\gg r_{\rm core}$), the velocity attenuation rate approaches theoretical expectations (either for a weak shock or strong shock, depending on the Mach number of the shock). However, we do not expect an exact match between simulations and expectations, like those shown in Fig. 4. There are a number of reasons for the differences. First, in our simulations the shocks have moderate Mach numbers, which corresponds to neither very strong nor very weak shocks. Secondly, the discussion presented in Section 2 is based on the assumption that the shock has sufficient time to evolve towards the self-similar solution, which is not necessarily achieved in simulations999The specific condition for the SS solution in our problem is that the shock-swept gas mass is larger than the initial gas mass within the core of the atmosphere. This is also the reason why the numerical solution very slowly approaches the self-similar regime when $\omega\sim 3$, as in runs W3P6 and W3P7, shown in Fig. 4. The shocks still need to propagate a large distance to approach the red line (note that it is just a coincidence that their positions are close to the blue line right now).. Therefore, detailed simulations of each particular merger configuration are needed, if a quantitative answer is required. However, a qualitative conclusion, that steep density profiles help to prevent attenuation of shocks, of course, holds. Accordingly, runaway shocks are expected to be “long-lived” outside $R_{\rm 500}$ where the logarithmic slope of the gas density profile of the atmosphere (i.e. $-\omega$) is steeper than $\sim-3$. The above discussion assumes isothermal gas over the entire radial range, so that the velocity and the Mach number are unambiguously related. In real clusters, the temperature declines beyond $\sim R_{500}$ (e.g. Vikhlinin et al., 2005). This leads to an additional increase of the Mach number for shocks propagating to the outskirts of clusters (e.g. Ha et al., 2018).

Finally, we consider the shock structure in the blast-wave model, which is important for interpretation of the synchrotron emission (i.e. radio relics) associated with the merger shocks in galaxy clusters (see also Section 3.2). The top panel in Fig. 5 shows the evolution of the gas density profile in the run W3P6 with $\omega=3$. The profiles corresponding to later phases (see e.g. the red curve) are of particular interest here. As expected, the shock strength only weakly declines with radius, the gas profile has a characteristic jump at the shock and a smooth density distribution before and after the shock. The density jump, at the shock, depends on the adiabatic index and is, of course, a function of the shock Mach number, which can be derived from X-ray data. The question arises if the same density jump (or, equivalently, the compression factor) applies to the entire downstream region. Two versions of the compression factors are explored in this work (shown in the bottom panel in Fig. 5).

• •

The first version, that we call $C_{1}$, is the ratio of the gas density at a given radius to the initial density at the same radius, i.e., $C_{1}=\rho_{\rm gas}(r)/\rho_{\rm gas,0}(r)$. This ratio (shown by the dotted lines in Fig. 5) characterizes the modification of the density profile by the runaway shock. Since $C_{1}$ slowly decreases towards the cluster center, the density profile seen in the top panel appears slightly shallower than the initial density profiles. This justifies the use, by X-ray astronomers, of an analytic model consisting of two power laws with slightly different slopes for fitting the shocks (Markevitch & Vikhlinin, 2007).

• •

The second version of the compression factor, that we call $C_{2}$, characterizes the compression in the Lagrangian sense, namely the ratio of the gas density to the initial gas density of the same gas lump, i.e., $C_{2}=\rho_{\rm gas}(r)/\rho_{\rm gas,0}(r^{\ast})$, where $r^{\ast}$ can be found from the following integral equation:

[TABLE]

The $C_{2}$ profiles in the run W3P6 are shown in the bottom panel of Fig. 5 (solid lines). While $C_{2}\equiv C_{1}$ at the shock front, the width of region with large $C_{2}$ is clearly much smaller than that for $C_{1}$. Therefore, relative to the initial state of gas, the strong compression is in a narrow shell with width $\sim R_{\rm s}/5$ just behind the shock. At larger distances from the shock, the gas density of a given lump is even smaller than it was before being swept by the shock. Fig. 6 shows that this result is not sensitive to $\omega$ rather, but depends on the shock Mach number (saturates for $\mathcal{M}_{\rm s}>3$) and reflects the spherical geometry of the problem. Indeed, unlike the case of a plane-parallel shock, in spherical geometry, the gas on the downstream side is expected to expand. Donnert et al. (2016) modelled the structure of the shock in the CIZA J2242.8+5301 cluster with an ad hoc model, that uses a parameter $t_{\rm exp}$ to account for adiabatic expansion (see their section 2.1). In their model, $t_{\rm exp}$ is essentially a free parameter. Figs. 5 and 6 show that, in fact, the compression ratio is tightly linked to the Mach number, removing the uncertainty associated with the choice of $t_{\rm exp}$. Moreover, the compression ratio drops below unity for $M\gtrsim 2$ at $r\sim 0.8R_{\rm s}$, when the expansion is accounted for self-consistently. An approximate comparison with the Donnert et al. (2016) model shows that the size of the compressed region is a factor of $2-3$ smaller than the value expected for $t_{\rm exp}=240{\rm\,Myr}$ used for the CIZA J2242.8+5301 shock.

3.2 Synchrotron radiation

Even though our simulations are purely hydrodynamic, we can assess the expected impact of propagating shocks on the synchrotron emission under the assumption that relativistic particles and magnetic fields have no direct impact on the dynamics of the gas.

Since the shocks considered here have modest Mach numbers $\lesssim 2-3$, the particle acceleration efficiency might be low (e.g. Kang, 2015). We therefore assume that no acceleration takes place and focus exclusively on the adiabatic compression of pre-existing (fossil) relativistic electrons as the mechanism to boost the synchrotron emissivity in the downstream region (see e.g. van der Laan, 1962; Enßlin & Gopal-Krishna, 2001; Enßlin & Brüggen, 2002; Pfrommer & Jones, 2011). Furthermore, we assume that these fossil electrons are uniformly distributed in the unshocked gas, and initially have a power-law energy distribution, i.e. ${\rm d}n_{e}/{\rm d}\gamma\propto\gamma^{-p}$, where $p=3.5$ and $\gamma$ is the particle Lorentz factor. The relativistic electrons are assumed to be isotropic in the upstream and downstream regions.

We assume that the magnetic field, on the upstream side of the shock, is fixed at $B_{0}=2{\,\rm\mu G}$. For the downstream side, the magnetic field strength is set to $B(r)\propto B_{0}C_{2}(r)^{2/3}$, which corresponds to the case of relativistic gas with adiabatic index $\Gamma=4/3$ with compression ratio $C_{2}$.

While the above assumptions are clearly idealized, they serve the purpose of illustrating the magnitude of the synchrotron emissivity boost factor $W_{\nu}$ at frequency $\nu$ for a given gas lump, after it passes through the shock. The time evolution of the synchrotron emission subject to the adiabatic compression and inverse Compton (IC) and synchrotron energy losses, can be easily calculated, if the time evolution of the compression factor is known, which in our case is given by function $C_{2}(t)$ (e.g. Matthews & Scheuer, 1990; Churazov et al., 2001; Enßlin & Gopal-Krishna, 2001). This allows us to recover the time evolution of the boost factor $W_{\nu}$.

Fig. 7 shows the synchrotron emissivity boost factor profiles for the simulation run W3P6. By design, the scale factor is equal to unity in the upstream region. The top panel shows the time evolution of the $W_{\nu}$ profiles, when the cooling processes (synchrotron radiation and IC scattering) are not included. In the absence of cooling losses, the adiabatic compression preserves the shape of the radiation spectrum but only modifies its amplitude by a frequency-independent factor $W_{\nu}=C_{2}(r)^{1+2p/3}$ (Markevitch et al., 2005). This plot shows that, even without cooling losses, the strongly enhanced synchrotron emission is confined to a relatively narrow shell of compressed gas and drops in the downstream region due to gas re-expansion. In projection, such structures might appear as a bright rim with a sharper outer edge and somewhat smoother inner one.

If the effects of cooling losses are included, the profiles become even narrower. This is illustrated in the bottom panel of Fig. 7. For simplicity, we turn on the cooling losses only on the downstream side of the shock. Since higher energy electrons have shorter lifetimes and the characteristic frequency of the synchrotron emission scales as $\propto\gamma^{2}$, the width of the enhanced emissivity region decreases with the emitted frequency $\nu$.

Thus, the structure of the runaway merger shocks can naturally produce peaked emissivity profiles even when the cooling losses can be neglected. This naturally follows from a sequence of compression at the shock and successive expansion of the gas. The above discussion is of course a qualitative one. We defer the discussion of real clusters, the analysis of the magnetic field topology, and projection effect to our forthcoming paper.

4 “Habitable zone” of runaway merger shocks

In Sections 2 and 3, we have shown that the fate of runaway merger shocks, traveling through the ICM, is sensitive to the steepness of the gas density profile of the atmosphere. These shocks could be long-lived if the gas density profile is as steep as $\sim r^{-2}-r^{-4}$. It is therefore of considerable interest to discuss the possible “Habitable zone” of runaway merger shocks existing in galaxy clusters.

Fig. 8 shows the logarithmic slope of the gas density profile as a function of radius for 12 galaxy clusters (solid lines, presented in Vikhlinin et al. 2006), derived from the cluster best-fit model (eq. 3 in Vikhlinin et al. 2006, but the second (inner) $\beta$-model component is not included, because we focus here on the outer regions). $Chandra$ and $ROSAT$ data were used in the fitting. These curves are truncated at the detection radius $r_{\rm det}$ (see table 1 in Vikhlinin et al. 2006), far beyond which the extrapolation becomes unreliable. We can see that, within $0.5R_{500}$, the gas density profiles are generally shallower than $r^{-2}$. If the shock becomes detached from the subcluster in this region, then its amplitude is expected to decrease with radius as $\propto r^{-1/2}$. When $r$ approaches $R_{500}$, the profile becomes progressively steeper, with the logarithmic slope between $-2$ and $-3$. Beyond $R_{500}$, the slope continues to steepen for most of the cluster sample (see also Morandi et al., 2015). We stress here that cool, dense clumps, and filaments complicate accurate measurements of the X-ray profiles (Nagai & Lau, 2011; Simionescu et al., 2011; Zhuravleva et al., 2013; Avestruz et al., 2014) and might lead to a shallower profile than that of the dominant (by volume) component, which has lower gas density. This problem can be overcome with sensitive, high-angular resolution observations. To this end, we show the results of deeper Chandra observation of A133 with total exposure time $\rm\approx 2\,Ms$ (red curve with an error band, A. Vikhlinin, priv. comm.; see also Morandi & Cui 2014), which suggest that the slope steepens to values between $-3$ and $-4$, once all resolved bright clumps and filaments are excised. Cosmological simulations also suggest very steep slopes of the cluster gas density profiles ($\sim-3$ or steeper) around the virial radius (see e.g. Roncarelli et al., 2006; Vazza et al., 2010; Lau et al., 2015). The profiles might be even steeper around the cluster splashback radius (as steep as $-4$; see e.g. Diemer & Kravtsov 2014; Lau et al. 2015).

All these results suggest that a “Habitable zone” populated by moderately strong shocks could exist outside $\simeq R_{500}$, where runaway merger shocks do not attenuate quickly. While these shocks might be difficult to detect in X-rays, they might be responsible for powering radio relics often found in the periphery of merging clusters (see discussions in Section 3.2). In this regard, radio relics, discovered in cluster peripheries (e.g. A746 van Weeren et al. 2011; A2744 Pearce et al. 2017; see also van Weeren et al. 2019 for a review), provide a unique opportunity to locate and study runaway merger shocks, whose X-ray signals are too weak to be detected. Studying such shocks is important to our understanding the merger configurations of galaxy clusters. For example, the southwestern radio relics, discovered in the Coma cluster (Brown & Rudnick, 2011; Akamatsu et al., 2013; Erler et al., 2015), support the post-merger scenario proposed by Lyskova et al. (2019), where a runaway shock (consistent with the position of observed radio relics) is generated by the infalling galaxy group associated with NGC 4839.

Finally, we suggest that the differences in the gas profile slopes might provide a natural explanation for the “segmented” appearance of merger shocks discussed in Paul et al. (2011). We know that, not very far from the virial radius, the matter distributions in galaxy clusters are highly inhomogeneous/asymmetric. The gas density profile is shallower along the direction of filaments than between them (e.g. see Sato et al., 2012). Along the filaments, runaway merger shocks would propagate more slowly and attenuate more quickly, while between the filaments, the shocks would move more quickly and weaken more slowly. As a consequence, these runaway shocks would visually break into separate segments/wedges and appear more prominent in the “valleys” between the high-density filaments. The “surviving” segmented shocks might further interact with virial shocks and form more complicated structures in cluster outskirts (e.g. cold fronts, see Birnboim et al. 2010; see also e.g. Paul et al. 2011; Ha et al. 2018 for discussions in the context of cosmological simulations).

5 Conclusion

In this paper, we have explored the propagation of runaway merger shocks in galaxy clusters. Our findings can be summarized as follows:

• •

The evolution of merger shocks can be described as a two-phase process: the “driven phase” and the “detached phase”. In the detached phase, the subcluster, which drives a bow shock, is decelerated by the gravity of the main cluster, but the shock continues on its way towards the main-cluster outskirts. Meanwhile, the steep gas density profiles help the runaway shocks maintain their strength or even amplify them (if the profile is as steep as, or steeper than, $\sim r^{-3}$).

• •

High angular-resolution X-ray observations show that, beyond $R_{500}$, the gas density profile of galaxy clusters could indeed be as steep as $\sim r^{-3}$. This suggests that a “Habitable zone” could exist in the cluster outskirts, where runaway merger shocks are relatively “long-lived”. Since the gas density in this region is very low, the detection of these shocks in X-ray is challenging.

• •

A characteristic feature of runaway shocks with $\mathcal{M}_{\rm s}\gtrsim 2$ is that the strongly compressed gas, relative to its initial state, is confined to a relatively narrow region just behind the shock (similar to the case of a spherical shock moving in a uniform medium). Namely, a strong compression at the shock front is followed by a strong decompression. If adiabatic compression plays a role in boosting the radio emissivity of radio relics, then even in the absence of synchrotron and IC losses, the radio emission is expected to be confined to this shell.

Acknowledgements

We are grateful to Alexey Vikhlinin for providing us the A133 density profile. CZ thanks Daisuke Nagai for helpful discussions and comments on the draft. EC and NL acknowledge support from the Russian Science Foundation grant 19-12-00369. WF acknowledges support from the Smithsonian Institution, NASA Contract NAS-08060, and NASA Grant GO2-13005X.