Shock Waves in Cluster Outskirts

# Shock Waves and Cosmic Ray Acceleration in the Outskirts of Galaxy Clusters

## Abstract

The outskirts of galaxy clusters are continuously disturbed by mergers and gas infall along filaments, which in turn induce turbulent flow motions and shock waves. We examine the properties of shocks that form within in sample galaxy clusters from structure formation simulations. While most of these shocks are weak and inefficient accelerators of cosmic rays (CRs), there are a number of strong, energetic shocks which can produce large amounts of CR protons via diffusive shock acceleration. We show that the energetic shocks reside mostly in the outskirts and a substantial fraction of them are induced by infall of the warm-hot intergalactic medium from filaments. As a result, the radial profile of the CR pressure in the intracluster medium is expected to be broad, dropping off more slowly than that of the gas pressure, and might be even temporarily inverted, peaking in the outskirts. The volume-integrated momentum spectrum of CR protons inside has the power-law slope of , indicating that the average Mach number of the shocks of main CR production is in the range of . We suggest that some radio relics with relatively flat radio spectrum could be explained by primary electrons accelerated by energetic infall shocks with induced in the cluster outskirts.

acceleration of particles — cosmic rays — galaxies: clusters: intracluster medium — methods: numerical — shock waves
3

## 1 Introduction

Clusters of galaxies are the largest gravitationally bound structures that emerged from hierarchical clustering during the large-scale structure (LSS) formation of the Universe. While the central part of many clusters looks relaxed into hydrostatic equilibrium, especially in X-ray observations (e.g., Markevitch et al., 1998; Vikhlinin et al., 2006), the outskirts around the virial radius, , are stirred by mergers of substructures and continuous infall of gas along adjacent filaments (e.g., Ryu et al., 2003; Skillman et al., 2008; Vazza et al., 2009a). Observational evidence for the deviation from equilibrium in the cluster outskirts can be seen in the entropy distribution. The radial profile of obtained in X-ray observations follows roughly in the inner part of clusters, but beyond it flattens off and turns down (Voit et al., 2005; George et al, 2009; Walker et al., 2012; Simionescu et al., 2013). Moreover, according to structure formation simulations, turbulent flow motions develop during the formation of clusters; the ratio of turbulence to gas pressure increases outwards and reaches order of 10 % in the outskirts of simulated clusters (e.g., Ryu et al., 2008; Vazza et al., 2009b; Lau et al., 2009).

Although flow motions are expected to be on average subsonic in the cluster outskirts4, shock waves have been observed in X-ray as well as in radio. In X-ray observations, some of sharp discontinuities in the surface brightness are attributed to shocks, while others are attributed to cold fronts or contact discontinuities. The physical properties of these shocks including the sonic Mach number, , can be determined using the deprojected temperature and density jumps (Markevitch & Vikhlinin, 2007). Since a shock was found in the so-called bullet cluster (1E 0657â56) (Markevitch et al., 2002), about a dozen of shocks have been detected with Chandra, XMM-Newton, and recently Suzaku (e.g., Russell et al., 2010; Akamatsu et al., 2012; Ogrean & Brüggen, 2013). The shocks identified so far in X-ray observations are mostly weak with .

Shocks have been identified also in radio observations through the so-called radio relics (e.g., Feretti et al., 2012; Brüggen et al., 2012, for reviews). Radio relics are the radio structures of megaparsec size observed within the virial radius. They often show elongated morphologies with sharp edges in one side, and occasionally come in pairs located in opposite sides of clusters. Radio emissions from these structures usually exhibit high polarization fractions. Radio relics are interpreted as shocks, where relativistic electrons emitting synchrotron radiation are accelerated or re-accelerated. The properties of radio relic shocks such as , magnetic field strength, and the age can be estimated from the spectral index and spatial profile of synchrotron emissions (van Weeren et al., 2010; Kang et al., 2012). So far several dozens of radio relics have been observed, and the Mach numbers of associated shocks are typically in the range of (e.g., Clarke & Enßlin, 2006; Bonafede et al., 2009; van Weeren et al., 2010, 2012).

In a few cases, shocks were detected both in X-ray and radio observations. Interestingly, however, the shock characteristics derived from X-ray observations are not always consistent with those from radio observations. For instance, the shock in the so-called sausage relic in CIZA J2242.8+5301 was estimated to have in the analysis of radio spectrum based on diffusive shock acceleration (DSA) model (van Weeren et al., 2010), but X-ray observations indicated (Akamatsu & Kawahara, 2013). And the shock in the so-called toothbrush relic in 1RXS J0603.3+4214 has according to the radio spectral analysis, but according to X-ray observations (van Weeren et al., 2012; Ogrean et al., 2013). In addition, the positions of shocks identified in radio are often spatially shifted from those found in X-ray (see the references above). Resolving these puzzles would require further observations as well as theoretical understandings of weak collisionless shocks in the intracluster medium (ICM), which is a high beta plasma with (e.g., Ryu et al., 2008). Here, and are the gas thermal and magnetic pressures, respectively.

With relatively low Mach numbers as well as elongated morphologies and occasional parings in opposite sides of clusters, shocks observed in the outskirts are often considered to be induced by mergers. The hypothesis of merger shocks was explored in simulated clusters, especially for the origin of radio relics (Skillman et al., 2011; Nuza et al., 2012; Skillman et al., 2013). In these studies, shocks in clusters are identified and the injection and acceleration of cosmic-ray (CR) electrons are modeled. Then, along with a model for the magnetic field in the intergalactic medium (IGM), synthetic radio maps are produced and examined. These studies suggested that merger shocks with sufficient kinetic energy flux are likely to be responsible for observed radio relics. However, it was also argued that typical mergers are expected to induce mostly weak shocks with and major mergers with similar masses, which are required to explain, for instance, the sausage relic (van Weeren et al., 2010), tend to generate very weak shocks with (e.g., Gabici & Blasi, 2003).

