Exploring the role of cosmological shock waves in the Dianoga simulations of galaxy clusters

The Dianoga set of simulations analyzed in this work has been performed with an improved version of the TreePM-SPH code GADGET-3 (Springel, 2005). This version of the code includes an updated hydrodynamical scheme that has been shown to ameliorate the performance of standard SPH algorithms in a number of issues (see Beck et al., 2016b, for details). Here we only provide a brief description of the simulations, while we defer the interested reader to previous works where different aspects of the same set of simulations were investigated (Rasia et al., 2015; Truong et al., 2018; Villaescusa-Navarro et al., 2016; Planelles et al., 2017; Biffi et al., 2017, 2018; Ragone-Figueroa et al., 2018; Truong et al., 2019).

The set of simulations consists in re-simulations of 29 Lagrangian regions extracted from a larger parent dark matter (DM) only simulation (see Bonafede et al., 2011, for details on the initial conditions). The 29 regions were identified at $z=0$ around 29 massive dark matter haloes with $M_{200}$²²2$M_{\Delta}$ is the mass within $R_{\Delta}$, which is the radius of a sphere enclosing an average density equal to $\Delta$ times the critical cosmic density, $\rho_{c}(z)$. In the following, we will refer to the virial overdensity (Bryan and Norman, 1998), that is $\Delta_{vir}\simeq 93$ at $z=0$ for our cosmology, or to $\Delta=180,200,500,2500$.$\sim 1-30\times 10^{14}h^{-1}\,M_{\odot}$. Each of the low-resolution regions has been re-simulated at a finer resolution and by incorporating the baryonic component (see Planelles et al., 2013, for details). The simulations are based on a flat $\Lambda$CDM cosmology with $\Omega_{\rm{m}}=0.24$, $\Omega_{\rm{b}}=0.04$, $H_{0}=$ 72 km s^-1 Mpc^-1, $\sigma_{8}=0.8$, and $n_{\rm{s}}=0.96$. For the Dianoga simulations analyzed in this paper, the mass resolution for the DM (gas) particles is $m_{\rm{DM}}=8.44\times 10^{8}\,\rm{\,M_{\odot}}$ ($m_{\rm{gas}}=1.53\times 10^{8}\,\rm{\,M_{\odot}}$). In the high-resolution region the gravitational force is computed with a Plummer-equivalent softening length of $\epsilon=3.75\,h^{-1}$ kpc for DM and gas and $\epsilon=2\,h^{-1}$ kpc for black hole and stellar particles.

These simulations have been run with three different prescriptions for the physics of baryons. The most complete version of these simulations, labelled as AGN, accounts for the effects of a number of gas physical processes such as gas radiative cooling, star formation, SNe feedback, metal enrichment and a novel AGN feedback scheme that accounts for both cold and hot gas accretion onto supermassive black holes (Steinborn et al., 2015). There are two additional sets of simulations that only differ with the AGN runs in the baryonic physics they include: the NR runs, that only account for non-radiative physics, and the CSF runs, that are like the AGN ones but without including the effect of AGN feedback.

In the following, unless otherwise specified, we will consider the AGN simulation as our reference set. The set of AGN cosmological simulations has been shown to be extremely successful in reproducing some of the main observational cluster properties. These simulations have produced, for the first time and for a statistically significant cluster sample, the CC/NCC dichotomy according to thermal and chemical cluster core properties (see Rasia et al., 2015). They have also shown an excellent agreement with cluster observations in terms of the X-ray and SZ scaling relations (both at low and high redshift; Planelles et al., 2017; Truong et al., 2018), the amount of hosted neutral hydrogen (Villaescusa-Navarro et al., 2016), the distribution of thermodynamical properties such as entropy, thermal pressure and temperature (Rasia et al., 2015; Planelles et al., 2017), the level of mass bias and hydrostatic equilibrium (Biffi et al., 2016), or the ICM chemical enrichment distribution (Biffi et al., 2017, 2018; Truong et al., 2019).

Overall, the AGN simulation comprises a sample of $\sim 100$ clusters and groups with $M_{500}>3\times 10^{13}h^{-1}M_{\odot}$ at $z=0$. However, in this work we will only consider the sample of the main 29 central haloes at $z=0$, that is, 24 clusters with $M_{200}>8\times 10^{14}h^{-1}M_{\odot}$ and 5 isolated groups with $M_{200}$ in the range $1-4\times 10^{14}h^{-1}M_{\odot}$.

