In a binary merger of two subclusters with comparable masses a pair of merger shocks are typically generated often manifesting as double radio relics. Using cosmological hydrodynamic simulations we identify major merger events with mass ratio mathcalM1mathcalM2lesssim4 and impact parameter brrm vir1lesssim1 where rrm vir1 is the virial radius of the larger subcluster. We analyze merger shock surfaces approximately 1 Gyr after the pericenter passage focusing on their morphology and the distribution of the Mach number Ms of their constituent shock zones. The shock surfaces exhibit an elongated shape with a minortomajor axis ratio of sim0.60.9 and cover the area of sim520 of the enclosed sphere. The area ratio of the two shock surfaces roughly scales with mathcalM1mathcalM2 typically positioning the larger shock ahead of the smaller subcluster. The axis connecting the two subclusters generally does not pass through the centers of the shock surfaces due to the nonzero impact parameter and the turbulent flows around them. The distribution of Ms of shock zones on each surface can be approximated by a lognormal function peaking at Msrmpeakapprox24.5 and extending up to sim10. The surfaceareaweighted and Xrayemissivityweighted average Mach numbers are comparable with langleMsranglermareaapprox2.34.4 and langleMsrangleXapprox24. In contrast the cosmicrayenergyfluxweighted average Mach numbers are higher with langleMsranglermCRapprox35. This discrepancy aligns with the differences between Mach numbers derived from Xray and radio observations of radio relic shocks. On the other hand we find that mostly langleMsrangleXgtrsim2 for simulated merger shocks although shocks with Mrm Xraylesssim2 are often reported in observations.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - High Energy Astrophysical Phenomena']
Morphology and Mach Number Distribution of Merger Shock Surfaces in Merging Galaxy Clusters
{'Morphology and Mach Number Distribution of Merger Shock Surfaces in Merging Galaxy Clusters': '<!-- image -->  \n1 Department of Physics, College of Natural Sciences, UNIST, Ulsan 44919, Korea  \n2 Department of Earth Sciences, Pusan National University, Busan 46241, Korea  \n(Received; Revised; Accepted)  \nSubmitted to The Astrophysical Journal', 'ABSTRACT': 'In a binary merger of two subclusters with comparable masses, a pair of merger shocks are typically generated, often manifesting as double radio relics. Using cosmological hydrodynamic simulations, we identify major merger events with mass ratio 1 < M 1 / M 2 ≲ 4 and impact parameter b/r vir , 1 ≲ 1, where r vir , 1 is the virial radius of the heavier subcluster. We analyze merger shock surfaces approximately 1 Gyr after the pericenter passage, focusing on their morphology and the distribution of the Mach number, M s , of their constituent shock zones. The shock surfaces exhibit an elongated shape with a minor-to-major axis ratio of ∼ 0 . 6 -0 . 9 and cover the area of ∼ 5 -20% of the enclosed sphere. The area ratio of the two shock surfaces roughly scales as A ss , 2 /A ss , 1 ∝ M 1 / M 2 , typically positioning the larger shock ahead of the lighter subcluster. The axis connecting the two subclusters generally does not pass through the centers of the shock surfaces, due to the nonzero impact parameter and ambient turbulent flows. The distribution of M s of shock zones on each surface can be approximated by a lognormal function, peaking at M s, peak ≈ 2 -4 . 5 and extending up to ∼ 10. The surface-area-weighted and X-ray-emissivity-weighted average Mach numbers are comparable, with ⟨ M s ⟩ area ≈ 2 . 3 -4 . 4 and ⟨ M s ⟩ X ≈ 2 -4. In contrast, the cosmic-ray-energy-flux-weighted average Mach numbers are higher with ⟨ M s ⟩ CR ≈ 3 -5. This discrepancy aligns with the differences between Mach numbers derived from X-ray and radio observations of radio relic shocks. On the other hand, we find that mostly ⟨ M s ⟩ X ≳ 2 for simulated merger shocks, although shocks with M X -ray ≲ 2 are often reported in observations.  \nKeywords: galaxies: clusters: general - methods: numerical - shock waves', '1. INTRODUCTION': "During the hierarchical formation of the large-scale structure (LSS) of the universe, galaxy clusters grow through the successive mergers of small subclusters and groups of galaxies. Such mergers are known to induce shock waves and turbulent flow motions in the intracluster medium (ICM). Energy dissipation via shocks and turbulence plays key roles not only in heating the ICM but also in amplifying magnetic fields and acceler-  \nCorresponding author: Dongsu Ryu  \ndsryu@unist.ac.kr  \nCorresponding author: Hyesung Kang  \nhskang@pusan.ac.kr  \nating cosmic rays (CRs) (e.g., Sarazin 2002; Ryu et al. 2003, 2008; Brunetti & Jones 2014; Porter et al. 2015; Roh et al. 2019; Wittor et al. 2021).  \nIn particular, binary mergers between two subclumps of comparable masses with small impact parameters are classified as major mergers and induce a pair of so-called 'merger shocks' after the pericenter passage of dark matter (DM) cores (e.g., Markevitch & Vikhlinin 2007; Ha et al. 2018). The overall picture of merger shocks has been extensively studied through numerical simulations. In simulations of idealized binary mergers with subclusters in hydrostatic equilibrium, it was shown that merger shocks exhibit different properties depending on the merger parameters, such as the mass ratio, impact parameter, viewing angle, etc (e.g., Gabici & Blasi 2003;  \nSpringel & Farrar 2007; van Weeren et al. 2011a; ZuHone 2011; Molnar & Broadhurst 2017).  \nThe properties of merger shocks were also described in the context of the LSS formation (e.g., Paul et al. 2011; Schmidt et al. 2017; Ha et al. 2018; glyph[suppress]Lokas 2023; Lee et al. 2024; Finner et al. 2024). For instance, analyzing merging clusters found in cosmological hydrodynamic simulations, Ha et al. (2018) described quantities, such as the time variations of the shock speed, v s , the shock Mach number, M s , and the energy dissipation at shock interfaces, in the inherently turbulent ICM 1 . While a variety of mergers with diverse merger parameters are present in these simulations, Ha et al. (2018) concentrated on the cases of nearly head-on collision with the mass ratio of M 1 / M 2 ≈ 2, replicating the characteristics of the well-studied Sausage relic in the cluster CIZA J2242.8+5301 (e.g., van Weeren et al. 2010; Hoang et al. 2017). As illustrated in Figure 1 of Ha et al. (2018), as two subclusters approach and gas clumps are compressed, 'equatorial shocks' first appear in the equatorial plane, propagating toward the direction perpendicular to the merger axis. After the DM core passage, two 'axial shocks' launch along the merger axis in front of subclusters (see also Figure 1 of van Weeren et al. 2011a).  \nMerger shocks often manifest as double radio relics in the outskirts of merging clusters (e.g., van Weeren et al. 2011a; Hoang et al. 2018; Golovich et al. 2019). With observations of around a hundred cases to date, it is widely accepted that radio relics exhibit diffuse synchrotron radiation emitted by CR electrons accelerated at merger shocks (see Feretti et al. 2012; Bruggen et al. 2012; Brunetti & Jones 2014; van Weeren et al. 2019, for reviews). Although it has been customary to associate a merger shock with a single Mach number, considering the curved, or often irregular, morphology of radio relics and the turbulent nature of the ICM, the correct picture should be that a merger shock surface is composed of many 'shock zones' characterized by different M s (e.g., Ha et al. 2018; Roh et al. 2019; Botteon et al. 2020; Wittor et al. 2021). Thus, the Mach number estimated from radio or X-ray observations should represent mean values, ⟨ M s ⟩ radio and ⟨ M s ⟩ X -ray , for an ensemble of shock zones associated with each merger shock surface (e.g., Hong et al. 2015; Rajpurohit et al. 2020; DominguezFernandez et al. 2021).  \nIn radio observations, radio relics frequently appear as arc-shaped patches at a distance of d s ∼ 1 -2 Mpc  \nfrom the cluster center, usually around the virial radius, r vir . The Mach numbers of radio relic shocks are estimated using either the 'injection spectral index', α inj (or α sh ), of the radio spectrum immediately behind the shock in the outer edge of observed relics, i.e., M radio ≈ [(3 + 2 α inj ) / (2 α inj -1)] 1 / 2 , or the 'integrated spectral index', α int , of the volume-integrated spectrum, i.e., M radio ≈ [( α int + 1) / ( α int -1)] 1 / 2 (e.g., Kang 2011). Radio observations indicate typically M radio ≈ 2 -5 (van Weeren et al. 2019). On the other hand, in X-ray observations, merger shocks are detected as discontinuities in the surface brightness or temperature profiles that can be used to estimate the Mach number, M X -ray (e.g., Markevitch et al. 2002; Ogrean et al. 2013; Dasadia et al. 2016; Sanders et al. 2022). These shocks are typically found at d s ≳ 1 Mpc and mostly weak with M X -ray ∼ 1 . 2 -4 (e.g., Eckert et al. 2016; Urdampilleta et al. 2018).  \nFor many observed radio relics, the Mach numbers derived from radio and X-ray observations do not coincide, exhibiting the so-called Mach number discrepancy (e.g., Wittor et al. 2021). In most cases, Mach numbers inferred from X-ray observations are lower than those inferred from radio observations, i.e., M X -ray ≲ M radio (e.g., Akamatsu & Kawahara 2013; van Weeren et al. 2016). For instance, in the so-called Toothbrush relic in 1RXS J0603.3, M radio ≈ 3 . 3 with α sh ≈ 0 . 7 and M radio ≈ 3 . 78 with α int ≈ 1 . 15 (Rajpurohit et al. 2018), whereas M X -ray ≲ 2 (Itahana et al. 2015; van Weeren et al. 2016). However, in some relics, the opposite is found; e.g., for a radio relic in Abell 521, M radio ≈ 2 . 27 and M X -ray ≈ 3 . 4 are reported (Giacintucci et al. 2006).  \nAs noted above, a merger shock is a surface comprising many shock zones with varying Mach numbers. Previous studies have suggested that radio observations tend to pick up regions with high CR acceleration and, consequently, higher M s , while X-ray observations preferentially reveal regions with higher densities and thus higher X-ray emissivity, which are typically associated with lower M s (e.g., Hong et al. 2015; Ha et al. 2018; Roh et al. 2019; Botteon et al. 2020). This proposition could explain why M radio ≳ M X -ray . However, complex flows surrounding merger shocks in the turbulent ICM could lead to the opposite case. These highlight the needs to understand the detailed properties of merger shocks and the dynamics of surrounding flows and to clarify the nature of observed radio relics.  \nIn this work, we study the morphological characteristics of merger shock surfaces and the Mach number distribution of shock zones in major binary merging clusters with different mass ratios and impact parameters using cosmological hydrodynamic simulation data. To  \nachieve this, we perform a set of structure formation simulations, identify synthetic merging clusters formed through major binary mergers, and track the time evolution of merger events such as the pericenter passage and the development of merger shocks. Using specific criteria, we try to isolate the surfaces of merger shocks from shocks generated by turbulence in the background ICM and the infall of the warm-hot intergalactic medium (WHIM) along filaments. We then analyze the shape, size, and orientation of the shock surfaces, and examine the probability distribution function (PDF) of M s for shock zones on these surfaces. Additionally, we estimate the average Mach numbers weighted by shock surface area, X-ray emissivity, and CR energy flux, and discuss the implications for observations of radio relics. Figure 1.  \nThe paper is organized as follows: In Section 2, we outline the numerical setup, including the generation of synthetic merging clusters and the method for identifying merger shock surfaces. Section 3 presents the results, discussing the morphological characteristics of merger shock surfaces and the Mach number distribution of shock zones. In Section 4, we compare the average Mach numbers of merger shocks with those observed in radio relic shocks. Finally, a brief summary is provided in Section 5.", '2.1. Simulation Setup': 'We generate samples of merging clusters and merger shocks from simulation data for the LSS formation of the universe. Simulations are performed using the particlemesh/Eulerian cosmological hydrodynamic code described in Ryu et al. (1993). We adopt the standard ΛCDM cosmology model with the following parameters: baryon density Ω b = 0.046, dark matter density Ω DM = 0.234, cosmological constant Ω Λ = 0.72, Hubble parameter h ≡ H 0 /(100 km s -1 Mpc -1 ) = 0.7, primordial spectral index n = 0.96, and rms density fluctuation σ 8 = 0.82, complying with WMAP 9 data (Hinshaw et al. 2013). A cubic box of the comoving size of 100 h -1 Mpc with periodic boundaries is employed and divided into 2048 3 uniform grid zones. This results in a uniform spatial grid resolution of ∆ l = 48.8 h -1 kpc.  \nIn our simulations, non-gravitational processes, such as radiative cooling and galaxy formation feedback, are not included, as their effects on the dynamical evolution of the ICM on megaparsec scales are expected to be insignificant. Magnetic fields are also not included; given the expected magnetic field strength of ≲ 1 µ G at cluster outskirts (e.g., Ryu et al. 2008; Vazza et al. 2014), the magnetic energy would be an order of magnitude smaller than the kinetic energy, and hence magnetic fields are  \n𝑑𝑑 ≈ 𝑟𝑟  \n(a) pre-merger approach  \nFigure 1 illustrates the schematic picture of a binary major merger, showing our definitions of the impact pa-  \n<!-- image -->  \n(b) after DM pericenter passage  \nFigure 1. Schematic picture of a binary merger with a nonzero impact parameter b during (a) the pre-merger approach and (b) far after the pericenter passage. Here, r 200 is the radius within which the mean total density of the cluster attains 200 times the critical density, ρ c , of the universe. The subscripts 1 and 2 denote heavier and lighter subclusters, respectively. A pair of merger shocks develop after the pericenter passage, and the one in front of subcluster 1 is labeled as relic shock 1 (blue), while the other in front of subcluster 2 is labeled as relic shock 2 (magenta).  \nunlikely to significantly affect the properties of merger shocks, such as their morphology and the Mach number distribution at the shock surface. On the other hand, on microscopic scales, the behavior of magnetic fields, such as amplification and dissipation through kinetic instabilities and wave-particle interactions, should play a crucial role in particle acceleration and subsequent synchrotron emission behind the shock (e.g., Kang et al. 2019; Ha et al. 2021, 2022), but such processes are beyond the scope of this study.', '2.2. Merging Cluster Sample': "From seven different sets of simulations, we identify clusters with an X-ray-weighted temperature T X ≳ 2 keV at redshift z = 0. Here, T X is calculated by weighting the gas temperature, T , with the bremsstrahlung emission, ε ff ∝ T 1 / 2 ρ 2 gas , across the cluster volume (see Kang et al. (1994) and Hong et al. (2015) for details), where ρ gas is the gas density. Among them, we find more than 20 clusters that have gone through major binary mergers. We construct a sample of 12 merging clusters satisfying the following three criteria: (1) mass ratio 1 < M 1 / M 2 ≲ 4; (2) impact parameter b ≲ r vir , 1 ; (3) T X ≳ 2 KeV at the redshift of z relic ≲ 0 . 5. Hereafter, we denote heavier and lighter subclusters with the subscripts 1 and 2, respectively; r vir , 1 is the virial radius of the heavier subcluster, and z relic is the redshift at ∼ 1 Gyr after the pericenter passage, the optimal time for merger shock surfaces to be observed as radio relics (see Ha et al. 2018).  \n+  \n𝑟𝑟  \nTable 1. Merging Galaxy Cluster Sample  \nrameter and two relic shocks. In panel (a), which depicts a pre-merger phase, two subclusters are about to touch with a distance of d ≈ r 200 , 1 + r 200 , 2 . Here, r 200 is the radius within which the mean total (gas plus DM) density of the cluster attains 200 times the critical density, ρ c , of the universe. The impact parameter, b , is defined as the distance perpendicular to the paths of the subcluster centers. Typically, this corresponds to ∼ 1 . 5 -2 Gyr before the pericenter passage. Panel (b) depicts a stage following the pericenter passage, showcasing two merger shocks, relic shock 1 in front of subcluster 1 and relic shock 2 in front of subcluster 2. We use the term 'relic shock' because, as noted in the introduction, merger shocks are often observed as radio relics. Given that most binary radio relics in major mergers are observed in the post-pericenter passage phase (Golovich et al. 2019), our analysis focuses on the stage shown in panel (b).  \nTable 1 lists the 12 merging clusters in our sample. Here, r vir ≈ 1 . 36 r 200 (e.g., Roncarelli et al. 2006; Reiprich et al. 2013), and M 1 and M 2 are the total (gas plus DM) mass contained within the subcluster's virial radius at the epoch of d ≈ r 200 , 1 + r 200 , 2 . The masses and temperatures of the merging clusters in this  \n<!-- image -->  \nFigure 2. 3D view of shock zones in the comoving volume of (4 . 88 Mpc) 3 that encompasses Cluster 3 (see Table 1). Color displays the Mach number of shock zones. (a) All shock zones within the volume. (b) Shock zones satisfying the density/temperature and distance criteria (see the main text). Here, only the shock zones within ⟨ d s ⟩ ± 1 σ d s are included, where ⟨ d s ⟩ is the mean distance of shock zones from the X-ray center of the cluster, and σ d s is its standard deviation. The images are captured at z relic .  \n<!-- image -->  \nstudy are somewhat lower than typical values for observed merging clusters (see Section 4) due to the limited computational volume in our simulations. Therefore, we present results primarily for normalized dimensionless quantities, assuming these are approximately scaleindependent.", '2.3. Shock Zones and Merger Shock Surfaces': 'In merged clusters, we identify shocks using the algorithm introduced in Ryu et al. (2003). Along the axis directions, grid zones are tagged as shock zones if they satisfy the following three conditions: (1) ∇ · v < 0, i.e., locally converging flow; (2) ∆ T × ∆ ρ > 0, i.e., the same sign of gas temperature and density gradients; and (3) | ∆log T | > 0 . 11, i.e., the temperature jump larger than that of Mach number M s = 1 . 3. Then, M s is calculated by solving the Rankine-Hugoniot temperature jump, T 2 /T 1 = (5 M 2 s -1)( M 2 s + 3) / (16 M 2 s ). Here, the subscripts 1 and 2 represent the preshock and postshock quantities, respectively. After applying these procedures for all three axis directions, the Mach number of a shock zone is defined as the maximum value of the three Mach numbers along the x , y and z axes, M s = max( M s,x , M s,y , M s,z ). Many grid zones are identified to contain weak shocks, but weak shocks are not energetically important. Thus, we only consider the shock zones with M s ≥ 1 . 5.  \nIn the ICM, shocks can be induced by dynamical activities other than mergers of subclusters, such as turbulent flow motions and the infall of the WHIM along filaments. To isolate shock zones that belong to the merger-driven shock surfaces, we employ the following steps. In the first step, considering merger shocks typically appear  \nFigure 3. Merging process in Cluster 10 with M 1 / M 2 ≃ 2 . 96 and b/r vir , 1 ≃ 0 . 3 (see Table 1). From top to bottom, 2D slice images of the gas density, ρ gas / ⟨ ρ gas ⟩ , the DM density, ρ DM / ⟨ ρ DM ⟩ , and the gas temperature, T , in the comoving area of (8 h -1 Mpc) 2 around the cluster are shown. Here, ⟨ ρ gas ⟩ and ⟨ ρ DM ⟩ are the mean gas and DM densities of the universe, and T is in units of Kelvin. The pericenter passage occurs at z ≈ 0 . 1, while a pair of merger shocks would appear as radio relics at z ≲ 0 . 05. Relic shock 2 in front of subcluster 2 is larger than relic shock 1 in front of subcluster 1. The black dot at z = 0 draws the circle of the virial radius, r vir .  \n<!-- image -->  \naround the virial radius (see the next section), we select shock zones satisfying (1) -1 . 4 ≤ log 10 ρ 1 / ⟨ ρ ⟩ ≤ 0 . 3 and (2) -1 . 0 ≤ log 10 T 1 / ⟨ T ⟩ ≤ 0 . 1. Here, ρ 1 and T 1 represent the preshock density and temperature, respectively, and ⟨ ρ ⟩ and ⟨ T ⟩ are the mean values calculated within the virial radius of each cluster. The upper bounds for ρ 1 and T 1 help eliminate turbulence-induced shocks in the dense and hot ICM inside the virial radius, while the lower bounds are imposed to exclude shocks induced by the infall of the WHIM along filaments in cluster outskirts. These criteria have been empirically determined by examining the merger shocks at z relic . In the second step, using the selected shock zones, we estimate the mean distance ⟨ d s ⟩ of each relic shock surface from the cluster center and its standard deviation, σ d s . Throughout the paper, the cluster center is defined as the peak of X-ray emission within the merged cluster. Only those within ⟨ d s ⟩ ± 1 σ d s are chosen as parts of the relic shock surface. We then iteratively update ⟨ d s ⟩ and σ d s and refine the list of constituent shock zones to achieve the best result.  \nFigure 2 illustrates three-dimensional (3D) views of the shock surfaces in Cluster 3 at z relic ≈ 0 . 1 (see Table 1). The left panel displays all shock zones within the volume, whereas the right panel shows the shock zones after the refinement according to the above criteria. Below we examine the properties of refined merger shock surfaces.', '3.1. Morphological Properties of Merger Shock Surfaces': "Mergers of subclusters take place during the hierarchical formation of the LSS of the universe; the resulting merger shocks, which develop in the inherently turbulent ICM, usually have complex morphology. Figure 3 displays two-dimensional (2D) slices for the gas density (top panels), DM density (middle panels), and gas temperature (bottom panels) in Cluster 10 (see Table 1). At z = 0 . 18, two prominent clumps approach each other with M 1 / M 2 ≃ 2 . 96 and b/r vir , 1 ≃ 0 . 3. As these clumps close in with a relative velocity of about 10 3 km s -1 , shocks develop along the merger axis and  \n(a)  \n(b)  \nFigure 4(b) plots the mean distance of merger shocks from the cluster center, ⟨ d s ⟩ , as a function of the virial radius of the merged cluster, r vir , at z relic . For our sample clusters with r vir ≈ 1 . 0 -2 . 3 Mpc (proper distance), ⟨ d s ⟩ /r vir ranges ∼ 0 . 6 -1 . 2 with an average value of ∼ 0 . 9. This indicates that at z relic , ∼ 1 Gyr after the pericenter passage, although there is a significant scatter, merger shocks are typically found around the virial radius, consistent with observations of radio relic locations (see the introduction). It is worth noting that Zhang et al. (2019) proposed a 'habitable zone' at distances ≳ r 500 , or ≳ 0 . 5 r vir , beyond which moderately strong shocks appear. Here, r 500 is the radius within which the mean total density is 500 times ρ c . Although their simulation setup and merger parameters differ from ours, the merger shocks in Figure 4(b) are indeed in the habitable zone, aligning with their suggestion. Additionally, the figure shows that there is almost no or only weak dependence of ⟨ d s ⟩ /r vir on r vir , implying that on average, ⟨ d s ⟩ tends to be larger in more massive clusters with larger r vir .  \n<!-- image -->  \n𝑣𝑣𝑣𝑣𝑣𝑣  \nFigure 4. (a) 3D view of merger shocks at z relic for four representative clusters. Color displays the Mach number of shock zones. The comoving volume of the boxes ranges approximately (4 -5 Mpc) 3 ; the box size varies to optimize the visualization of shock surfaces. Merger shocks for all 12 sample clusters in Table 1 are provided in Figure A1. (b) Mean distance, ⟨ d s ⟩ , between the merger shock surface and the X-ray center of the merged cluster, normalized to the virial radius of each cluster, r vir , versus r vir for 24 merger shocks in Figure A1. Here, r vir is given as the proper distance. The vertical error bars represent the ± 1 σ d s deviation. The horizontal dashed line indicates the average value of ⟨ d s ⟩ /r vir ≈ 0 . 9. Blue dots and lines represent relic shock 1, while red dots and lines represent relic shock 2, the same colors as in Figure 1. The same color scheme is used in subsequent figures.  \n= 2.2, = 0.46  \n= 3.08,  \n= 0.64  \nbecome visible at z = 0 . 14. The epoch near the pericenter passage occurs around z = 0 . 1. Following this, two merger shocks are observed propagating outward over z = 0 . 1 -0 . 0. The shock surfaces and the turbulent nature of the ICM are most clearly depicted in the temperature images. On the other hand, the gas density images reveal multiple minor mergers and secondary infalls that contribute to turbulent flow motions in the ICM. Although not shown here, the shock speed, v s , and also M s initially increase with time and then tend to converge, though significant fluctuations are observed (see Figure 5 of Ha et al. 2018).  \nFigure 4(a) presents 3D views of merger shocks at z relic for four of our sample clusters (Clusters 1, 3, 6, and 11), highlighting pairs of shock surfaces moving in opposite directions. Merger shocks for all 12 sample clusters are shown in Figure A1. In Clusters 6 and 11 with relatively large M 1 / M 2 , a larger shock is in front of the lighter dark matter clump (relic shock 2), while a smaller shock appears in front of the heavier dark matter clump (relic shock 1). Additionally, the figure illustrates the asymmetry, caused by nonzero impact parameters and the turbulent ICM.  \n𝑟𝑟  \n(Mpc)  \nThe left panel of Figure 5 schematically depicts the geometry of merger shock surfaces as elliptical sections  \n0.4  \n/ 𝑟𝑟 𝑣𝑣𝑣𝑣𝑣𝑣 , 1 𝑏𝑏 Figure 5. Left Panel: Schematic picture showing the geometry of relic shocks in a binary merger far after the pericenter passage. The shock surfaces are approximated as ellipses with major axis a ss and minor axis b ss . A ss denotes the surface area of shocks and is calculated as A ss ≈ 1 . 19(∆ l ) 2 N sh , where N sh is the number of shock zones associated with a particular shock surface. Right Panel: (a) Axial ratio of merger shock surfaces, b ss /a ss , for the 24 merger shocks in Figure A1. (b) Ratio of the area of merger shock surfaces, A ss , to the area of the shock spheres of radius ⟨ d s ⟩ , A sphere = 4 π ⟨ d s ⟩ 2 . (c)-(d) Ratio of A ss , 2 /A ss , 1 versus the mass ratio of two subclusters, M 1 / M 2 , and the normalized impact parameter, b/r vir , 1 , for 12 sample clusters. The horizontal dashed lines in (a)-(b) denote the average values for the 24 merger shocks. The solid line in panel (c) draws y = x .  \n<!-- image -->  \nof spherical shells at distance ⟨ d s ⟩ , described by their semi-major axis a ss and semi-minor axis b ss . The major axis ranges 2 a ss ≈ 1 . 2 -3 . 6 Mpc, while the minor axis ranges 2 b ss ≈ 0 . 7 -2 . 8 Mpc, in terms of the proper distance. As shown in Figure 5(a), the aspect ratio is generally b ss /a ss ≳ 0 . 6 with an average of ∼ 0 . 77, except for an outlier with b ss /a ss ≈ 0 . 4 observed for relic shock 1 in Cluster 5. So, while most merger shock surfaces normally appear moderately elongated, highly elongated, stripe-like surfaces can also form depending on the merger history of clusters as well as on the flow dynamics in the turbulent ICM.  \nWe define the 'shock sphere' as a sphere with radius ⟨ d s ⟩ and surface area A sphere = 4 π ⟨ d s ⟩ 2 . For merger shock surfaces, the area is calculated as A ss ≈ 1 . 19(∆ l ) 2 N sh , where N sh is the number of shock zones associated with a given merger shock. The factor of 1.19 accounts for the mean projected area within a 3D zone with random shock normal orientations. In Figure 5(b), the normalized area of merger shock surfaces, A ss /A sphere , is plotted. For most relic shocks, A ss /A sphere ranges ∼ 0 . 05 -0 . 2, with the exception of relic shock 2 in Cluster 8, which has A ss /A sphere ≈ 0 . 36. The average value is A ss /A sphere ∼ 0 . 15, indicating that  \nmerger shock surfaces typically occupy a small fraction of the surface area of shock spheres.  \nThe morphological properties of merger shock surfaces can vary depending on merger parameters. In idealized binary mergers with subclusters in hydrostatic equilibrium, it was shown that the mass ratio influences the relative size of two merger shock surfaces, while the impact parameter affects the degree of asymmetry around the merger axis (e.g., van Weeren et al. 2011a; Molnar & Broadhurst 2017). Figure 5(c) and (d) show the ratio of the surface areas of two merger shocks, A ss , 2 /A ss , 1 , as a function of merger parameters, M 1 / M 2 and b/r vir , 1 , for our sample clusters. The ratio A ss , 2 /A ss , 1 tends to increase with increasing M 1 / M 2 , consistent with the findings of previous studies on idealized binary mergers mentioned above. In contrast, while there is considerable scatter, no clear dependence of A ss , 2 /A ss , 1 on b/r vir , 1 is found.", '3.2. Mach Numbers as a Function of Angular Distance from the Merger Axis': "As an effort to describe the morphological characteristics of relic shocks in the turbulent ICM, we examine the distribution of M s in merger shock surfaces as a function of the position angle, θ , which is defined as the angle  \nFigure 6. (a) 3D view of merger shocks, superimposed on 3D volume-rendered images of the DM density, ρ DM / ⟨ ρ DM ⟩ , for the four clusters shown in Figure 4. The color of merger shocks displays the 'position angle' of shock surfaces, θ , which is the angle between the axis passing through the two DM density peaks (indicated by black line segments) and the line connecting the midpoint of these peaks to shock zones. The images are at z relic as in Figure 4, but with different viewing angles. The corresponding images for all 12 sample clusters are provided in Figure A2. (b) Average Mach number of shock zones associated with relic shocks 1 (blue) and 2 (red) as a function of θ with bins of ∆ θ = 5 · for the four sample clusters. Error bars represent ± 1 σ deviations of M s . Similar plots for all 12 sample clusters are shown in Figure A3.  \n<!-- image -->  \nlog  \n𝐷𝐷𝐷𝐷  \n𝜃𝜃  \n(°)  \n𝜌𝜌  \nbetween the axis connecting two DM density peaks and the line extending from the midpoint of these peaks to the given shock zone. Figure 6(a) depicts the 3D distribution of θ for merger shocks, on top of the DM density distribution, for the four clusters shown in Figure 4(a). The corresponding images for all 12 sample clusters are available in Figure A2. Figure 6(b) plot the average and standard deviation of M s in bins of ∆ θ as a function of θ for shock zones that belong to relic shocks 1 (blue) and 2 (red) in the four sample clusters. Similar plots for all 12 sample clusters are provided in Figure A3.  \nIn the idealized scenario of a spherical bow shock propagating through a uniform medium, the Mach number distribution is expected to vary as ⟨ M s ⟩ ∝ cos θ (see, e.g., Petrinec & Russell 1997, for further discussion). Russell et al. (2022) investigated the angular dependence of the Mach number for a merger shock in Abell 2146 using X-ray data, finding that it roughly follows the cosine behavior. In contrast, our study examines the angular dependence of the Mach number in 3D. As shown in Figures 6(b) and A3, while some merger shocks, such as relic shock 2 of Cluster 3, exhibit behaviors close to the idealized pattern, many do not follow this trend. Specifically, some shocks exhibit peaks at large nonzero θ values, rather than around θ ∼ 0 · , with a roughly symmetrical decline on either side of the peak (e.g., relic shock 1 of Cluster 11). Others display significant fluctuations in the Mach number distribution with no distinct peak. As a matter of fact, in most cases, the axis connecting DM peaks does not go through the center of  \nshock surfaces. This is partly due to nonzero impact parameters in our sample clusters and also due to the fact that merger shock surfaces are disturbed by turbulent flow motions and infalling clumps in the ICM as they propagate from the core to outskirts. Consequently, the center of the shock surface fluctuates and becomes misaligned with the axis connecting the DM peaks. For example, for relic shock 1 of Cluster 1, the shock surface is significantly tilted, so there is no shock zone for θ ≲ 5 · . We point out that Cluster 1 has experienced a major merger at z = 0 . 22, followed by a minor merger with a mass ratio of ∼ 6 in the propagation direction of relic shock 1 at z = 0 . 18. And the shock properties are analyzed at z relic = 0 . 14. In this relic shock, the shock surface morphology is affected by the late minor merger (see the shock on the right hand side of Cluster 1 in Figure 4(a)) and offset from the axis connecting DM peaks.", '3.3. Mach Number Distribution on Merger Shock Surfaces': "The PDFs for M s of shock zones across merger shock surfaces are shown in Figure 7(a) for the four representative clusters. The corresponding PDFs for all 12 sample clusters are provided in Figure A4. In general, the PDFs exhibit positively skewed distributions. We fit the PDFs to a log-normal function, with the fittings shown as gray lines in the figures. The log-normal fits peak around M s, peak ≈ 2 -4 . 5, with a high Mach number tail extending up to M s ∼ 10. This tail is primarily due to  \nFigure 7. (a) PDF of the Mach number of shock zones associated with relic shocks 1 (blue) and 2 (red) in the four clusters shown in Figure 4. The vertical dashed lines denote the shock-surface-area-weighted average Mach number, ⟨ M s ⟩ area , while the gray solid lines show the log-normal fitting. The corresponding PDFs for all 12 sample clusters are provided in Figure A4. (b) Standard deviation normalized to the mean value, σ M s / ⟨ M s ⟩ area , as a function of ⟨ M s ⟩ area for the PDFs of 24 merger shocks in 12 sample clusters. The subscript 'area' is omitted for simplicity. The horizontal dashed line represents the average value of σ M s / ⟨ M s ⟩ ≈ 0 . 37. (c) σ ln M s /µ ln M s versus µ ln M s , where µ ln M s and σ ln M s are the mean and standard deviation for the fitted log-normal distribution of the Mach number.  \n<!-- image -->  \nintermittent accretions of the WHIM, as merger shocks are often attached to infalling filaments. These Mach number distributions agree with the findings of previous numerical studies (e.g., Botteon et al. 2020; Wittor et al. 2021).  \nFor the further analysis of their statistical characteristics, we calculate (1) the mean and standard deviation of M s , ⟨ M s ⟩ area and σ M s , using the actual PDFs (blue and red lines in Figures 7(a) and A4, and (2) the mean and standard deviation, µ ln M s and σ ln M s , using the fitted log-normal distributions (gray lines). The subscript 'area' is added to ⟨ M s ⟩ area to specify that it is actually the shock-surface-area-weighted average Mach number. Figures 7(b) and (c) demonstrate that σ M s / ⟨ M s ⟩ area and σ ln M s /µ ln M s are largely independent of ⟨ M s ⟩ area and µ ln M s , with average values of σ M s / ⟨ M s ⟩ area ≈ 0 . 37 and σ ln M s /µ ln M s ≈ 0 . 33, respectively. This implies that the PDFs of M s are more or less universal, although there is scatter owing to diverse merger histories of clusters, including the infall of small clumps, as well as turbulent flow motions in the ICM.", '3.4. Average Mach Numbers of Merger Shock Surfaces': "In addition to M s, peak and ⟨ M s ⟩ area , the Mach number distribution of shock zones on merge shock surfaces can be characterized by other representative values, such as the average Mach numbers weighted with observable quantities. Specifically, we examine the Xray-emissivity-weighted average Mach number, ⟨ M s ⟩ X , and the CR-energy-flux-weighted average Mach number, ⟨ M s ⟩ CR . Here, ⟨ M s ⟩ X is obtained by weighing  \nM s with the bremsstrahlung emission, ε ff , over shock zones on the given shock surface. For ⟨ M s ⟩ CR , M s is weighted by the CR energy flux, f CR = η e ( M s ) · f ϕ , where f ϕ = (1 / 2) ρ 1 v s 3 represents the shock kinetic energy flux. The CR electron acceleration efficiency is approximated as η e ≈ 0 . 01 η , where the CR proton acceleration efficiency ranges η ( M s ) ≈ 3 . 6 × 10 -3 -0 . 01 for 2 . 25 ≲ M s ≲ 5 . 0 (Ryu et al. 2019). In this study, we regard ⟨ M s ⟩ CR as a proxy for radio-emissivity-weighted average Mach number, ⟨ M s ⟩ radio .  \nIn general, M s, peak tends to be slightly smaller than ⟨ M s ⟩ area ≈ 2 . 3 -4 . 4, as expected with the log-normal shape of the Mach number PDFs, as shown in Figure 8(a). The figure also illustrates that ⟨ M s ⟩ X ∼ 2 -4 is comparable to ⟨ M s ⟩ area . On the other hand, ⟨ M s ⟩ CR ∼ 3 -5 is greater than both ⟨ M s ⟩ area and ⟨ M s ⟩ X . This trend reflects our CR acceleration model, η ( M s ), which predicts higher efficiencies for higher M s , as noted in the introduction. These results are consistent with those of previous studies. For example, using simulations with a cosmological magnetohydrodynamic code, ENZO, Wittor et al. (2021) reported average Mach numbers of ⟨ M s ⟩ X -ray ≈ 1 -4 and ⟨ M s ⟩ radio ≈ 3 -5 for merger shocks; in their work, ⟨ M s ⟩ X -ray and ⟨ M s ⟩ radio were estimated mimicking observations, including projection effects and synchrotron modeling. In summary, M s, peak ≲ ⟨ M s ⟩ area ≈ ⟨ M s ⟩ X ≲ ⟨ M s ⟩ CR , indicating that while ⟨ M s ⟩ X provides a reasonable estimate of ⟨ M s ⟩ area , ⟨ M s ⟩ CR is typically higher. However, there are exceptions; in relic shock 1 of Cluster 1, where a secondary minor merger has disturbed the shock surface, ⟨ M s ⟩ CR is  \nFigure 8. (a) Peak of the Mach number PDF, M s, peak (magenta), X-ray-emissivity-weighted average Mach number, ⟨ M s ⟩ X (red), and CR-energy-flux-weighted average Mach number, ⟨ M s ⟩ CR (blue), versus shock-surface-area-weighted average Mach number, ⟨ M s ⟩ area , for 24 merger shocks. (b) ⟨ M s ⟩ X (red) and ⟨ M s ⟩ CR (blue) as a function of the normalized distance of merger shocks from the cluster center, ⟨ d s ⟩ /r vir .  \n<!-- image -->  \nclose to ⟨ M s ⟩ area , while ⟨ M s ⟩ X is smaller (see the points with ⟨ M s ⟩ area = 4 . 4 in the figure).  \nFigure 8(b) presents ⟨ M s ⟩ X and ⟨ M s ⟩ CR , both of which are linked to observable quantities, as a function of the normalized mean distance of merger shocks from the cluster center, ⟨ d s ⟩ /r vir . They generally tend to increase with ⟨ d s ⟩ /r vir . This trend is partly due to the fact that, on average, the gas density decreases and the shock speed increases outwardly in merged clusters (Ha et al. 2018). There is substantial scatter, which would be again attributed to the merger histories of clusters and the influence of turbulent flow motions in the ICM.  \nFigure 9. ⟨ M s ⟩ X (top panels) and ⟨ M s ⟩ CR (bottom panels) versus mass ratio M 1 / M 2 (left panels) and impact parameter b/r vir , 1 (right panels) for 24 merger shocks in 12 sample clusters.  \n<!-- image -->  \nThe Mach number distribution of merger shocks, and hence the representative Mach numbers, could depend on merger parameters, similar to the morphological properties discussed in Section 3.1. Figure 9 shows ⟨ M s ⟩ X and ⟨ M s ⟩ CR as functions of M 1 / M 2 and b/r vir , 1 . The average Mach numbers tend to decrease with increasing M 1 / M 2, although there is relatively large scatter, but the dependence on b/r vir , 1 appears to be only marginal. In Ha et al. (2018), where cases with b ≈ 0 and M 1 / M 2 ≈ 2 were examined, relic shock 1 ahead of the more massive subcluster tends to be stronger with higher Mach numbers than relic shock 2, albeit with large fluctuations (see their Figure 5(a)). However, this tendency is not obvious in our merger shock samples from clusters with nonzero impact factors. As can be seen in Figures 9(a) and (b), in only 5 of 12 clusters, ⟨ M s ⟩ X for relic shock 1 is higher than for relic shock 2. Similarly, in Figures 9(c) and (d), in only 6 of 12 clusters, ⟨ M s ⟩ CR for relic shock 1 exceeds that for relic shock 2. Given the complex interplay of various factors such as merger parameters, merger history, and ICM turbulence, however, we argue that a larger cluster sample would be necessary to establish a more robust conclusion on the relative strengths of M s between relic shock 1 and 2.  \nTable 2. Properties of Observed Merger shocks in Galaxy Clusters Hosting Double Radio Relics  \nNote - Column1: Name of galaxy clusters with double radio relics. Column2: Redshift of cluster. Column3: X-ray temperature of cluster. Column4: Projected distance of merge shock from the cluster center. Column5: Integrated radio spectral index. Column6: Radio Mach number calculated with α int . Column7: X-ray Mach number derived from temperature jump. Column8: X-ray Mach number derived from density jump. The # symbol indicates that M radio is calculated using α int given in the reference. The * symbol means that M radio can not be estimated, because α int < 1. The † symbol indicates that ⟨ kT X ⟩ is estimated using the L X -T relation in Pratt et al. (2009). For cases where shocks are not detected in X-rays or X-ray Mach numbers are not reported, n/a is indicated.  \nReferences -Abell 3376: Akamatsu & Kawahara (2013), George et al. (2015), Urdampilleta et al. (2018); Abell 3667: Akamatsu & Kawahara (2013), Sarazin et al. (2016), Storm et al. (2018), de Gasperin et al. (2022); Abell 3365: Duchesne et al. (2021a), Urdampilleta et al. (2021); ZwCl 0008.8+5215: van Weeren et al. (2011b), Kierdorf et al. (2017), Di Gennaro et al. (2019); Abell 3186: Nesci & Norci (1997), Duchesne et al. (2021a); Abell 1240: Barrena et al. (2009), Hoang et al. (2018), Sarkar et al. (2024); Abell 2345: George et al. (2017), Stuardi et al. (2021); 8C 0212+703: Hoang et al. (2021), Tumer et al. (2023); CIZA J2242.8+5301: Akamatsu & Kawahara (2013), Akamatsu & Kawahara (2013), Di Gennaro et al. (2018), Loi et al. (2020); Abell 2146: Russell et al. (2011), Russell et al. (2012), Hoang et al. (2019); RXC J1314.4-2515: Valtchanov et al. (2002), Venturi et al. (2007), George et al. (2017), Stuardi et al. (2019), Botteon et al. (2020); ZwCl 2341.1+0000 : Zhang et al. (2021), Stuardi et al. (2022); SPT-CL J2032-5627: Bulbul et al. (2019) Duchesne et al. (2021b) ; PSZ1 G096.89+24.17: Finner et al. (2021), Jones et al. (2021); PSZ1 G108.18-11.53: de Gasperin et al. (2015); PSZ2 G233.68+36.14: Ghirardini et al. (2021); MACS J1752.0+4440: van Weeren et al. (2012), Finner et al. (2021); ZwCl 1447.2+2619: Lee et al. (2022); PLCK G287.0+32.9: Bagchi et al. (2011), George et al. (2017); El Gordo: Menanteau et al. (2012), Botteon et al. (2020), Stuardi et al. (2022)  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 10. (a) Comparison of ⟨ M s ⟩ CR (black circles) of simulated merger shocks in our sample clusters with M radio (red squares) of observed radio relics in Table 2, plotted as a function of T X . (b) Comparison of ⟨ M s ⟩ X (black circles) of simulated merger shocks with M X -ray (red and blue squares) of observed radio relics, plotted as a function of T X . (c) and (d) Plots similar to (a) and (b), but shown as a function of the radio relic's distance from the cluster center, d s . Here, for simulated radio relics, d s is the mean 3D distance of the shock surface, ⟨ d s ⟩ . In contrast, for observed radio relics, d s refers to the distance 'projected' onto the sky, listed in Table 2. The vertical lines denotes the errors in observation.  \n<!-- image -->", '4. COMPARISON WITH OBSERVATIONS': "In this section, we attempt to compare the properties of merger shocks reproduced in our simulations with those inferred from observed double radio relics. So far, about a hundred clusters hosting radio relics have been identified, and of those, approximately two dozen have double relics (see, e.g., van Weeren et al. 2019). Table 2 provides a summary of the observational data and the Mach numbers inferred from radio and X-ray observations for the double radio relics reported in the literature (see Stuardi et al. 2022, and the references listed below the table). Here, M radio is the radio Mach number calculated using the integrated radio spectral index, α int , while the X-ray Mach number, M X -ray , is estimated from either the temperature jump, T jump , or the density jump, ρ jump , across the relic shock. And d s is the distance from the cluster center, projected onto the sky.  \nFigure 11. ⟨ M s ⟩ CR versus ⟨ M s ⟩ X (black circles) for simulated relic shocks in our sample clusters, and M radio versus M X -ray ( T jump ) (red square) for observed radio relics in Table 2. The vertical and horizontal lines denotes the errors in observation. The tilted line draws x = y .  \n<!-- image -->  \nIn Figure 10, we plot ⟨ M s ⟩ CR and ⟨ M s ⟩ X as functions of the X-ray temperature, T X , and the average distance from the center, ⟨ d s ⟩ , for our 24 simulated merger shocks, along with M radio and M X -ray as functions of T X and the 'projected' distance, d s , for the observed radio relics listed in Table 2. We compare ⟨ M s ⟩ CR with M radio , assuming that ⟨ M s ⟩ CR serves as a proxy of radio-emissivity-weighted average Mach number, since radio synchrotron emissions are produced by CR electrons accelerated via diffusive shock acceleration (DSA). Similarly, we compare ⟨ M s ⟩ X with M X -ray , noting that ⟨ M s ⟩ X is derived from 3D shock properties in our simulations, whereas M X -ray is estimated from the observed X-ray surface brightness map projected onto the sky.  \nA couple of points are notable. First, Figures 10(a) and (b) show that T X for simulated merger shocks (filled black circles) is, on average, lower than T X for observed radio relics (filled red and blue squares). This discrepancy arises because the merged clusters in our simulations have relatively low T X ≲ 4 keV due to the limited volume of our simulation box, as pointed in Section 2.2. In contract, X-ray and radio observations preferentially pick up clusters with heavier mass and higher T X . Hence, we present results in dimensionless quantities that would be roughly scale-independent, such as M s and ⟨ d ⟩ /r vir . Second, Figures 10(c) and (d) indicate that ⟨ d s ⟩ for simulated merger shocks is, on average, larger than d s for observed radio relics. This should be partly because the distances projected onto  \nthe sky are less than the true 3D distances for observed radio relics. We point out that ⟨ d s ⟩ ∼ r vir for simulated merger shocks, although there is considerable scatter, as shown in Figure 4(b); radio relics are also observed around r vir , as discussed in the introduction.  \nGiven the differences between simulated and observed quantities, the comparison should proceed with caution. Nevertheless, we point out the following: (1) The simulations predict ⟨ M s ⟩ CR ∼ 3 -5 for merger shocks, whereas the observed relics exhibit M radio ∼ 2 -5. Considering the limitations inherent in our simulations, this level of agreement would be considered reasonable. On the other hand, it was argued that the critical Mach number for electron preacceleration and subsequent injection to DSA is M crit ≈ 2 . 3 (e.g., Kang et al. 2019; Ha et al. 2021, 2022). Several observed relics have low M radio close to ∼ 2; however, these relics likely have regions with M s ≳ M crit where electrons DSA may occur. (2) In contrast, ⟨ M s ⟩ X ∼ 2 -4 for simulated merger shocks looks greater than M X -ray ∼ 1 . 2 -4 for observed relics in Figure 10. This discrepancy may partly arise from projection effects in X-ray observations (Hong et al. 2015; Hoang et al. 2018). Moreover, Wittor et al. (2021) demonstrated that M X -ray derived from the X-ray surface brightness map projected onto the sky plane can vary significantly depending on the relic's orientation.  \nIn Figure 11 we show the relation between the radio and X-ray Mach numbers for both simulated merger shocks and observed radio relics. In both samples, M radio ≳ M X -ray , confirming the earlier interpretation in Section 3.4 that electron acceleration is more efficient in portions of shock surfaces with higher M s , and consequently, M radio tends to be more heavily weighted by those higher M s regions. However, there are exceptions, such as Abell 3667 (NW) and ZwCl 0008.8+5215 (W), where M radio ≲ M X -ray ; these cases require more careful interpretation of radio and X-ray Mach numbers.", '5. SUMMARY': 'A major merger between two subclusters of comparable masses is expected to produce a pair of megaparsecscale bow shocks that propagate into the outskirts of the merged cluster over gigayears following pericenter passage (e.g., van Weeren et al. 2011a). When projected onto the sky, these merger shocks often appear as double radio relics in radio observations or as temperature and surface brightness discontinuities in X-ray observations (see Table 2).  \nIn this paper, we analyzed the data from a set of cosmological hydrodynamic simulations to investigate the properties of such shocks induced by major binary mergers. We focused on 12 merging clusters with mass ra-  \ntio 1 ≤M 1 / M 2 ≲ 4 and normalized impact parameter b/r , 1 ≲ 1, resulting in a total of 24 merger shocks examined at the optimal redshift for radio relic observation, z relic (see Table 1). Relic shock 1 forms in front of the heavier subcluster 1, while relic shock 2 forms in front of the lighter subcluster 2. We isolated the shock zones associated with the shock surfaces, distinguishing them from background ICM shocks produced by turbulence and infall, based on the criteria outlined in Section 2.3.  \nWe inspected the morphological characteristics of these merger shock surfaces, including their shape and area. We then quantified the Mach number distribution across the shock surfaces. Since merger shocks consist of numerous zones with varying Mach numbers, we defined representative average Mach numbers. We explored how these representative Mach numbers manifest in simulated merger shocks and compared them with the Mach numbers inferred in X-ray and radio observations.  \nOur main results are summarized as follows:  \n- 1. In our sample, merger shock surfaces are located at ⟨ d s ⟩ /r vir ≈ 0 . 6 -1 . 2 from the X-ray center of the clusters at z relic (see Figure 4). The shock surfaces can be approximated as elliptical sections of spherical shells, with an axial ratio of b ss /a ss ≳ 0 . 6. Each shock surface covers approximately ∼ 5 -20 % of the surface area of the shock sphere with r sphere = ⟨ d s ⟩ . The surface area ratio of relic shock 2 to relic shock 1 scales roughly with the mass ratio of subclusters as A ss , 2 /A ss , 1 ∝ M 1 / M 2 (see Figure 5).\n- 2. Due to nonzero impact parameters and turbulent ICM flows, merger shock surfaces become distorted, with their centers offset from the axis connecting the two DM density peaks (see Figures 6(a) and A2). As a result, the distribution of the average Mach number of shock zones, as a function of the position angle θ between the axis connecting the two DM density peaks and the line extending from the mid-point of these peaks to the shock zone, deviates from the typical cos θ pattern. Instead, these distributions often peak at fairly large θ , or exhibit significant fluctuations without a distinct peak (see Figures 6(b) and A3).\n- 3. The PDFs for M s of shock zones across merger shock surfaces are positively skewed and well-fitted by a log-normal function (see Figures 7(a) and A4). The fitted log-normal distributions peak around M s, peak ≈ 2 -4 . 5 with tails extending up to M s ∼ 10. The ratio of the standard deviation to the mean of the PDFs is approximately σ M s / ⟨ M s ⟩ area ≈ 0 . 37 (see Figures 7(b)). Here, the mean of the PDF is the shock-surface-areaweighted average Mach number, ⟨ M s ⟩ area .\n- 4. The area-weighted average Mach numbers, ⟨ M s ⟩ area ≈ 2 . 3 -4 . 4 are comparable to the X-ray-  \nvity-weighted average Mach numbers, ⟨ M s ⟩ X ≈ 2 -4. In contrast, the CR-energy-flux-weighted average Mach numbers are higher with ⟨ M s ⟩ CR ≈ 3 -5 (see Figure 8). This discrepancy aligns with observations of radio relics, where generally the Mach numbers inferred from radio data are greater than those derived from Xray data, M radio ≳ M X -ray .  \n- 5. In our simulated merger shocks, both ⟨ M s ⟩ X and ⟨ M s ⟩ CR exhibit a weak, decreasing trend with increasing M 1 / M 2 , albeit with significant scatter. On the contrary, the dependence of these Mach numbers on b/r vir , 1 appears to be only marginal (see Figure 9).\n- 6. For simulated merger shocks, we find ⟨ M s ⟩ X ≳ 2, while X-ray observations frequently report shocks with M X -ray ≲ 2 (see Figure 11). As suggested by several previous papers (e.g., Hong et al. 2015; Wittor et al. 2021), this discrepancy may be partly attributed to projection effects in X-ray observations, which can complicate accurate estimation of the true Mach numbers.  \nIn this paper, we focused on major binary mergers with 1 ≤M 1 / M 2 ≲ 4 and b/r , 1 ≲ 1, which are expected to present a relatively simple merger scenario. How-  \ndue to the complex interplay of factors such as merger history and turbulence, in addition to merger parameters, the structure of merger shocks turns out to be fairly complicated. Our aim was to characterize the properties of shock surfaces, yet a larger sample of simulated merging clusters appears necessary to establish more definitive relationships, such as the dependence of average Mach numbers on merger parameters.  \nFinally, some observed radio relics appear to have formed under merger conditions different from those considered here and also may result from multiple merger events. We leave investigations of these cases for future work.  \nThe authors would like to thank the anonymous referee for constructive comments and suggestions. This work was supported by the National Research Foundation (NRF) of Korea through grants 2020R1A2C2102800, 2023R1A2C1003131, RS-202200197685, and 2023K2A9A1A01099513. We thank Dr. J.-H. Ha for discussions.', 'REFERENCES': "Akamatsu, H., & Kawahara, H. 2013, PASJ, 65, 16, doi: 10.1093/pasj/65.1.16  \nBagchi, J., Sirothia, S. K., Werner, N., et al. 2011, ApJL, 736, L8, doi: 10.1088/2041-8205/736/1/L8 Barrena, R., Girardi, M., Boschin, W., & Das'ı, M. 2009, A&A, 503, 357, doi: 10.1051/0004-6361/200911788 Botteon, A., Brunetti, G., Ryu, D., & Roh, S. 2020, A&A, 634, A64, doi: 10.1051/0004-6361/201936216 Bruggen, M., Bykov, A., Ryu, D., & Rottgering, H. 2012, SSRv, 166, 187, doi: 10.1007/s11214-011-9785-9 Brunetti, G., & Jones, T. W. 2014, International Journal of Modern Physics D, 23, 1430007, doi: 10.1142/S0218271814300079 Bulbul, E., Chiu, I. N., Mohr, J. J., et al. 2019, ApJ, 871, 50, doi: 10.3847/1538-4357/aaf230 Dasadia, S., Sun, M., Sarazin, C., et al. 2016, ApJL, 820, L20, doi: 10.3847/2041-8205/820/1/L20 de Gasperin, F., Intema, H. T., van Weeren, R. J., et al. 2015, MNRAS, 453, 3483, doi: 10.1093/mnras/stv1873 de Gasperin, F., Rudnick, L., Finoguenov, A., et al. 2022, A&A, 659, A146, doi: 10.1051/0004-6361/202142658 Di Gennaro, G., van Weeren, R. J., Hoeft, M., et al. 2018, ApJ, 865, 24, doi: 10.3847/1538-4357/aad738 Di Gennaro, G., van Weeren, R. J., Andrade-Santos, F., et al. 2019, ApJ, 873, 64, doi: 10.3847/1538-4357/ab03cd  \nDominguez-Fernandez, P., Bruggen, M., Vazza, F., et al. Duchesne, S. W., Johnston-Hollitt, M., & Bartalucci, I. 2021a, PASA, 38, e053, doi: 10.1017/pasa.2021.45  \n2021, MNRAS, 500, 795, doi: 10.1093/mnras/staa3018 Duchesne, S. W., Johnston-Hollitt, M., Bartalucci, I., Hodgson, T., & Pratt, G. W. 2021b, PASA, 38, e005, doi: 10.1017/pasa.2020.51 Eckert, D., Jauzac, M., Vazza, F., et al. 2016, MNRAS, 461, 1302, doi: 10.1093/mnras/stw1435 Feretti, L., Giovannini, G., Govoni, F., & Murgia, M. 2012, The Astronomy and Astrophysics Review, 20, 1. https://api.semanticscholar.org/CorpusID:119258284 Finner, K., Jee, M. J., Cho, H., et al. 2024, arXiv e-prints, arXiv:2407.02557, doi: 10.48550/arXiv.2407.02557 Finner, K., HyeongHan, K., Jee, M. J., et al. 2021, ApJ, 918, 72, doi: 10.3847/1538-4357/ac0d00 Gabici, S., & Blasi, P. 2003, ApJ, 583, 695, doi: 10.1086/345429 George, L. T., Dwarakanath, K. S., Johnston-Hollitt, M., et al. 2015, MNRAS, 451, 4207, doi: 10.1093/mnras/stv1152 -. 2017, MNRAS, 467, 936, doi: 10.1093/mnras/stx155 Ghirardini, V., Bulbul, E., Hoang, D. N., et al. 2021, A&A, 647, A4, doi: 10.1051/0004-6361/202039554 Giacintucci, S., Venturi, T., Bardelli, S., et al. 2006, NewA, 11, 437, doi: 10.1016/j.newast.2005.11.001  \nStorm, E., Vink, J., Zandanel, F., & Akamatsu, H. 2018, MNRAS, 479, 553, doi: 10.1093/mnras/sty1435 Stuardi, C., Bonafede, A., Lovisari, L., et al. 2021, MNRAS, 502, 2518, doi: 10.1093/mnras/stab218 Stuardi, C., Bonafede, A., Rajpurohit, K., et al. 2022, A&A, 666, A8, doi: 10.1051/0004-6361/202244179 Stuardi, C., Bonafede, A., Wittor, D., et al. 2019, MNRAS, 489, 3905, doi: 10.1093/mnras/stz2408 Tumer, A., Wik, D. R., Zhang, X., et al. 2023, ApJ, 942, 79, doi: 10.3847/1538-4357/aca1b5 Urdampilleta, I., Akamatsu, H., Mernier, F., et al. 2018, A&A, 618, A74, doi: 10.1051/0004-6361/201732496 Urdampilleta, I., Simionescu, A., Kaastra, J. S., et al. 2021, A&A, 646, A95, doi: 10.1051/0004-6361/201937160 Valtchanov, I., Murphy, T., Pierre, M., Hunstead, R., & L'emonon, L. 2002, A&A, 392, 795, doi: 10.1051/0004-6361:20020940 van Weeren, R. J., Bruggen, M., Rottgering, H. J. A., & Hoeft, M. 2011a, MNRAS, 418, 230, doi: 10.1111/j.1365-2966.2011.19478.x van Weeren, R. J., de Gasperin, F., Akamatsu, H., et al. 2019, SSRv, 215, 16, doi: 10.1007/s11214-019-0584-z van Weeren, R. J., Hoeft, M., Rottgering, H. J. A., et al. 2011b, A&A, 528, A38, doi: 10.1051/0004-6361/201016185  \nvan Weeren, R. J., Rottgering, H. J. A., Bruggen, M., & Hoeft, M. 2010, Science, 330, 347,  \ndoi: 10.1126/science.1194293  \nvan Weeren, R. J., Rottgering, H. J. A., Intema, H. T., et al. 2012, A&A, 546, A124,  \ndoi: 10.1051/0004-6361/201219000  \n- van Weeren, R. J., Brunetti, G., Bruggen, M., et al. 2016, ApJ, 818, 204, doi: 10.3847/0004-637X/818/2/204  \nVazza, F., Bruggen, M., Gheller, C., & Wang, P. 2014, MNRAS, 445, 3706, doi: 10.1093/mnras/stu1896  \nVenturi, T., Giacintucci, S., Brunetti, G., et al. 2007, A&A, 463, 937, doi: 10.1051/0004-6361:20065961  \n- Wittor, D., Ettori, S., Vazza, F., et al. 2021, MNRAS, 506, 396, doi: 10.1093/mnras/stab1735\n- Zhang, C., Churazov, E., Forman, W. R., & Lyskova, N. 2019, MNRAS, 488, 5259, doi: 10.1093/mnras/stz2135\n- Zhang, X., Simionescu, A., Stuardi, C., et al. 2021, A&A, 656, A59, doi: 10.1051/0004-6361/202141540  \nZuHone, J. A. 2011, ApJ, 728, 54,  \ndoi: 10.1088/0004-637X/728/1/54", 'A. FIGURES OF ALL SAMPLE CLUSTERS': 'In the main part, we present figures for Clusters 1, 3, 6, and 11 to maintain conciseness. In this appendix, we provide figures for all 12 sample clusters for comprehensive reference.  \nFigure A1. Same as in Figure 4(a) for all 12 sample clusters.  \n<!-- image -->  \nFigure A2. Same as in Figure 6(a) for all 12 sample clusters.  \n<!-- image -->  \nlog  \n𝐷𝐷𝐷𝐷  \n𝜌𝜌  \n𝜃𝜃  \n(°)  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝜃𝜃  \n(°)  \n𝜃𝜃  \n(°)  \nFigure A3. Same as in Figure 6(b) for all 12 sample clusters.  \n<!-- image -->  \n𝜃𝜃  \n(°)  \n𝜃𝜃  \n(°)  \n𝜃𝜃  \n(°)  \n𝜃𝜃  \n(°)  \n)  \n𝑠𝑠  \n𝑀𝑀  \n(  \n𝑃𝑃𝑃𝑃  \n𝑃𝑃  \n)  \n𝑠𝑠  \n𝑀𝑀  \n(  \n𝑃𝑃𝑃𝑃  \n𝑃𝑃  \n)  \n𝑠𝑠  \n𝑀𝑀  \n(  \n𝑃𝑃𝑃𝑃  \n𝑃𝑃  \n)  \n𝑠𝑠  \n𝑀𝑀  \n(  \n𝑃𝑃𝑃𝑃  \n𝑃𝑃  \nFigure A4. Same as in Figure 7(a) for all 12 sample clusters.  \n<!-- image -->  \n0.12  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀  \n𝑠𝑠  \n𝑀𝑀'}