The nature and origin of cosmological shocks have been studied extensively, using numerical simulations for the LSS formation of the Universe (Minati et al., 2000; Ryu et al., 2003; Pfrommer et al., 2006; Kang et al., 2007; Skillman et al., 2008; Hoeft et al., 2008; Vazza et al., 2009a; Brüggen et al., 2012). Shocks are induced as a consequence of hierarchical clustering of nonlinear structures and can be classified into two categories. External shocks form around clusters and filaments of galaxies, when the cool ( K), tenuous gas in voids accretes onto them. So the Mach number of external shocks can be very high, reaching up to or so. Internal shocks, which form inside nonlinear structures, on the other hand, have lower Mach numbers of or so, because they form in much hotter gas that was previously shocked. It was shown that while a large fraction of internal shocks have , those with are most important in dissipating the shock kinetic energy into heat in the ICM. Internal shocks are induced by mergers of substructures, as well as by turbulent flow motions and by infall of warm gas from filaments to clusters (e.g., Ryu et al., 2003).

Turbulent shocks, induced by turbulent flow motions, are expect to be weak with at most , because the root-mean-square (rms) flow motions are subsonic. Inflall shocks5, on the contrary, can have Mach numbers as large as , since they form by the infall of the so-called warm-hot intergalactic medium (WHIM) with K to the hot ICM with K. We note that a continuous infall containing density clumps would be difficult to be differentiated from a stream of minor mergers with small mass ratios. But infall shocks clearly differ from merger shocks, which are generated by major mergers, in the sense that they do not appear as a pair in opposite sides of clusters. In addition, infall shocks should be found mostly in the cluster outskirts, since the gas infall from filaments normally stops around the virial radius and does not penetrate into the core. So it would be reasonable to conjecture that while weak shocks with in clusters are mostly merger or turbulent shocks, stronger shocks with found in the cluster outskirts are likely to be infall shocks.

Shocks that can be categorized as infall shocks were identified in observations before. For instance, the radio relic 1253+275 in the Coma cluster was interpreted as an infall shock formed by a group of galaxies along with the intra-group medium accreting into the ICM (Brown & Rudnick, 2011). Also the radio structure of NGC 1265 in the Perseus cluster was modeled as the passage of the galaxy through a shock with formed by the infalling WHIM (Pfrommer & Jones, 2011). However, the properties such as the frequency, spatial distribution, and energetics of infall shocks have not been studied in simulations before, partly because the automated distinction of infall shocks from merger shocks or other types of shocks in simulated clusters is not trivial.

It is well established that CRs are produced via DSA process at collisionless shocks, such as interplanetary shocks, supernova remnant shocks, and shocks in clusters (Bell, 1978; Blandford & Ostriker, 1978; Drury, 1983). Shocks in the LSS of the universe are the primary means through which the gravitational energy released during the structure formation is dissipated into the gas entropy, turbulence, magnetic field, and CR particles (e.g., Ryu et al., 2008, 2012). Post-processing estimations with simulation data for the amount of CR protons produced in clusters showed that the CR pressure in the ICM may reach up to a few % of the gas thermal pressure (Ryu et al., 2003; Kang et al., 2007; Skillman et al., 2008; Vazza et al., 2009a). Observationally, on the other hand, the CR-to-thermal pressure ratio in clusters was constrained to be less than a few %, with the upper limits on -ray fluxes set by Fermi-LAT and VERITAS (Ackermann et al., 2010, 2013; Arlen et al., 2012).

In some simulations for the LSS formation, the injection/acceleration of CR protons at shocks and the spatial advection of the CR pressure were followed self-consistently in run-time (e.g., Miniati et al., 2001; Pfrommer et al., 2007; Vazza et al., 2012). Pfrommer et al. (2007) and Vazza et al. (2012), adopting specific DSA efficiency models, showed that the CR acceleration occurs mostly in the cluster outskirts. Because of long lifetime and slow particle diffusion (e.g., Berezinsky et al., 1997), the CR protons accelerated in the outskirts are likely to be contained in clusters and accumulated in the ICM over cosmological time scales. But they can be advected with turbulent flows toward the central part of the cluster (see, e.g., Enßlin et al., 2011). For simplicity, let us assume that the transport of CR protons due to flow motions can be approximated by turbulent diffusion, then it could be be described by , where is the density of CR protons and is the turbulent diffusion coefficient. If only the radial diffusion is considered and the diffusion coefficient is approximated as , where is the average flow speed at , then the advection time scale can be estimated rather roughly as . In the cluster outskirts, typically and a few km/s, so a few yrs. This is a substantial fraction of the age of the universe, implying that it would take a while for CR protons produced at energetic shocks in the outskirts to reach the core region. As a result, the radial distribution of the CR pressure would be broad, dropping off more slowly than that of the gas thermal pressure in the outskirts. Vazza et al. (2012) also showed that the CR pressure distribution could be temporarily inverted, that is, the CR pressure can increase outwards.

Brunetti et al. (2012), on the other hand, attempted to constrain the radial distributions of nonthermal components (including the CR proton energy density) in the Coma cluster by combining radio observations with recent Fermi-LAT -ray observations and with Faraday rotation measure (RM) data. They argued that the model based on the turbulent acceleration of secondary electrons would best reproduce the radio halo of the Coma cluster with the CR energy density that scales with the thermal energy density as with to , implying that is higher at lower . The outer region of the Coma cluster is strongly disturbed by ongoing mergers and infalls (e.g., Simionescu et al., 2013), so it would be probably one of rare cases with this kind of inversion of the profile. But these indicate that the partitioning of thermal and CR energies (and possibly turbulent and magnetic field energies too) could be very different in different parts of clusters.