As described in Rasia et al. (2015), according to the core thermal properties of the clusters, the main 29 objects in the AGN simulations were classified, at $z=0$, in CC and NCC systems. In particular, those systems with a central entropy $K_{0}<60$ keV cm² and a pseudo-entropy $\sigma<0.55$ were considered as CCs, whereas those systems not fulfilling these conditions were instead classified as NCCs. In this way, 11 out of 29 halos are classified as CC clusters. The samples of CC and NCC clusters have been shown to be different in terms of metallicity and thermodynamical profiles, showing a remarkably consistency with observations of these two cluster populations. In a similar way, based on the global dynamical state of the 29 main systems in the AGN simulations, we have classified them in 6 regular, 8 disturbed and 15 intermediate systems (e.g. Biffi et al., 2016; Planelles et al., 2017). In order to do this classification we have considered a system to be regular when $\delta r<0.07$ and $f_{sub}<0.1$, where $\delta r$ and $f_{sub}$ represent, respectively, the centre shift between the cluster minimum potential and the cluster centre of mass, in units of its radius $R_{200}$, and the mass fraction contributed by substructures within the same radius. Systems with larger values for both $\delta r$ and $f_{sub}$ are classified as disturbed, whereas those systems not satisfying concurrently both criteria are labeled as intermediate cases. The samples of dynamically regular and disturbed systems have been also shown to differ in terms of the levels of hydrostatic equilibrium and gas clumpiness, especially in outer cluster regions (see Biffi et al., 2016; Planelles et al., 2017, for further details).

Note that the classification of clusters in dynamically regular or disturbed depends on the radius at which the dynamical state is defined (e.g. De Luca et al., 2021, and references therein). Therefore, although the choice of this radius is quite arbitrary, it is important to keep this aspect in mind when discussing our results about the subsamples of dynamically regular and disturbed systems.

In general, the development of a shock in a cosmological simulation generates a discontinuity in all the hydrodynamical gas quantities (namely, gas density, temperature, pressure and entropy). The pre- and post-shock values of these quantities are connected to the strength of the shock through the Mach number, $\mathcal{M}=v_{s}/c_{s}$, where $v_{s}$ and $c_{s}$ are the shock speed and the sound speed ahead of the shock, respectively.

Numerical methods developed to identify hydrodynamical shocks in cosmological simulations commonly rely on the Rankine-Hugoniot jump conditions (Landau and Lifshitz, 1966), which provide all the necessary information to calculate the shock Mach number. However, depending on the particular quantity and the particular way in which these conditions are applied, a number of shock-finding algorithms have been suggested. Thus, attending to the particular jump condition from which the shock Mach number is ultimately derived, there are methods based on the temperature, the density, the entropy or the velocity jump condition (e.g. Vazza et al., 2009a). Likewise, depending on the way in which the direction of shock propagation is defined, there are two main numerical approaches: the coordinate or directional-splitting methods, that follow the jump in a given hydrodynamical quantity along the coordinate axes of the computational domain (e.g. Ryu et al., 2003; Vazza et al., 2009b; Planelles and Quilis, 2013; Hong et al., 2014, 2015), and the methods based on the local temperature or pressure gradients (e.g. Skillman et al., 2008; Schaal and Springel, 2015; Schaal et al., 2016; Beck et al., 2016a). Despite the obvious and unavoidable deviations between different shock-finding schemes and the distinct nature and properties of the cosmological simulations on which they are applied, all algorithms tend to produce similar properties of the shock waves in the LSS (see, e.g., Vazza et al., 2011b, for a comparison study with different cosmological codes).

In this work, in order to analyze the distribution, the evolution and the strength of shock waves in the simulations described in Section 2.1, we will employ the shock-finding algorithm presented in Planelles and Quilis (2013) (see also Martin-Alvarez et al., 2017, for further details). The definition of cells for these Lagrangian simulations will be described in Section 2.3. Here we only outline the main procedure of the algorithm and defer the interested reader to the previous references for further details.

The shock finder algorithm presented in Planelles and Quilis (2013) was developed to measure the strength of shocks in grid-based cosmological simulations. The code relies on a directional-splitting approach along the three coordinate axes and evaluates the shock Mach number from the temperature jump equation:

where $T_{1}$ and $T_{2}$ refer to the pre- and post-shock temperatures.

In our scheme, we first mark grid cells as tentative shocked if they fulfill the following conditions:

where (i) ensures the local gas velocity field ${\bf v}$ to be convergent and (ii) guarantees the same direction for the gradients of gas temperature $T$ and entropy $S$. Then, among all the tentative shocked cells, the cell where $\nabla\cdot{\bf v}$ is minimum is identified as the first shock centre. Moving outwards from the shock centre, we define the extension of the shock along each of the three coordinate axes by looking for adjacent tentative shocked cells where the pre- and post-shock temperature and density satisfy $T_{2}>T_{1}$ and $\rho_{2}>\rho_{1}$. Once the furthest shocked cells ahead (pre-shock) and behind (post-shock) of the shock centre along each coordinate axis are identified, the Mach number along each coordinate direction is computed substituting the corresponding temperatures $T_{1}$ and $T_{2}$ in Eq. 1. Following this approach shock discontinuities are spread, typically, over a few cells. The Mach numbers obtained along the three coordinate directions are then combined to get the final strength of the first shocked cell: $\mathcal{M}=(\mathcal{M}_{x}^{2}+\mathcal{M}_{y}^{2}+\mathcal{M}_{z}^{2})^{1/2}$, thus minimizing projection effects in case of diagonal shocks. After this step, this shocked cell is removed from the list of tentative shocked cells. Then, the whole process is iteratively repeated focusing on the cell with the minimum value of $\nabla\cdot{\bf v}$.

With this procedure we can identify all the shocked cells within the computational volume, obtaining their Mach number and the typical shock surfaces associated to cosmological shock waves. As an additional condition, in order to prevent noisy shock regions, it is a common practice to neglect all shocked cells with a Mach number lower than 1.3 (e.g. Ryu et al., 2003; Planelles and Quilis, 2013; Schaal and Springel, 2015).

The Rankine-Hugoniot jump conditions also provide an estimation of the amount of kinetic energy flux crossing a shock ($f_{K}$) converted into thermal energy flux ($f_{th}$) in the post-shock region³³3Note that conversion of kinetic to thermal energy at shocks in SPH simulations is described by the prescription of artificial viscosity, which not necessarily coincides with the Rankine-Hugoniot jump conditions.. The incoming kinetic energy flux at a shock is given by $f_{K}={\frac{1}{2}\rho_{1}(c_{s,1}\mathcal{M})^{3}}$, where $\rho_{1}$ and $c_{s,1}$ are the gas density and the sound speed in the pre-shock region, respectively. Written in this way, $f_{K}$ has units of energy flux, that is energy per unit time and per unit area. Therefore, the kinetic energy flux across shock surfaces could be computed as $F_{K}=f_{K}(\Delta x)^{2}$, being $\Delta x$ the spatial grid resolution. Given the kinetic energy flux, the generated thermal energy flux can be estimated as follows:

where $\delta({\mathcal{M}})$ is the gas thermalization efficiency at shocks. Following Kang et al. (2007) for the case without a pre-existing CR component, $\delta({\mathcal{M}})$ can be computed as a function of the Mach number:

where $\gamma$ is the adiabatic exponent and $R$ depends on the jump in density:

Besides gas thermalization, part of the energy flux at a shock can also be transferred into cosmic rays. Indeed, by means of the DSA mechanism, cosmic rays can also be accelerated and injected at shocks. The efficiency of this injection depends on the shock Mach number, on the degree of Alfvén turbulence and on the direction of the magnetic field relative to the shock surface. Kang et al. (2007) inferred limits for the efficiency of CR acceleration, $\eta({\mathcal{M}})$, providing fitting functions for the cases with and without a pre-existing CR component. In our analysis, once shocks are identified we will use these fitting formulae, for the case without a pre-existing CR component, to estimate the flux of energy transferred into CRs:

Note that the treatment of energy dissipation at shocks by means of an approximated post-processing approach, like the one described above, presents a number of simplifications and caveats. However, it will allow us to get some upper limits for the efficiency of energy dissipation at shocks. It is also important to keep in mind that the Rankine-Hugoniot jump conditions provide a description of the strength of shock waves in the case of ideal shocks. Since shocks in simulations are far from being ideal, any numerical procedure based on these equations will suffer from inevitable uncertainties.