In this paper, we study shocks within the virial radius in a sample of clusters taken from LSS formation simulations. Specifically, we examine the properties of energetic shocks with relatively high Mach number and high shock kinetic energy flux that can produce large amounts of CR protons via DSA. The plan of this paper is as follows. In Section 2, numerical details are presented. In Section 3, the properties and nature of shocks in the cluster outskirts are described. In Section 4, the properties of CRs produced at energetic shocks are described. Discussion is given in Section 5, and summary follows in Section 5.

## 2 Numerical Details

### 2.1 Cluster Sample

To produce a sample of galaxy clusters, we performed simulations of the LSS formation, using a Particle-Mesh/Eulerian cosmological hydrodynamics code described in Ryu et al. (1993). A standard CDM cosmological model was assumed with the following parameters: baryon density , dark matter density , cosmological constant , Hubble parameter , rms density fluctuation , and primordial spectral index . These parameters are consistent with the WMAP7 data (Komatsu et al., 2011). Cubic boxes of comoving sizes of and with periodic boundaries were employed and divided into grid zones with uniform spatial resolutions of and , respectively. Nongravitational processes such as radiative cooling, star formation and feedback, and reionization of the IGM were not considered. Instead, a temperature floor was set to be K for the gas in voids, assuming that the unshocked gas outside nonlinear structures was uniformly heated by reionization. To compensate the cosmic variance and acquire an enough number of clusters, 16 runs with different realizations of initial condition were performed for each of and boxes (so the total number of runs is 32).

In addition, we used a higher resolution simulation with grid zones in box of comoving size (), to mainly examine the resolution effects. This simulation was performed with a numerical code described in Li et al. (2008), adopting the same set of cosmological parameters except . It is basically the same simulation reported in Cen & Chisari (2011), but with the box size of instead of . The simulation includes a mild feedback from star formation (low galactic superwind feedback of Cen & Chisari (2011)) and cooling/heating processes. Kang et al. (2007) examined the effects of a similar feedback and radiative processes on the properties of shocks in the LSS. They showed that the dynamics and energetics of shocks are governed primarily by the gravity of matters, so mild feedback and cooling do not significantly affect the statistics of the shocks in the ICM (see Pfrommer et al., 2007, for the case that the feedback is stronger and its effects are more important).

In the simulation data, we identified clusters as the volumes with high X-ray luminosity (see Kang et al., 1994, for details). For each identified cluster, we calculated the gas mass, , and the X-ray emission-weighted average temperature, , inside from the cluster center that locates at the peak of X-ray emissivity. Here is defined as the radius within which the gas overdensity is 200 times the mean gas density (not the critical density) of the universe.6 We built our sample with clusters of keV from box simulations and those of keV from box simulations, by optimizing the resolution limitation and the size of cluster sample; 125 clusters were identified from 16 simulations of box with zones, 94 clusters from 16 simulations of box with zones, and 9 clusters from one simulation of box with zones. Figure 1 shows the radius-mass relation and the radius-temperature relation of the total 228 clusters in our sample. The simulated clusters have , , and keV. From the virial theorem, the mass, temperature, and radius of relaxed clusters are expected to follow and or (e.g., Peebles, 1980). As can be seen from Figure 1, overall , , and of our clusters follow these relations, but there are scatters, a part of which are caused by dynamical activities in the cluster outskirts.

We point that with uniform grids of , our simulated clusters have poorer resolution than those generated using adaptive mesh refinement (AMR) or SPH codes (see Introduction for references), especially in the core regions. This means that shocks in the inner regions of clusters may not be fully reproduced. But those shocks in high density regions are expected to be weak with low Mach numbers of (e.g., Vazza et al., 2011), so their contribution to the production of CRs would not be significant (see the discussion in Section 2.3). Since this paper focuses on relatively strong, CR-producing shocks in the outskirts, having an uniform resolution throughout the entire simulation volume should be actually an advantage.

### 2.2 Shock Identification

A number of algorithms that can be applied to the identification of shocks in structure formation simulation data have been suggested (e.g., Ryu et al., 2003; Pfrommer et al., 2006; Skillman et al., 2008; Vazza et al., 2009a). They all employed the Rankine-Hugoniot jump conditions but in slightly different ways. Although there are some differences, the properties of identified shocks in the LSS by different algorithms are overall consistent with each other (e.g., Vazza et al., 2011, for a comparison study). Here we adopted the algorithm suggested by Ryu et al. (2003).

A series of the following one-dimensional procedures are first gone through for three primary directions. The grid zones are tagged as ‘shocked’ if they fulfill all of the following conditions: (1) , i.e., the local flow is converging, (2) , i.e., the Mach number is greater than 1.3, and (3) , i.e., the gradients of temperature and density have the same sign. The central difference is defined as for the quantity in the zone . A shock in simulation usually spreads over several grid zones, and the ‘shock center’ is defined as the grid zone with minimum . While Ryu et al. (2003) used the condition (where is the entropy), here we used the condition in order to exclude the possible misidentification of contact discontinuities. We found that the current method with the condition may miss some of weak shocks, but the overall statistics of identified shocks are not significantly affected. The Mach number at the shock center, , is calculated by solving the relation for the gas temperature jump along the three primary directions: . Hereafter, the subscripts 1 and 2 indicate the preshock and postshock quantities, respectively. Then the Mach number of the shock center is assigned as the maximum value of the three Mach numbers, i.e., . Because of complex flow patterns and shock surface topologies, very weak shocks are difficult to be identified unequivocally, so only shocks with are considered. Hereafter, we refer to a grid zone with assigned as “a shock”, which represent a small patch with an area of . A shock surface normally consists of many of these shocks (or shock zones), so the total number of identified shocks multiplied by is effectively equal to the total area of shock surfaces contained in a given volume.

Once the Mach number is determined, the shock speed and the shock kinetic flux are estimated as and , where is the gas adiabatic index.

### 2.3 Energy Dissipation at Shocks

If no CRs are produced at a shock, the gas thermalization efficiency can be calculated directly from the Rankine-Hugoniot relation as , where is the internal energy density. Note that the second term inside the brackets subtracts the effect of adiabatic compression that occurs at a shock as well as the contribution of the thermal energy flux entering the shock. Then the generation of heat can be estimated with the thermal energy flux, . However, if CRs are accelerated via DSA, a fraction of the shock kinetic energy is transferred to the CR component and the resulting thermalization efficiency is reduced, i.e., . With defined as the CR-acceleration efficiency (e.g., Ryu et al., 2003; Kang & Jones, 2007), the acceleration of CR protons at shocks can be quantified with the CR energy flux, .

At the moment it is not possible to predict and from first principles, because complex wave-particle plasma interactions governing the CR injection and acceleration at collisonless shocks are not fully understood. It has been recognized that the magnetic field amplification (MFA) due to CR streaming instabilities and the Alfvénic drift of scattering centers in the amplified field play significant roles in DSA at astrophysical shocks such as supernova remnant shocks (e.g., Lucek & Bell, 2000; Bell, 2004; Schure et al., 2012; Caprioli, 2012; Kang, 2012). Through numerical simulations of nonlinear DSA for shocks expected in the LSS, Kang & Ryu (2013) has shown that if self-amplification of magnetic fields and fast Alfvénic drift in the shock precursor are implemented into the standard DSA theory, the CR energy spectrum is steepened and the CR-acceleration efficiency is reduced, compared to the cases without including those processes. Here, we adopted and of Kang & Ryu (2013).

Figure 2 shows the curves that fit the values of and for shocks that form in a weakly magnetized plasma with and . The high value of is expected for plasmas in the ICM, as noted in Introduction. Both and increase as increases, and asymptote to and , respectively, for strong shocks with . Compared to the previous estimate of at strong shocks, given in Kang et al. (2007) where MFA and Alfvenic drift were not considered, the newly estimated CR-acceleration efficiency is smaller by a factor of for . Also the CR-acceleration is inefficient at weak shocks with in our new estimation.

The steepening of CR spectrum due to Alfvénic drift and the ensuing reduction of become important, only if the magnetic field is strong enough so that the Alfvénic speed is substantial, that is, (), where is the Alfvénic speed in the amplified magnetic field in the shock precursor (Caprioli, 2012; Kang, 2012). At weak shocks (), however, MFA would be inefficient, so , where is the Alfvénic speed in the background magnetic field. For , then , and so the Alfvénic drift effect would be only mildly important at weak shocks. At strong shocks, on the other hand, the diffusive CR pressure induces a precursor, in which the upstream flow is decelerated and adiabatically compressed, and the streaming CRs amplify significantly the turbulent magnetic fields (Bell, 2004). According to the MFA prescription adopted in Kang & Ryu (2013), the MFA factor increases with and can be approximated as , where and are the magnetic field strengths in the background medium and immediately upstream of the shock, respectively. Then, the Alfvénic Mach number defined by the amplified magnetic field becomes (independent of the plasma ), so the Alfvénic drift is expected to be important at strong CR modified shocks. We note, however, that relevant plasma physical processes, such as the injection of CRs as well as MFA and Alfvénic drift, are not well understood, so any attempts to predict the DSA efficiency involve large uncertainties, especially for weak shocks. So the dissipation efficiencies given in Figure 2 should be taken as rough estimates.

In this paper, we do not consider the re-acceleration of pre-existing CRs. We note, however, that it could be important especially for weak shocks with (Kang & Ryu, 2011, 2013). Kang et al. (2012) and Pinzke et al. (2013), for instance, argued that the re-acceleration of pre-existing CR electrons would be operating at radio relics associated with weak structure formation shocks.

## 3 Properties and Nature of Shocks in Cluster Outskirts

We first examine the Mach number and energetics of shocks within and around the virial radius, specifically in , of simulated clusters. Figure 3 shows the frequency distribution of shocks (i.e., zones with assigned ) as a function of Mach number and CR energy flux for the shocks found in 134 clusters from box simulations with and grid zones. The energetics of shocks is quantified with the CR energy flux, . Most of these shocks are internal shocks which are produced in the hot ICM of clusters or in the WHIM of filaments, according to the classification of Ryu et al. (2003). As previously shown, they are mostly weak with . The fractions of shocks with relatively high Mach numbers of , , and are , , and , respectively, among all the shocks with in and box simulations. We categorize shocks or shock zones with as energetic shocks. We note that the shocks responsible for observed radio relics are estimated to have the total kinetic energy flux of over the entire shock surface of or so (e.g., van Weeren et al., 2010; Kang et al., 2012). Considering at , shocks with can be considered to be energetic enough to be parts of observable radio relics. The fraction of energetic shocks is among all the shocks. And the fraction of shocks with and is .

Shocks with lower form on average close to the core with higher gas density (see Figure 4 below) and so have larger , but they have lower . Stronger shocks, on the other hand, form mostly in the outskirts and have lower , but they have higher . Such tendencies are reflected in the relation between and in Figure 3. For weak shocks with , the CR acceleration efficiency increases steeply with , while the shock kinetic energy flux varies only mildly. So increases strongly with , resulting in a relatively robust correlation between the two quantities. For shocks with , the dependence of on becomes much softer, while the variance of increases. So the correlation between and substantially weakens. We find that shocks with the largest have typically , which interestingly coincides with the Mach numbers of strong radio relic shocks (see Introduction).