Given the characteristics of our shock-finder algorithm, in order to obtain the distribution of shock waves and their associated Mach numbers, the code needs to be applied onto a computational grid where the gas main quantities are stored in cells of any given resolution. Therefore, before running the shock-finder we need to convert all the information provided by our simulations from the distribution of particles to a proper computational grid. In order to compute the Mach number within different radial apertures from the cluster centre we proceed in the following way. For each selected simulated cluster, with virial radius $R_{vir}$, we build a cubic box, with centre coinciding with the cluster centre, and with a side length of $8\times R_{vir}$. The box is then discretized with a given number of cubical cells ($N=128^{3},256^{3}$, 5$12^{3}$ cells). On average, for the whole sample of 29 clusters in the AGN run at $z=0$, when we consider a number of grid cells equal to $N=128^{3},256^{3}$, $512^{3}$ we obtain, respectively, a mean spatial resolution of $\Delta x\simeq 144,72$ and $36\,\,h^{-1}\rm{\,kpc}$. Then, we select all the gas particles within a distance of $4\times R_{vir}$ from the halo centre and we interpolate the main gas properties onto the regular grid (see Röttgers and Arth, 2018, for a recent comparison of different smoothing methods) by spreading them accordingly to the same Wendland C4 kernel used for the SPH computations in the simulations (see Beck et al., 2016b). After this, we apply the shock finder in order to compute the shock Mach number for several radial apertures from the cluster centre. In particular, for each halo we have considered 4 different radial apertures: $R_{core}\equiv 0.05\times R_{180},R_{2500},R_{500}$ and $R_{vir}$, although throughout this work we will mainly focus on $R_{500}$ and $R_{vir}$. Once this procedure is applied to the main 29 clusters, we analyze the distribution of shock Mach numbers as a function of redshift, cluster mass, cluster cool-coreness, cluster dynamical state, radial aperture, spatial grid resolution and physics included in the simulations.

Unless otherwise specified, in the following we will employ as our reference results those obtained for the AGN simulations when smoothing each computational domain onto a regular grid with $128^{3}$ cells. In Appendix A we explore the effects of grid resolution on our results and we justify our choice. The selected spatial grid resolution is in line with the resolution employed in some recent works on the topic (e.g. Vazza et al., 2010, 2011b, 2013).

We start by analysing the global shock cell distribution around D21, one of the most massive clusters in our AGN simulation. Figure 1 shows 2-D projections along the z-axis of the gas overdensity (upper panels) and the Mach number distribution (lower panels) at different redshifts (panels from left to right show the projections at $z=2,1,0.5$ and 0, respectively). Each map has an extent of $8\times R_{vir}$ at each redshift and a projection depth of 7 cells centred on the cluster position. Cluster D21, with $M_{vir}=14.26\times 10^{14}\,h^{-1}\rm{\,M_{\odot}}$ and $R_{vir}=2.37\,h^{-1}\rm{\,Mpc}$, has been classified as dynamically relaxed and it is a CC cluster. As can be seen in the density maps, the cluster is already well identified at $z=2$. The mean Mach number within the whole region is around $\sim 4$ throughout the redshift evolution (see also Fig. 12). However, as evolution proceeds the Mach number of shocks around the central cluster increases, developing a significant high-Mach number accretion shock around the cluster centre at $z=0$. At the same time, as redshift decreases, the distribution of shocks becomes somewhat more volume-filling (see Fig. 12). Moreover, according to these maps, high-Mach number external shocks seem to be clearly dominant, in terms of strength and extension, over low-Mach number internal ones. To understand this behaviour it is interesting to note that accretion shocks from voids and low-density regions onto the cosmic web are very volume-filling. In addition, as evolution proceeds void regions expand and increase their size, making the distribution of low-Mach number shocks within them to appear somewhat more diluted. Concurrently, the temperature of the gas in voids decreases, producing a significant difference between the temperature in voids and that in collapsed regions and, therefore, generating stronger shocks around collapsed structures.

In Appendix C we compare the distribution at $z=0$ shown in Fig. 1 for cluster D21 with the same distribution as obtained in the NR and CSF simulations (see Fig. 13). On average, the shock configurations are similar with comparable Mach numbers and locations, confirming the fact that the high-Mach number external shocks shown in Fig. 1 are indeed accretion shocks, with the exception of small differences in the Mach number distribution between the NR and the AGN simulations, as discussed more in Appendix B. As inferred from previous simulations (e.g. Schaal and Springel, 2015; Schaal et al., 2016), outer accretion shocks show typically quasi-spherical shapes and are found at distances of $\sim 1.5-2\times R_{vir}$. According to Figs. 1 and 13 (see also Fig. 3), the accretion shock developed at $z=0$ around cluster $D21$ also shows a quasi-spherical pattern, it wraps the LSS around the cluster and it is located at $\sim 2.5\times R_{vir}$ from the cluster centre.