Figure 4 shows two-dimensional slices of three sample clusters with the X-ray emission-weighted temperatures of keV (left), keV (middle) and keV (right), respectively, at present (). The slices were chosen to highlight the shock structures, so they pass through short comoving distances of from the cluster centers. The CR luminosity, shown in the bottom panels, is at the comoving surfaces of shocks. Hereafter the shock with the largest among shocks in each cluster will be called the most energetic shock (MES). Thick arrows in the and panels of Figure 4 point the MESs of each cluster. The MESs are located at (left), (middle) and (right) from the cluster centers, which correspond to (left), (middle) and (right), respectively. The spatial distribution of tells that strong shocks are found in the outer regions of clusters. Obviously, the strongest shocks are external shocks that form in the accreting gas from voids (K) (see Ryu et al., 2003). But as noted in Introduction, owing to low density, they are energetically unimportant, so we are not concerned about those external shock in this paper. Energetic shocks that produce large amounts of CRs are internal shocks and they reside mostly in the outskirts of the clusters, as shown in the distribution of .

The distributions of , , , and indicate that the structures including the MESs in Figure 4 look like infall shocks that form by the infall of gas from filaments. Those infall shocks are energetic enough to penetrate into the region inside the virial radius where the gas density is relatively high, indicating that their shock kinetic energy flux is large. They are also relatively strong with , so they are efficient CR accelerators. These characteristics make the infall shocks the MESs in the clusters shown here. We point that not all filaments induce infall shocks inside the virialized regions of clusters. Also the cross sectional areas of penetrated filaments are small, compared to the surface area of virialized regions . So the energetic infall shocks inside the virial radius should account for a small fraction of internal shocks in the ICM. The distribution of in Figure 4, however, indicates that infall shocks in the outskirts could be responsible for a significant fraction of CR production in the clusters (see below).

Figure 5 shows the time evolution of the cluster shown in the right column of Figure 4, demonstrating how the cluster has evolved dynamically and how various types of shocks have been induced. Strong external shocks persist at the comoving distance away from the cluster center, while numerous internal shocks form and disappear in a dynamical timescale of Gyr. One can see that the cluster experienced a merger at the look-back time of years, producing several merger shocks. Then infall from an attached filament in the south-west direction followed, and penetrated into the core region. It was halted by an energetic infall shock that developed around the look-back time of 1 Gyr and lasted to the present epoch.

To quantify the statistics of infall shocks, we separated them from other shocks. Infall shocks are defined as those that decelerate the WHIM accreting from filaments to a cluster (see Introduction). So we employed the following criteria for infall shocks in using the entropy and density of the preshock gas and the sonic Mach number: (1) , (2) , and (3) . Note that the first criterion is an entropy condition. Figure 6 shows the volumetric distribution of the gas from box simulations with grid zones. The domain demarcated by the first and second criteria does not coincide with the conventional definition of the WHIM, K K. But visual inspections indicated that the above three criteria pick up infall shocks in the region of best among different criteria we have tried. Figure 7 shows a slice displaying the positions of shocks identified as infall shocks according to the above criteria as well as those which are not infall shocks, for the cluster shown in the right column of Figure 4. We note that the morphology of shock surfaces could be quite complicated, depending on the dynamical history of clusters. In the cluster shown, for instance, the infall shock including the MES is connected to a larger shock surface, a portion of which is a (not-infall) shock expanding from the core to the outskirt. In general a connected shock surface can consist of a number of infall and not-infall shocks.

With the above criteria, we found that, among all the shocks with located in of the sample clusters, about are classified as infall shocks. So most of ICM shocks are merger shocks or turbulent shocks (see Introduction). As noted in Figure 3, of the identified shocks have , so about a half () of them are infall shocks. Identifying merger shocks would require the examination of the time evolution of cluster dynamics, which we did not attempt here.

Among the energetic shocks with , the fraction of infall shocks is . And among the MESs in 228 sample clusters, 177 are infall shocks; i.e., of the MESs are infall shocks in our sample. We expect the MESs that are not infall shocks to be merger shocks. The top panels of Figure 8 show the distributions of the radial position, , and the Mach number, , of the MESs. The MESs that are infall shocks (red histogram, MES-ISs hereafter) are distributed over all radius, peaking at . On the other hand, the MESs that are not infall shocks (blue histogram, MES-NISs hereafter) are mostly found at inner parts of clusters. For the MES-ISs, the Mach number distribution peaks around , and decreases sharply for smaller but extends to larger . For the MES-NISs, the Mach number distribution is mostly limited to . The figure demonstrates that the MES-ISs are found mainly at outer parts of clusters and they have higher Mach numbers than the MES-NISs, as expected. We attempt to find correlations among , , and in the bottom panels of Figure 8. It appears that there is no noticeable correlation between and . But tends to be larger at larger , confirming that higher Mach number shocks form at outer regions of clusters.

## 4 Cosmic Ray Production at Shocks in Cluster Outskirts

We next examine the spatial characteristics of the CR proton production in clusters. Figure 9 shows the radial distributions of the maximum Mach number, , the CR luminosity per unit radius, , and the fraction of CR luminosity due to infall shocks, , in four sample clusters. Here, is defined as the highest Mach number of shocks located in the bin of , while is the sum of of shocks in the bin divided by the width . Note that represents the amount of CR energy produced per unit time per unit length and has a rather unusual unit, 7. The vertical dashed lines mark the radial bins that contain the MESs.

The distribution of demonstrates that on average shocks tend to be stronger in the outer regions of clusters, as expected. The CR luminosity is dominated by energetic shocks in a given radial bin, and energetic shocks are found mostly in the outskirts, so is higher there. The MESs and the peaks of are located in (although the two do not necessarily coincide), indicating that CRs are produced mostly at the outer regions of clusters. These findings are consistent with the previous studies using the structure formation simulations in which the production of CR protons was explicitly followed in runtime (Pfrommer et al., 2007; Vazza et al., 2012). With high CR production at the outskirts, we expect the the radial profile of the CR pressure is flatter than that of the gas pressure and could be even inverted (see Brunetti et al., 2012), as mentioned in the Introduction.