To have an idea of the way in which shocked cells of different strengths are distributed through the computational volume, in the left panel of Fig. 2 we show the mean shock cell distribution function (SDF, from now on) obtained for the whole sample of 29 clusters in the AGN simulation at different redshifts. From the inspection of this figure at $z=0$ (red continuous line), we can extract several broad conclusions. As shown in previous studies (e.g. Skillman et al., 2008; Vazza et al., 2009b; Planelles and Quilis, 2013), on average most of the computational volume is filled with low-Mach number shocks ($\mathcal{M}\raise-2.0pt\hbox{\hbox to0.0pt{\hbox{$\sim$}\hss}\raise 5.0pt\hbox{$<$}\ }3$, blue-violet colors in Fig. 1), which tend to be located within the virial radius. In contrast, stronger shocks with higher Mach numbers (orange-yellow colors in Fig. 1) are less numerous and tend to be found in outer cluster regions (see also Fig. 3). According to this behaviour, the SDF shows a peak at around $\mathcal{M}\sim 2$ and decreases beyond that value. As already discussed by some other authors (e.g. Ryu et al., 2003), internal shocks within clusters result from a number of processes such as substructure mergers, accretion of warm-hot gas along filaments, or internal turbulent flow motions. On the contrary, external shocks around clusters are associated to large-scale accretion phenomena, especially when the cold gas from lower-density regions accretes onto the hot ICM. The contrast, in terms of density and temperature, between shocked and unshocked regions within the different environments where internal and external shocks are developed makes them to have very different strengths.

As for the dependence of the SDF on redshift, we show that the general shape of the mean SDF remains quite stable already from $z=2$. However, the amount of shocked cells within the computational domain slightly increases with decreasing redshift, reaching a mean fraction of shocked cells of $\sim 13$ per cent at $z=0$. In this sense, as the redshift decreases the number of shocked cells within the considered regions augments, giving raise, however, to a mild increase in the amount of weak internal shocks. In a similar way, as the cosmic evolution proceeds, the temperature of the IGM, especially in low-density regions, decreases, explaining the increasing high-Mach number end of the SDF (e.g. Vazza et al., 2009b; Schaal and Springel, 2015; Schaal et al., 2016).

For completeness, the right panel of Fig. 2 shows the mean z-evolution of the distribution of the thermal energy flux through shocks for the sample of 29 clusters in the AGN simulation. Overall, at any redshift, the maximum of the gas thermalization by shocks is found at $\mathcal{M}\sim 2$, consistent with previous studies (e.g. Ryu et al., 2003; Pfrommer et al., 2006; Skillman et al., 2008; Vazza et al., 2009b). In addition, independently of the considered redshift, the shape of the thermal energy distribution is quite alike, showing a steep negative dependence with the Mach number. In broad agreement with previous results (e.g. Pfrommer et al., 2007; Vazza et al., 2009b), we find the slope of this distribution to be within $\alpha\simeq-4$ ($z=0$) and $\alpha\simeq-3$ ($z=2$). Despite the high Mach number of external shocks, given the low-density environment where they are usually developed, the associated kinetic energy flux is very small, making them energetically less important. On the contrary, low-Mach number shocks, mainly developed within collapsed and denser structures, have larger values of $f_{K}$ and, therefore, they are more relevant in thermalizing the ICM and producing CRs. As a consequence, most of the total thermal energy flux is processed by shocks with relatively low Mach numbers ($\mathcal{M}\raise-2.0pt\hbox{\hbox to0.0pt{\hbox{$\sim$}\hss}\raise 5.0pt\hbox{$<$}\ }10$; e.g. Vazza et al., 2011b).

In order to explore the radial distribution of shocks around the main clusters of each simulated region, Fig. 3 shows, at $z=0$, the averaged radial profiles of the mean Mach number for the whole sample of 29 objects in the AGN simulation (black continuous line). The mean radial profiles corresponding to the subsamples of CC/NCC and regular/disturbed systems are also shown for comparison (left and right panels, respectively).