The distribution of shows that infall shocks contribute to the CR production by a large fraction, especially in the outskirts. As a matter of fact, we estimate that infall shocks produce of CRs in , when summed for all 228 clusters in our sample. Recall that the fraction of infall shocks among those with inside the volume of is .

We also examine the momentum distribution of the CR protons, which are expected to be produced in the outskirts and then mixed in the ICM via turbulent flow motions (see Introduction). For each shock with , we adopted the test-particle power-law distribution function, , where (Drury, 1983). Note that the CR acceleration at cluster shocks with is reasonably described by the test-particle solution (Kang & Ryu, 2010). In our sample, of shocks have and shocks with the largest have typically (see Section 3). So the test-particle solution should give reasonable results. For each sample cluster, the volume-integrated, momentum distribution function of CR protons within the radius , , is estimated as follows: (1) At each shock zone , the power-law function, , is normalized with the CR energy flux at the shock as . Here, is expressed in unit of , and and are assumed. 2) The volume-integrated momentum distribution function is calculated by adding up for all shocks inside , i.e., . 3) Then, the slope of the integrated momentum distribution function is estimated by fitting to a power-law form with the slope .

The left panel of Figure 10 shows the values of calculated for some of sample clusters8. The average value of decreases with , reflecting the fact that shocks are stronger on average in the outer regions of clusters, and the variation of its distribution also decreases with . The average values, at and at , correspond to the DSA power-law slopes for shocks with and 3.5, respectively. At , the values spread over a narrow range of , corresponding to , as shown in the right panel of Figure 10. The plot indicates that does not have any noticeable correlation with the cluster temperature (nor with the cluster mass and radius although not shown). These imply that the averaged Mach number of shocks, weighted with CR production, in our sample clusters is in the range of , regardless of the properties of clusters.

## 5 Discussion

The existence of CR protons in galaxy clusters remains to be confirmed directly from observations. CR protons produce -ray radiation through collisions with background thermal protons, but so far only upper limits on -ray fluxes from clusters have been set by Fermi-LAT and VERITAS, as noted in the Introduction (Ackermann et al., 2010, 2013; Arlen et al., 2012). On the other hand, secondary electrons are also produced through collisions, and they along with G-level magnetic fields emit the synchrotron radiation in radio over the cluster scale. The observed radio emission is produced typically by electrons with energy of several GeV corresponding to the Lorentz factor of (Kang et al., 2012). For this energy range, the secondary electrons at production have the momentum distribution similar to the proton spectrum, that is, with (e.g., Dermer, 1986). For the proton power-law of estimated in the previous section, the secondary electron slope is also . Relativistic electrons suffers radiative coolings, dominantly by synchrotron and inverse Compton losses. The cooling time scale of electrons of is yrs with cluster magnetic fields of a few G (e.g., Kang et al., 2012). The momentum distribution function of the secondary electrons at the steady-state governed by the production through collisions and the cooling processes is given as with (Dolag & Enßlin, 2000). So for , the secondary electron slope becomes .

So far, our discussions on DSA at energetic shocks have been focused mostly on CR protons and secondary electrons resulted from collisions. As for supernova remnant shocks, primary CR electrons can be accelerated at ICM shocks in the same manner as CR protons, although the injection (pre-acceleration) of electrons into DSA process remains rather uncertain. We point that the projected surfaces of energetic shocks would have morphologies of partial shells or elongated arcs (see Figure 4), so diffuse synchrotron emissions from primary CR electrons accelerated at these shocks could produce radio structures that resemble radio relics discovered in the cluster outskirts (e.g. van Weeren et al., 2010; Kang et al., 2012). Moreover, the average radial distance of the MESs in Figure 8 is , which is comparable to the average distance of observed radio relics from the cluster center (e.g., van Weeren et al., 2009). So, for instance, radio relics with flat radio spectrum such as the sausage relic in CIZA J2242.8+5301 () could be explained by primary electrons accelerated by energetic shocks (a substantial fraction of which are infall shocks) in the cluster outskirts.

In addition, pre-existing CR electrons in the ICM (previously produced at shocks or through collisions) can be re-accelerated at energetic shocks (Kang et al., 2012; Pinzke et al., 2013). To explain the observed properties of radio relics with flat radio spectra, Kang et al. (2012), for instance, proposed a model in which pre-existing CR electrons with the momentum distribution corresponding to the observed radio spectral index are re-accelerated at weak shocks with . The sausage relic in CIZA J2242.8+5301 then requires pre-existing CR electrons with . It is interesting to note that this is close to the slope of the secondary electrons at production () due CR protons accelerated at energetic shocks in the cluster outskirts. The electrons with have the cooling time longer than the age of the universe. So we may conjecture a scenario in which the secondary electrons produced with are boosted to at shocks in the cluster outskirts, producing radio-relic-like structures. The acceleration or re-acceleration of CR electrons at shocks in clusters, compared to those of CR protons, involve additional complications such as injection, pre-existing CR population, and cooling.

## 6 Summary

The outskirts of galaxy clusters are dynamically active, reflecting disturbances due to mergers of substructures and continuous infall of gas along filaments. A manifestation of such activities is the formation of shock waves, which can be observed in radio and X-ray. In this paper, we have studied structure formation shocks in the cluster outskirts. Specifically, in a sample of 228 clusters from numerical simulations of the LSS formation with uniform grid resolution, we have identified the ICM shocks located in and studied their properties and the CR proton production via DSA there.

Main results are summarized as follows.

1. As previously known (e.g., Ryu et al., 2003), the ICM shocks existing in are mostly weak. But there are a number of shocks that are strong and energetic enough to produce substantial amounts of CR protons via DSA. Shocks with produce the largest amount of CR protons. Among shocks with , the fractions of shocks (actually shock zones) with , , and are , , and , respectively, in our sample. Shocks with the CR energy flux , were categorized as energetic shocks. The fraction of energetic shocks (again shock zones) is of the identified shocks. The shock with the largest in each cluster was designated as the most energetic shock (MES).

2. Infall shocks, which form by the infall of the WHIM from filaments, were separated from other shocks by employing the entropy, density, and Mach number criteria. In , of shocks with and about a half of shocks with are classified as infall shocks. Infall shocks are not as common as merger or turbulent shocks. But with relatively high Mach numbers () and large kinetic energy fluxes, they contribute to a significant fraction of CR production in clusters. We found that of energetic shocks and of the MESs are classified as infall shocks. And infall shocks produce of CRs in , when summed for all clusters in our sample.

3. Strong energetic shocks, including infall shocks, reside mostly in the cluster outskirts. Hence, CR protons are expected to be produced mostly in the outskirts and then advected into the core regions along with flow motions. The advection time scale is a substantial fraction of the age of the universe. Consequently, the radial profile of the CR pressure is expected to be broad, dropping off more slowly than that of the gas pressure, and might be even temporarily inverted peaking in the outskirts, as shown in previous simulations (e.g., Pfrommer et al., 2007; Vazza et al., 2012).

4. We have estimated the momentum distribution of the CR protons produced at shocks in . The volume-integrated momentum spectrum has the plower-law slope of . So the average Mach number of shocks in our sample clusters, weighted with CR production, is in the range of . It is greater than the characteristic Mach numbers of merger shocks (). This confirms that a substantial fraction of CRs are produced by energetic infall shocks in the cluster outskirts.

5. We suggest that some radio relics with flat radio spectrum such as the sausage relic in CIZA J2242.8+5301 () could be explained by primary electrons accelerated by energetic infall shocks with induced in the cluster outskirts.

The implications of our results on radio observations of clusters were briefly discussed, leaving detailed studies for future works.

We thank the anonymous referee for constructive comments. We also thank G. Brunetti and T. W. Jones for comments on the manuscript. HK thanks V. Petrosian and KIPAC for their hospitality during the sabbatical leave at Stanford University, where a part of work was done. SEH was supported by the National Research Foundation of Korea through grant 2007-0093860. DR was supported by research fund of Chungnam National University. HK was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2012-001065). RC was supported in part by grant NASA NNX11AI23G. Computing resources were in part provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

### Footnotes

1. affiliationmark:
2. affiliation: Corresponding author.
3. slugcomment: draft of March 21, 2018
4. The ratio of turbulent to gas pressure of %, for instance, corresponds to the turbulent Mach number of .
5. Here we distinguish infall shocks from external accretion shocks that decelerate never-shocked gas accreting onto clusters and filaments from void regions. Infall shocks are by nature also accretion shocks that stop previously shocked gas accreting onto clusters from filaments.
6. The relation between and is rather complicated and depends on cosmological parameters (e.g., Nakamura & Suto, 1997; Bryan & Norman, 1998; Eke et al., 1998). For the parameters we employed, approximately .
7. The volume-averaged CR energy production rate per unit volume at given radius is . Due to the large range of , the radial distribution of looks similar to that of .
8. The profile of , the slope of the momentum distribution of CR protons produced at shocks in the bin of , is similar to that of , since the CR production is dominated by shocks in the outer regions of clusters.

### References