In general, when we analyze the mean global profile, we find that the mean Mach number radial distribution is quite flat out to $\sim 0.8\times R_{vir}$, showing values within the range $\sim 2.2-2.5$. At radial distances above $\sim 0.8\times R_{vir}$ the profiles raise significantly, reaching values as high as $\mathcal{M}\sim 5$ in the outer cluster regions. This trend is in broad agreement with results obtained in previous studies (e.g. Vazza et al., 2009a, 2010, 2011b). In particular, Vazza et al. (2011b) performed a comparison of the radial Mach number distribution in two massive clusters as obtained with two grid-based cosmological codes and with the SPH code GADGET-3. Comparing with our results, they also found a mean Mach number value below $\sim 3$ inside the virial radius (see also, e.g. Ha et al., 2018) and a very similar radial shock cell distribution, at least within $\sim 0.5R_{vir}$. However, beyond the virial radius, they showed that grid codes tend to identify thinner and stronger shocks than GADGET-3, producing a much sharper increase of the profiles. This abrupt rise in the strength of shocks clearly marks the transition between the outskirts of clusters and the low-density IGM. The particular shape of the profiles in these outer cluster regions could be affected by a number of factors such as the simulation numerical scheme, the specific shock-detecting approach, the different cluster environment derived from the intrinsic difference between cosmological and zoom-in simulations, or the grid resolution and/or the gas sampling procedure (especially in outer and underdense cluster regions where, by construction, SPH is not very good at sampling). For all these reasons, a precise comparison between different simulations is not straightforward. Nevertheless, even with these weaknesses, results are in confortable agreement with those obtained using grid codes.

As for the connection between the distribution and strength of shocks and the cluster cool-coreness and global dynamical state, some previous studies attempted to analyze the relation between relaxed and unrelaxed (merging) clusters with the strength of their associated shock waves (e.g. Vazza et al., 2011a, 2012). Overall, these studies found similar shock cell distributions from cluster to cluster with only small deviations connected to their particular dynamical state, especially above the virial radius.

In our case, according to the results shown in the left panel of Fig. 3, we do not obtain any difference between the mean Mach radial profiles when we look separately to the samples of CC and NCC clusters. This result is not surprising, since the classification of clusters according to their cool-coreness is based on a much smaller radial aperture ($R_{core}\equiv 0.05\times R_{180}$; Rasia et al., 2015), which is completely smoothed out at the considered resolution. As discussed in Appendix A, given the numerical approach we follow to smooth the particles’ distribution onto a computational grid, it would be risky to artificially increase the spatial grid resolution to try to ‘resolve’ the core of clusters.

When we analyze instead the mean radial profiles obtained for the subsamples of relaxed and disturbed systems (right panel of Fig. 3) we obtain a larger difference between them: disturbed clusters clearly show a higher mean Mach number than regular ones throughout the radial range. A similar behaviour was also obtained for our sample of regular and disturbed systems when we compared their degree of deviation from the hydrostatic equilibrium condition or their levels of gas clumpiness (see Biffi et al., 2016; Planelles et al., 2017, respectively). The trend shown in the right panel of Fig. 3, which was somehow expected in our sample, is also in broad agreement with previous analyses (e.g. Vazza et al., 2009a).

Galaxy cluster scaling relations are extremely useful to estimate cluster masses from cluster observables. In the left panel of Fig. 4 we show the mean Mach number within $R_{vir}$ as a function of the cluster mass $M_{500}$ for the whole sample of 29 central clusters in the AGN simulation at redshifts $z=0,0.5,1$ and 2. The value of $\mathcal{M}(<R_{vir})$ is shown for the populations of dynamically regular, disturbed and intermediate systems as specified in the legend. Moreover, CC clusters are marked with an additional circle around the main symbols, being all the clusters without this circle NCC clusters. Similarly, in the right panel of Fig. 4 the mean thermal energy flux, $F_{th}$, within $R_{vir}$ is shown. Table 1 shows the mean values of these quantities at $z=0$ for the different cluster samples. For the sake of completeness, the mean energy flux inverted in CR acceleration is also given.

On average, when all redshifts are considered, the mean Mach number value within the virial radius is almost flat as a function of cluster mass, with a significant dispersion for all redshifts ($\mathcal{M}(<R_{vir})\sim[2.2,3.8])$. However, if we consider the whole cluster sample at $z=0$, the Pearson’s correlation coefficient between the Mach number values within the virial radius and the mass $M_{500}$ is $\sim 0.47$, indicating a moderate correlation. Instead, when we consider the subsamples of regular and disturbed systems, we obtain a positive correlation with values of the Pearson’s correlation coefficient above $\sim 0.7$. We have checked that this trend holds from $R_{200}$ and beyond. However, given that the considered mass range is quite narrow and the cluster sample is relatively limited, we should be cautious when deriving conclusions on this regard. As suggested by the profiles shown in Fig. 3, mild differences between dynamically regular and disturbed systems are observed throughout the radial range, contributing to a difference between the mean $\mathcal{M}(<R_{vir})$ of regular and disturbed systems which is not found between that of CC and NCC clusters (see also Table 1).

On the contrary, as shown in the right panel of Fig. 4, the mean thermal energy flux within $R_{vir}$ shows, at fixed redshift, some correlation with cluster mass although somewhat scattered. It is only combining different redshifts that a stronger correlation shows up. However, this seems to be due to the redshift dependence. Indeed, as discussed by other authors (e.g. Ryu et al., 2003), massive galaxy clusters are supposed to represent the main sites for energy dissipation at shocks. Therefore, since the total energy processed by shocks in a simulated volume depends on the number of clusters, we should be cautious when comparing these results with preceding works (e.g. Ryu et al., 2003; Pfrommer et al., 2006; Vazza et al., 2009a). According to the results shown in Table 1, at $z=0$ our sample of CC clusters have larger mean values of $F_{th,vir}$ than the NCC systems. In a similar way, as for the global dynamical state of clusters, disturbed objects show larger values of the mean thermal energy flux than regular ones. Since all these cluster subsamples have similar mean virial Mach numbers, differences could be also connected with a different amount of gas density contributing to the dissipation of energy.

Therefore, according to the global analysis of Fig. 4, our results suggest that, provided that we have a significant cluster statistics, these relations could be used as potential cluster mass proxies.

A complementary way to analyze the radial distribution of shocks around clusters is given in the left panel of Fig. 5, that shows the radial shock surface distribution, out to $4\times R_{vir}$, for the whole sample of clusters in the AGN run at $z=0$ (black continuous line). This quantity provides an approximated estimation of the surface of shock cells, $S$, at their locations, $r$. For the sake of comparison, we also report the mean shock surface distribution obtained for the subsamples of CC/NCC and regular/disturbed systems. As shown in previous studies (Schaal and Springel, 2015; Schaal et al., 2016), this distribution provides information on the location of the main shock surfaces within and around the cluster virial region and it is, therefore, very useful to estimate the position and shape of external accretion shocks. Within the virial radius, the different mean shock surface distributions show a number of peaks or bumps which make the profiles quite irregular. These internal peaks could be connected with a number of internal shock events, such as shocks developed as a consequence of merging substructures, flow motions or feedback processes. However, above $R_{vir}$ all mean profiles tend to smoothly increase out to $\sim 2.5-3\times R_{vir}$, where the curves show their maximum value. The position of this maximum could be connected with the radius (provided that accretion shocks are assumed to be spherical) of the main external accretion shock around the central clusters. It is interesting to note that there is a clear distinction between the mean shock surface distribution associated to regular and disturbed systems, with the latter showing the highest values throughout the radial range. On the contrary, we have not found such a strong dependence of $d(S/r^{2})/dR$ on the clusters’ core properties. Interestingly, we have obtained a moderate correlation (with a Pearson’s correlation coefficient of $0.6$) between the maximum value of the distribution and the mass of the associated systems. At $z=0$, by fitting linearly the relation max$(d(S/r^{2})/dR)=f(log(M_{500}))$ (see right panel of Fig. 5), we have obtained for the whole sample of clusters in the AGN run a slope $\alpha=1.1\pm 0.2$, although with some scatter⁴⁴4Similarly, we have obtained for the reduced sample of clusters in the CSF and NR simulations a value of $\alpha=1.1\pm 0.9$ and $\alpha=1.3\pm 0.9$, respectively.. This trend points towards a connection between the extension (or surface filling factor) of external accretion shocks developed during the collapse of structures and the mass of the final formed systems. Assuming that future observational radio facilities, such as the Square-Kilometre Array⁵⁵5https://www.skatelescope.org/ (SKA; e.g. Acosta-Pulido et al., 2015), will be able to routinely detect large-scale shock waves, correlations like the one presented here, which is rather independent on the physics included, could become a complementary tool to estimate cluster masses (e.g. Planelles and Quilis, 2013).