1. Ackermann, M. et al. 2010, ApJ, 717, L71
2. Ackermann, M. et al. 2013, preprint, arXiv:1308.5654
3. Akamatsu, H., & Kawahara, H. 2013, PASJ, 65, 16
4. Akamatsu, H., Takizawa, M., Nakazawa, K., Fukazawa, Y., Ishisaki, Y., & Ohashi, T. 2012, PASJ, 64, 67
5. Arlen, T. et al. 2012, ApJ, 757, 123
6. Bell, A. R. 1978, MNRAS, 182, 147
7. Bell, A. R. 2004, MNRAS, 353, 550
8. Berezinsky, V. S., Blasi, P., & Ptuskim, V. S. 1997, ApJ, 487, 529
9. Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29
10. Bonafede, A., Giovannini, G., Feretti, L., Govoni, F., & Murgia, M. 2009, A&A, 494, 429
11. Brown, S., & Rudnick, L. 2011, MNRAS, 412, 2
12. Brüggen M., Bykov A., Ryu D., & Röttgering H., 2012, Space Sci. Rev., 166, 187
13. Brunetti, G., Blasi, P., Reimer, O., Rudnick, L., Bonafede, A., & Brown, S. 2012, MNRAS, 426, 956
14. Brunetti, G., Setti, G., Feretti, L., & Giovannini, G. 2001, MNRAS, 320, 365
15. Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80
16. Caprioli, D. 2012, J. Cosmology Astropart. Phys, 07, 038
17. Cen, R., & Chisari, N. E. 2011, ApJ, 731, 11
18. Clarke, T. E., & Enßlin, T. A. 2006, AJ, 131, 2900
19. Dermer, C. D. 1986, A&A, 157, 223
20. Dolag, K., & Enßlin, T. A. 2000, A&A, 362, 151
21. Drury, L. O’C. 1983, Rept. Prog. Phys., 46, 973
22. Eke, V. R., Navarro, J. F., & Frenk, C. S. 1998, ApJ, 503, 569
23. Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, A&A, 527, A99
24. Feretti, L., Giovannini, G., Govoni, F., & Murgia, M. 2012, A&A Rev., 20, 54
25. Gabici, S., & Blasi, P. 2003, ApJ, 583, 695
26. George, M. R., Fabian, A. C., Sanders, J. S., Young, A. J., & Russell, H. R. 2009, MNRAS, 395, 657
27. Hoeft, M., Brüggen, M., Yepes, G., Gottlober, S., & Schwope, A. 2008, MNRAS, 391, 1511
28. Kang, H. 2012, J. Korean Astron. Soc., 45, 127
29. Kang, H., Cen R., Ostriker, J. P., & Ryu, D. 1994, ApJ, 428, 1
30. Kang, H., & Jones, T. W. 2007, Astropart. Phys., 28, 232
31. Kang, H., & Ryu, D., 2010, ApJ, 721, 886
32. Kang, H., & Ryu, D., 2011, ApJ, 734, 18
33. Kang, H., & Ryu, D., 2013, ApJ, 764, 95
34. Kang, H., Ryu, D., Cen, R., & Ostriker, J. P. 2007, ApJ, 669, 729
35. Kang, H., Ryu, D., & Jones, T. W. 2012, ApJ, 756, 97
36. Komatsu, E., et al. 2011, ApJS, 192, 18
37. Lau E. T., Kravtsov A. V., & Nagai D., 2009, ApJ, 705, 1129
38. Li, S., Li, H., & Cen, R. 2008, ApJS, 174, 1
39. Lucek, S. G., & Bell, A. R. 2000, MNRAS, 314, 65
40. Markevitch, M., Forman, W. R., Sarazin, C. L., & Vikhlinin, A. 1998, ApJ, 503, 77
41. Markevitch, M., Gonzalez, A. H., David, L., Vikhlinin, A., Murray, S., Forman, W., Jones, C., & Tucker, W. 2002, ApJ, 567, 27
42. Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1
43. Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, ApJ, 559, 59
44. Miniati, F., Ryu, D., Kang, H., Jones, T. W., Cen, R., & Ostriker, J. 2000, ApJ, 542, 608
45. Nakamura, T. T., & Suto, Y. 1997, Prog. of Theor. Phys., 97, 49
46. Nuza, S. E., Hoeft, M., van Weeren, R. J., Gottlöber, S., & Yepes, G. 2012, MNRAS, 420, 2006
47. Ogrean, G. A., & Brüggen, M. 2013, MNRAS, 433, 1701
48. Ogrean, G. A., Brüggen, M., van Weeren, R. J., Röttgering, H., Croston, J. H., & Hoeft, M. 2013, MNRAS, 433, 812
49. Peebles P. J. E. 1980, The Large-Scale Structure of the Universe (Princeton: Princeton Univ. Press)
50. Pfrommer, C., & Enßlin, T. A., 2004, MNRAS, 352, 76
51. Pfrommer, C., Enßlin, T. A., Springel, V., Jubelgas, M., & Dolag, K. 2007, MNRAS, 378, 385
52. Pfrommer, C., & Jones, T. W. 2011, ApJ, 730, 22
53. Pfrommer, C., Springel, V., Enßlin, T. A., & Jubelgas, M. 2006, MNRAS, 367, 113
54. Pinzke, A., Oh, S. P., & Pfrommer, C. 2013 MNRAS, 435, 1061
55. Russell, H. R., Sanders, J. S., Fabian, A. C., Baum, S. A., Donahue, M., Edge, A. C., McNamara, B. R., & O’Dea, C. P. 2010, MNRAS, 406, 1721
56. Ryu, D., Kang, H., Cho, J., & Das, S. 2008, Science, 320, 909
57. Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, ApJ, 593, 599
58. Ryu, D., Ostriker, J. P., Kang, H., & Cen, R. 1993, ApJ, 414, 1
59. Ryu, D., Schleicher, D. R. G., Treumann, R. A., Tsagas, C. G., & Widrow, L. M. 2012, Space Sci. Rev., 166, 1
60. Schure, K. M., Bell, A. R, Drury, L. O’C., &. Bykov, A. M. 2012, Space Sci. Rev., 173, 491
61. Simionescu, A., Werner, N., Urban, O., Allen, S. W., Fabian, A. C., Mantz, A., Matsushita, K., Nulsen, P. E. J., Sanders, J. S., Sasaki, T., Sato, T., Takei, Y., & Walker, S. A. 2013, ApJ, 775, 4
62. Skillman, S. W., Hallman, E. J., O’Shea, B. W., Burns, J. O., Smith, B. D., & Turk, M. J. 2011, ApJ, 735, 96
63. Skillman, S. W., O’Shea, B. W., Hallman, E. J., Burns, J. O., & Norman, M. L. 2008, ApJ, 689, 1063
64. Skillman, S. W., Xu, H., Hallman, E. J., O’Shea, B. W., Burns, J. O., Li, H., Collins, D. C., & Norman, M. L. 2013, ApJ, 765, 21
65. van Weeren, R. J., Röttgering, H. J. A., Brüggen, M., & Cohen, A. 2009, A&A, 508, 75
66. van Weeren, R. J., Röttgering, H. J. A., Brüggen, M., & Hoeft, M. 2010, Science, 330, 347
67. van Weeren, R. J., Röttgering, H. J. A., Intema, H. T., Rudnick, L., Brüggen, M., Hoeft, M., & Oonk, J. B. R. 2012, A&A, 546, A124
68. Vazza, F., Brüggen, M., Gheller, C., & Brunetti, G. 2012 MNRAS, 421, 3375
69. Vazza, F., Brunetti, G, & Gheller, C. 2009a, MNRAS, 395, 1333
70. Vazza, F., Brunetti, G., Kritsuk, A., Wagner, R., Gheller, C., & Norman, M. 2009b, A&A, 504, 33
71. Vazza, F., Dolag, K., Ryu, D., Brunetti, G., Gheller, C., Kang, K., & Pfrommer, C. 2011, MNRAS, 418, 960
72. Vikhlinin, A., Kravtsov, A., Forman, W., Jones, C., Markevitch, M., Murray, S. S., & Van Speybroeck, L. 2006, ApJ, 640, 691
73. Voit, G. M., Kay, S. T., & Bryan, G. L. 2005, MNRAS, 364, 909
74. Walker, S. A., Fabian, A. C., Sanders, J. S., & George, M. R. 2012, MNRAS, 427, L45
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters