N IFT Y galaxy cluster simulations – VI. The dynamical imprint of substructure on gaseous cluster outskirts

Galaxy cluster outskirts mark the transition region from the mildly non-linear cosmic web to the highly non-linear, virialized, cluster interior. It is in this transition region that the intracluster medium (ICM) begins to inﬂuence the properties of accreting galaxies and groups, as ram pressure impacts a galaxy’s cold gas content and subsequent star formation rate. Conversely, the thermodynamical properties of the ICM in this transition region should also feel the inﬂuence of accreting substructure (i.e. galaxies and groups), whose passage can drive shocks. In this paper, we use a suite of cosmological hydrodynamical zoom simulations of a single galaxy cluster, drawn from the N IFT Y comparison project, to study how the dynamics of substructure accreted from the cosmic web inﬂuence the thermodynamical properties of the ICM in the cluster’s outskirts. We demonstrate how features evident in radial proﬁles of the ICM (e.g. gas density and temperature) can be linked to strong shocks, transient and short-lived in nature, driven by the passage of substructure. The range of astrophysical codes and galaxy formation models in our comparison are broadly consistent in their predictions (e.g. agreeing when and where shocks occur, but differing in how strong shocks will be); this is as we would expect of a process driven by large-scale gravitational dynamics and strong, inefﬁciently radiating, shocks. This suggests that mapping such shock structures in the ICM in a cluster’s outskirts (via e.g. radio synchrotron emission) could provide a complementary measure of its recent merger and accretion history.

The imprint of this cosmic web can be inferred from the spatial and kinematical distribution of cluster galaxies (e.g. Knebe et al. 2004;Bailin & Steinmetz 2005;Pimbblet 2005;Hao et al. 2011;Song & Lee 2012;Tempel et al. 2015). However, its influence should also be evident in the presence of accretion-driven shocks in the intracluster medium (ICM), especially at a large clustercentric radius (e.g. Tozzi, Scharf & Norman 2000;Ryu et al. 2003;Vazza et al. 2011). Such shocks arise naturally from the cosmological accretion of gas on to the cluster from the warm (T 10 4 -10 6 K) diffuse intergalactic medium (IGM) that resides in its lower density environs (see, e.g. Davé et al. 2001, and more recently Cui et al. 2018 andMartizzi et al. 2019). This gas will be accelerated gravitationally to peculiar velocities of the order of ∼ 1000 km s −1 during infall, which will produce strong shocks (e.g. Zinger et al. 2018) that give rise to various forms of non-thermal emission, ranging from high-energy X-rays and γ -rays (e.g. Keshet et al. 2003;Keshet, Waxman & Loeb 2004;Zandanel, Pfrommer & Prada 2014) to radio continuum emission in the form of synchrotron radiation as electrons spiral along magnetic fields (Hoeft et al. 2008;Vazza, Brunetti & Gheller 2009).
Both observations and simulations show trends between a cluster's merging and mass accretion history -which, observationally, are inferred from estimates of its dynamical state -and properties of its ICM, such as gas ellipticity (e.g. Chen et al. 2019), turbulent motions (e.g. Nagai et al. 2013;Avestruz, Nagai & Lau 2016;Lau et al. 2017;Shi, Nagai & Lau 2018), and the presence of density and pressure inhomogeneities (e.g. Nagai, Daisuke;& Erwin 2011;Zhuravleva et al. 2013;Roncarelli et al. 2013;Lau et al. 2015;Tchernin et al. 2016). This is of particular interest because these trends contribute to biases in cluster mass estimates (e.g. Zhuravleva et al. 2013); for example, gas turbulence and accelerations undermine assumptions of hydrostatic equilibrium, which can underestimate a cluster's true mass by between ∼5 and 20 per cent (e.g. Nelson et al. 2014), while inhomogeneities drive clumping factors that boost estimates of the X-ray emissivity (e.g. Planelles et al. 2017;Ghirardini et al. 2018), which can be up to 40 per cent in the cluster outskirts, leading to cluster masses being underestimated by up to ∼30 per cent (e.g. Rasia et al. 2014).
However, examining how substructure impacts properties of the ICM is of more general, astrophysical, interest. Examining the nature and origin of these trends between a cluster's assembly history, its dynamical state, and evidence of inhomogeneities in its ICM can help to clarify our understanding of how different physical processes can shape the ICM. For example: What are the relative contributions to turbulence of, for example, active galactic nucleus (AGN) jets and outflows compared to infalling substructures, and how does this change with position within the cluster and its dynamical age? Developing an understanding in this way can provide theoretical insights into potentially new and powerful observable signatures and tests of cluster assembly.
Previous theoretical studies have investigated how merging and accretion influence the thermodynamical structure of the ICM (e.g. Avestruz et al. 2016) using a statistical approach to infer the likely contribution of substructure (such as, e.g. using the density distribution in radial bins to decompose the ICM into smooth and clumpy components; cf. Zhuravleva et al. 2013), and comparing to the degree to which a cluster is considered dynamically relaxed (e.g. Zhuravleva et al. 2013;Rasia et al. 2014;Ota, Nagai & Lau 2018;Zinger et al. 2018) or by quantifying its assembly history by measuring mass growth over time (e.g. Nelson et al. 2014;Lau et al. 2015). In this study, we take a complementary approach; we look in detail at how the passage of substructures through the ICM, tracked from infall from the cosmic web, influences its thermodynamical properties, as inferred from radial profiles of density, temperature, and radial velocity. In particular, we focus our analysis on ICM properties in the cluster outskirts, between 0.3R 200 and ∼ (2-3)R 200 . It is here that we expect the signature of mergers and accretion to be most pronounced (see, e.g. Lau et al. 2015) -we reason that infalling substructures are likely to drive shocks into their surroundings when they are at their most gas-rich, subsequent to their accretion on to the cluster from the cosmic web -and we expect the effects of feedback to be secondary.
We base our analysis on a suite of cosmological zoom simulations of a single massive galaxy cluster, run with a range of state-ofthe-art astrophysical codes -spanning mesh-and particle-based hydrodynamics solvers, as well as galaxy formation prescriptions -as part of the NIFTY galaxy cluster comparison project (see Sembolini et al. 2016a, and subsequent papers). Previous papers in this series (e.g. Cui et al. 2016;Sembolini et al. 2016a,b) have demonstrated that there are significant differences between code predictions in the inner parts of clusters, reflecting both the choice of hydrodynamics solver and the adopted prescription for galaxy formation. This allows us to assess if ICM properties in the outskirts are similarly affected. If shocks driven by infalling substructures are driving turbulence and inhomogeneities in density and pressure, then we might expect systematic differences in the strength and duration of shocks because of the choice of hydrodynamics solver -as has been suggested by the work of Rasia et al. (2014) -but broad agreement in the timing and locations of shocks because substructure orbits will be governed by the large-scale gravitational field.
The structure of this paper is as follows. In Section 2, we briefly describe the simulations we have used, including how they were set up and their bulk properties at z= 0. In Section 3, we present our main results -connecting features evident in the radial profiles of ICM properties (gas density, temperature, and kinematics) with substructure in Section 3.1; exploring how the dynamics of substructure are imprinted on the ICM in Section 3.2; and demonstrating the influence of the cosmic web in Section 3.3. Finally, we summarize our results in Section 4.

The simulations
Our analysis focuses on cosmological hydrodynamical zoom simulations of a single galaxy cluster -'Cluster 19' from the MUSIC suite 1 (Sembolini et al. 2013Biffi et al. 2014) and identified originally in the MULTIDARK 2 2048 3 particle parent cosmological N-body simulation (Prada et al. 2012). The adopted cosmology assumes total matter, baryon, and dark energy density parameters ( m , b , )=(0.27, 0.0469, 0.73); a power-spectrum normalization of σ 8 = 0.82; a primordial spectral index of n = 0.95; and a dimensionless Hubble parameter of h = 0.7, all in accord with the WMAP7+BAO+SNI data set of Komatsu et al. (2011).
Initial conditions for all of the MUSIC clusters were generated using the approach of Klypin et al. (2001), which can be summarized as follows: (i) All particles within a sphere of radius 6 h −1 Mpc at z = 0 centred on the cluster in the parent MULTIDARK simulation are found in the low-resolution 256 3 particle version of the parent.
(ii) These particles are mapped back to the parent's initial conditions to identify the Lagrangian region from which they originated, from which a mask is created.
(iii) The initial conditions of the parent simulation are regenerated on a finer mesh of size 4096 3 , a factor of 8 increase in mass resolution relative to the parent simulation.
(iv) The mask is then applied to the 4096 3 particle data set, such that particles within the masked region are retained and those outside of the mask are binned to produce coarser mass resolution tidal particles, equivalent to the low-resolution 256 3 version of the parent.
The result is a set of initial conditions in which particles in the high-resolution patch have a mass resolution a factor of 8 higher than in the parent run, corresponding to dark matter and gas particle masses of m DM = 9.01 × 10 8 h −1 M and m gas = 1.9 × 10 8 h −1 M respectively.
We used a selection of astrophysical codes -AREPO, G2-X, G3-MUSIC, G3-OWLS, and G3-X -and their associated galaxy formation models (non-radiative and full physics) in this paper; in the cases of AREPO and G3-MUSIC, we included two variants of their galaxy formation models, with (AREPO-IL) and without AGN feedback -(AREPO-SH, G3-MUSIC-SH, and G3-MUSIC-PS). Key features are presented in the Appendix, and we refer the interested reader to Sembolini et al. (2016b) for more details. We use only codes that have used the same aligned parameters (see table 4 in Sembolini et al. 2016a) to resimulate the selected cluster; these govern the accuracy of the gravity solvers, and by requiring alignment, we ensure that cluster features that are driven by gravity -for example, large-scale filamentary structure, orientation, timing of merger and accretion events -should be consistent between runs. This allows us to focus on differences between results that are driven by either the hydrodynamics solver or the adopted galaxy formation prescription.

Structure finding
We have used the phase-space structure-finder VELOCIRAPTOR 3 (Elahi et al. 2019) to identify the main cluster and its substructure. VELOCIRAPTOR identifies particle groups using a threedimensional friends-of-friends (3D-FOF) algorithm and then refines each 3D-FOF group using an FOF algorithm applied to the full 6D phase space. This 6D-FOF group catalogue is cleaned to correct for sets of haloes that are combined artificially into a more massive system, arising from spurious bridges of particles between haloes that occur during the early stages of mergers; this is done using the velocity dispersions of the 3D-FOF groups. Each 6D-FOF halo is then decomposed into a smooth background and its substructures by applying the phase-space FOF algorithm recursively, locating sets of particles that are dynamically distinct from the background. This approach is capable of finding both subhaloes and streams associated with tidally disrupted systems (Elahi et al. 2013). Cluster 19 has been studied extensively as part of the NIFTY galaxy cluster comparison project (cf. Sembolini et al. 2016a, and subsequent papers), and its z = 0 properties in the fiducial G3-MUSIC full-physics run are given in Table 1. We follow the standard convention of defining the cluster mass as M 200 , the mass enclosed within a spherical region of radius R 200 encompassing an overdensity of 200 times the critical density ρ crit at the given redshift z, i.e.
Following previous studies (e.g. Thomas et al. 1998;Power, Knebe & Knollmann 2012), we use the centre-of-mass offset to quantify the cluster's dynamical state, where R cm is the centre of mass of the material within R 200 and R cen corresponds approximately to the centre of the densest substructure within R 200 , as calculated via the iterative method of Power et al. (2003). The measured value of R = 0.04 indicates that the cluster is relatively dynamically relaxed, which is consistent with a visual impression of the cluster (see below).

Visual impression
In Fig. 1, we show projected gas density (upper panels) and temperature maps (lower panel) for the G3-MUSIC run at z = 0. On the left-hand side, we show results for the fiducial non-radiative run, while on the right-hand side, we show one of its full-physics counterparts. Positions are plotted relative to R cen and normalized by R 200 , both evaluated for the particular run; we project from within a cube of side 7R 200 , or equivalently out to a radius of 3.25R 200 . The projected density maps were calculated using the PY-SPHVIEWER package (Benitez-Llambay 2015), which computes smoothing lengths on a per particle basis and so naturally adapts to regions of local density, allowing the fine structure to be discerned. For projected temperature maps, we use MATPLOTLIB's HEXBIN with a fine grid of 1536 2 . Because the simulations were run with aligned parameters, we find excellent agreement between the runs on intermediate-to-large scales; the structure of the filamentary network, as traced by the gas, in which the cluster is embedded is in very good agreement, and there is excellent correspondence between the positions of the more massive substructures (cf. Sembolini et al. 2016a). The differences we see within the core of the cluster have been studied previously (cf. Cui et al. 2016), and are driven by hydrodynamical shocks and the physics of galaxy formation (i.e. cooling and feedback).
In contrast, we see clear differences between the runs on small scales. In the projected gas density maps, we see more high-density small-scale structure in the full-physics run than in the galaxy formation run throughout the volume, as we would expect in the Here we show projected gas densities (upper panels) and temperatures (lower panels) in the fiducial non-radiative and full-physics G3-MUSIC runs at z = 0. The inner soliddashed/outer solid circles indicate 1/2/3 times R 200 . In plotting maps of projected gas densities with PY-SPHVIEWER (Benitez-Llambay 2015), smoothing lengths are estimated on a per particle basis before projection, allowing for the fine detail in gas density to be discerned; temperature maps are constructed using the MATPLOTLIB HEXBIN function.
presence of radiative cooling. However, in the projected temperature maps, we see more sharply defined features in the non-radiative run -for example, the ridge extending from the upper left-hand side to lower right-hand side in the cluster core, as well as several cool, dense knots with R 200 -than in the full-physics run. We shall return to these observations in the next section.

R E S U LT S
We now investigate the thermodynamical properties of the ICM in the outskirts of the cluster at z = 0. In particular, we focus on radial profiles of gas density, temperature, and radial velocity in the regime 0.3 ≤ R/R 200 3, and investigate in detail the relationship between structures evident in these profiles, which we interpret as arising from shocks, and the orbital motions of substructures (Section 3.2) and accretion from the cosmic web (Section 3.3). Because we have a range of astrophysical codes and galaxy formation prescriptions, we can estimate the degree of variation that arises. Our reference simulation is the non-radiative G3-MUSIC run, and residuals (when shown) are with respect to this run. All figures include data from both non-radiative and full-physics simulations, unless otherwise stated. For reference, we show as the heavy dot-dashed curve the dark matter density profile (as found in the dark-matter-only G3-MUSIC reference run) scaled by the cosmic baryon density, bar / m . In the lower panel, we show the residuals of each of these profiles with respect to our G3-MUSIC nonradiative reference run. Note that there are two sets of AREPO and G3-MUSIC full-physics runs -the AREPO-IL run with AGN (heavy-dashed curves with filled circles), the AREPO-SH run with stellar feedback but without AGN (light-dashed curves), and the G3-MUSIC-SH and G3-MUSIC-PS runs (heavyand light-dashed curves) with stellar feedback and without AGN.

Radial profiles of ICM gas properties
We begin our analysis by considering spherically averaged gas density and temperature profiles measured in the non-radiative (solid curves) and full-physics (long-dashed curves) simulations in our sample. Here we estimate the gas density at radius R measured with respect to R cen as the total gas mass (i.e. summation over gas particle masses) within a spherical shell, divided by the volume of the shell; gas temperature is deduced from the average, massweighted, specific internal energy of particles in the shell, multiplied by 3μm p /2k B , where μ is the mean molecular weight, m p is the proton mass, k B is the Boltzmann constant, and we have assumed an adiabatic index of γ = 5/3. We show these spherically averaged profiles in Figs 2 (density) and 3 (temperature). All radii are normalized to R 200 , the value in the fiducial G3-MUSIC non-radiative run; densities are expressed in units of ρ crit , the critical density at z = 0; and temperatures are expressed in units of kelvin. Lower panels show residuals with respect to the fiducial. Note that there are two curves for the AREPO and G3-MUSIC full-physics runs, corresponding to the two sets of galaxy formation models provided -AREPO-IL (heavy dashed curve with filled circles), which has been run with AGN feedback, while the AREPO-SH (light-dashed curve), G3-MUSIC-SH (heavydashed curve), and G3-MUSIC-PS (light-dashed curve) runs have been run without AGN feedback (see Appendix for more details). For reference, in Fig. 2, we overlay the density profile of the dark- . Spherically averaged gas temperature profiles at z = 0. In the upper panel, we show the spherically averaged gas density profiles measured in the non-radiative runs (solid curves) and full-physics runs (dashed curves). Radii are again normalized to R 200 . In the lower panel, we show the residuals of each of these profiles with respect to the G3-MUSIC non-radiative run. We follow the same convention as in Fig. 2 to distinguish between the different full-physics AREPO and G3-MUSIC runs.
matter-only version of the cluster (dot-dashed curve), scaled by the universal baryon fraction, bar / m . We find general agreement in the amplitudes and shapes of these spherically averaged profiles measured in the different sets of runs, with differences no greater than 0.1 dex (∼25 per cent) over the radial range we consider, as indicated by the residuals. The magnitude of the variations is substantially smaller than seen for the central parts of the same cluster (see figs 8 and 9 in Sembolini et al. 2016b). There are small enhancements in the spherically averaged gas density -for example, at ∼ 2R 200 and ∼ 0.7R 200 -and they are broadly consistent with the dark matter density profile scaled by the universal baryon fraction, which provides a good approximation to both the amplitude and the shape of the gas density profiles, which is consistent with the findings of earlier hydrodynamical simulations (e.g. Lewis et al. 2000;Pearce et al. 2000) and the assumptions of analytical modeling (e.g. Makino, Sasaki & Suto 1998;Komatsu & Seljak 2001).
Enhancements in temperature are more pronounced. Within R 200 , there is a noticeable bump at ∼0.5-0.6R 200 , although its precise location and shape vary between runs. We find that the nonradiative runs produce temperature inhomogeneities that are more sharply pronounced than in their full-physics counterparts, which is consistent with our observations in Section 2.3. Although we expect that cooling in the full-physics runs should exacerbate existing inhomogeneities in temperature, the net effect is to remove lower entropy, colder, denser, gas from the diffuse phase, and weaken rather than enhance temperature inhomogeneities. We also note that the variety of feedback processes that are active within a galaxy cluster -for example, AGN and cosmic rays -will tend to smooth out temperature inhomogeneities (e.g. Planelles et al. 2017).
An instructive alternative to spherically averaging profiles is to consider measures of the distribution of values within each spherical dist , upper curves) at z = 0, deduced from the distribution of gas densities within the spherical shells used in Fig. 2. As before, we show results for non-radiative and full-physics runs (solid and dashed curves respectively), and radii are normalized to R 200 . In the lower panel, for clarity, we focus on the residuals of each of the median profiles with respect to the G3-MUSIC non-radiative run. Note that we distinguish between the two sets of full-physics AREPO and G3-MUSIC runs as in Fig. 2. shell, as explored in Zhuravleva et al. (2013) and subsequent papers. In Fig. 4, we show how the median (lower curves) and 90th percentile values (upper curves) of gas density, ρ med dist and ρ 90 dist , vary with radius. Here the gas density is the intrinsic value (for a particle or cell) tracked in the simulation. We find that the median profile (ρ med dist ) and the corresponding spherically averaged profile in Fig. 2 are in very good agreement within R 200 , at the 5-10 per cent level.
The enhancement in ρ med dist apparent at R 2R 200 is consistent across codes and galaxy formation models, but its magnitude differs between non-radiative and full-physics runs by ∼30 per cent. We associate this with substructure infalling along the cosmic web, evident in the upper right-hand quadrants of the upper panels in Fig. 1. At R ≤ R 200 , pronounced enhancements are apparent in ρ 90 dist ; the locations of these peaks are broadly consistent between codes and galaxy formation models, albeit with a significant scatter in magnitude. We find the stronger enhancements in the non-radiative runs (for the reasons already discussed), but the larger scatter (by a factor of ∼10 in the enhancement at ∼ 0.6R 200 ) in the fullphysics runs. As noted by Zhuravleva et al. (2013) and, for example, Avestruz et al. (2016), these are naturally interpreted as signatures of substructure. Fig. 5 shows that the density enhancement in Fig. 4 at R 2R 200 has a corresponding temperature enhancement, which is most easily discernible in the median profile (T med dist , lower curves) and is most sharply defined in the non-radiative runs (again, for the reasons already discussed). The enhancement evident in the spherically averaged temperature profile at ∼0.5−0.6R 200 is also apparent in the T med dist and 90th percentile temperature profiles (T 90 dist , upper curves), while smaller features are evident at the locations of the density enhancements in Fig. 4. The sharpness of these features varies systematically between non-radiative runs, which are in very good  Fig. 4, we consider only the residuals of the median profiles with respect to the G3-MUSIC non-radiative run, and we distinguish between the two sets of full-physics AREPO and G3-MUSIC runs using the same convention. agreement with each other, and the full-physics runs, which tend to produce less sharply defined features than in the non-radiative runs.
We have used snapshots finely spaced over the redshift range z = 0.5 to 0 to verify that the temperature enhancements within R 200 evident in these radial profiles occur frequently during the cluster's assembly. This strengthens our assertion that they are transient features whose position moves to smaller radii over a dynamical time, consistent with the motions of substructures whose passage shocks the ICM. In contrast, the feature at ∼ (2-3)R 200 identified at z = 0 varies in amplitude but has remained at approximately ∼ (2-3)R 200 since z = 0.5, which we argue is associated with the gas inflow and substructure motion in the cosmic web (see, e.g. Burns, Skillman & O'Shea 2010).
We note that Avestruz et al. (2016) considered the mean massweighted temperature of the ICM at a larger radius (defined relative to R 500 , the radius enclosing a mean density of 500 times critical density) of an ensemble average of clusters and found that the presence of substructure tended to reduce it compared to the temperature when substructure was excluded, as a result of the cooler gas associated with substructures bringing down the average temperature. This is consistent with our results -we consider the temperature distribution of all gas within the shell, and so the 90th percentiles will include the contributions of shocked gas, some of which can be associated with the 90th percentiles in density. We have verified that excluding overdense gas leads to an increase in the mass-weighted average temperature within shells, albeit our result is noisier than that of Avestruz et al. (2016).
How might different physical processes modelled in the simulations impact radial density and temperature profiles? In Fig. 6, we investigate how these profiles are influenced by the physical processes modelled in the simulation, focusing on the 90th percentile values of the density and temperature distributions within spherical shells where we expect differences between the simulations to be Figure 6. Radial profiles of the 90th percentile gas density (top panels) and temperature (bottom panels), computed as in Figs 4 and 5. We distinguish between non-radiative runs (upper panel) and full-physics runs "Without AGN" and "With AGN" (middle and lower panels, respectively). As before, there are two full-physics G3-MUSIC runs without AGN -G3-MUSIC-SH and G3-MUSIC-PS -which we distinguish by heavy-and light-dashed curves, while there is one AREPO full-physics run with AGN and one without. most marked. Here we group the simulations into (1) non-radiative runs and full-physics runs (2) with and (3) without AGN. Note that we zoom in on the radial range 0.1 ≤ R/R 200 ≤ 1 to assess the extent to which the physics of galaxy formation -particularly AGN feedback -might affect the ICM outside of where the central cluster galaxy resides.
(1) We see a sharply peaked density enhancement at ∼0.6R 200 , evident in all but one of the runs, which is most pronounced in the non-radiative runs. There is excellent agreement between the classic smoothed particle hydrodynamics (SPH) codes (G3-MUSIC, G2-X, and G3-OWLS) -unsurprising, given that they are variants of GAD-GET2/3 (Springel 2005) with the same gravity solver and the treatment of SPH -while AREPO underpredicts the magnitude of the density enhancement in comparison, but, otherwise, tracks the predictions of the classic SPH codes. Similarly, there is excellent consistency between the temperature profiles predicted by the classic SPH codes, which show an enhancement at ∼0.5R 200 , whereas AREPO underpredicts the temperature across most of the radial range, while the enhancement is offset to ∼0.4R 200 . Note that the AREPO density and temperature profiles begin to peel away from their classic SPH counterparts at ∼0.1R 200 , as we expect because of the difference in hydrodynamics solvers (e.g. Sembolini et al. 2016a).
(2) The trends evident in the non-radiative runs are also apparent in the full-physics runs without AGN -AREPO-SH, G3-MUSIC-SH, and G3-MUSIC-PS. The amplitude of the density enhancement at ∼0.6R 200 is reduced relative to the non-radiative runs, and we see from G3-MUSIC that the assumed stellar feedback model has an impact.
(3) The amplitude of the density enhancement at ∼0.6R 200 shows a considerable variation between the different full-physics runs that implement AGN feedback, with the two classic SPH codes G2-X and G3-OWLS bracketing the range of variation -the enhancement is peaked in G2-X, whereas it is effectively smoothed away in G3-OWLS. A similar range of variation is evident in the temperature profiles, with the amplitude, sharpness, and, in particular, the location of the peak of the temperature enhancement at ∼0.5R 200 ; for example, the enhancement in the G3-OWLS run peaks at ∼0.35R 200 , compared to ∼0.45R 200 in the G2-X and G3-X runs. Interestingly, the classic SPH runs predict systematically higher temperatures across most of the radial range when compared to the modern SPH run with G3-X and AREPO.
These observations show that the physical processes modelled in a simulation can influence the ICM density and temperature radial profiles at larger radii, but the most significant differences are associated with the density and temperature enhancements associated with substructure. Where we see systematic differences -for example, in the temperature profiles in the 'With AGN' runs -these can be attributed to expected differences between classic SPH, on the one hand, and modern SPH and AREPO, on the other.
To isolate the particular effect of feedback on the ICM density and temperature radial profiles, in Fig. 7 we look at the profiles in the AREPO-SH and AREPO-IL runs, which are run without and with AGN feedback, and in G3-MUSIC-SH and G3-MUSIC-PS, which are run without AGN feedback but utilize two different stellar feedback models (cf. Piontek & Steinmetz 2011). We show logarithmic differences with respect to the nonradiative counterparts, and plot medians (light curves) and 90th percentiles (heavy curves) in the radial range 0.3 ≤ R/R 200 ≤ 1; models with AGN are indicated by dotted curves overlaid with filled circles. The light-dashed horizontal lines indicate ±25 per cent with respect to the non-radiative run. As we would expect, the logarithmic difference is dominated by the density enhancement at ∼0.6R 200 ; logarithmic differences in the temperature profiles are small ( 0.05 dex), while the corresponding differences in density are predominantly 0.1 dex, although these increase where we see enhancements that we associated with substructure.
In Fig. 8, we shift our focus to the radial velocity profile of the gas within the cluster, reasoning that features in the temperature profile are a response to gravitational and hydrodynamical influences, namely shocks. Here we show the variation of the spherically averaged radial velocity of gas with radius, with radii normalized to R 200 while radial velocities are in units of km s −1 ; v rad < 0 (> 0) implies   inflow (outflow). For a system in hydrodynamic equilibrium, we would expect v rad 0, but, instead, we see a pronounced dip to ∼ −500 km s −1 at ∼ 0.7R 200 before a sharp rise to ∼ 400 km s −1 peaking at ∼ 0.5R 200 . There is an appreciable variation between the runs, but the trend is systematic.
Interestingly, we see a feature in the dark matter radial velocity (heavy dot-dashed curve), similar to the feature evident in the gas radial velocity, albeit smaller in amplitude and displaced outwards in radius. At radii R r 200 , we see substantial differences between runs of ∼200 km s −1 . It is noteworthy that the unstructured mesh code, AREPO, and the modern SPH code, G3-X, predict systematically more negative radial velocities than the other codes, which are based on classic SPH, which are in broad agreement with one another. Strikingly, the dark matter radial velocity profile is more negative than the gas radial velocity profiles by ∼200−400 km s −1 , but there are common features across the various profiles (for example, the bump at ∼2R 200 ).
In Fig. 9, we provide an estimate for the strength of the shocks in the cluster gas by calculating the spherically averaged Mach number of gas as a function of cluster-centric distance. There are several approaches that we could adopt to estimate where and when shocks occur, which include tracking spatially localized jumps in temperature and/or velocity (e.g. Vazza et al. 2009), which are relatively straightforward in mesh-based codes where there are well-defined boundaries between cells, or tracking jumps in entropy or internal energy over time (Keshet et al. 2003), which is particularly straightforward in particle-based codes. However, we follow Mohapatra & Sharma (2019) and use the root-mean-square (rms) Mach number, where v 2 is the rms velocity computed for all gas elements (particles or cells) within a spherical shell and c s is the sound speed, which we compute directly from the internal energies.
As we would expect, the locations of features in the spherically averaged temperature and radial profiles correlate with features in the rms Mach number radial profile. The magnitude of M rms is larger in the moving-mesh AREPO runs -both the non-radiative and fullphysics variants -and in the modern SPH G3-X run, compared to the classic SPH runs, in the regime of stronger shocks, where M rms 2 and R 1.5R 200 , by between 10 per cent and 25 per cent. The trend is for the average rms Mach number to increase with cluster-centric distance, with values of M rms 1 within 0.5R 200 Vazza et al. (2009) note that 'internal shocks' within the virialized region of a cluster tend to have Mach numbers of M 3, with occasional spikes as a result of mergers, while 'external shocks' tend to be much stronger and larger (M 10), arising from the accretion on to filaments and sheets (see also Miniati et al. 2001).
We have demonstrated that enhancements in gas density and temperature profiles within R 200 , as well as features in gas and dark matter radial velocity profiles, plausibly correspond to substructures evident in projected gas density and temperature maps (Fig. 1). We now investigate how substructure gives rise to the features in detail.

The impact of substructure dynamics on the ICM
So far we have identified the imprint of substructure on the radial profiles of the ICM gas density, temperature, and radial velocity; we now shift our focus to the substructure population and investigate how its dynamical effects impact the ICM.
In Fig. 10, we show the projected spatial distribution of haloes of mass M 200 ≥ 10 12 h −1 M and within 1.25R 200 of the cluster's centre in the full-physics simulations. Symbol size is proportional to halo mass, while symbol colour indicates gas fraction. The solid outermost circle represents the cluster's radius R 200 , while the dashed circle indicates the position of the peak of the circular velocity profile. The shaded regions mark the outer and inner bounds of the features identified in the spherically averaged temperature and radial velocity profiles (cf . Figs 3 and 8).
The most striking feature of Fig. 10 is the presence of two distinct massive substructures whose radial distances coincide with features in the gas radial velocity profile. As a case study, we investigate these substructures in more detail. Their trajectories are consistent with the trends we see in the radial velocity profiles -the outer, gas-rich substructure is infalling with a velocity ∼−1700 km s −1 , whereas the inner, gas-poor substructure is moving outwards with a velocity ∼1500 km s −1 -although their velocities are much higher than the spherically averaged velocities at these radii, and supersonic. Their velocity vectors are relatively well aligned, with an angle between them of ∼20 • , and their direction of motion aligns with the orientation of the large-scale filaments within which the cluster is embedded (see Section 3.3 for further discussion).
These results are consistent with non-radiative and dark-matteronly versions of these simulations, 4 which we check by cross- Figure 10. Projected spatial distribution of haloes within 1.25R 200 of the cluster. We include only those haloes with masses ≥ 10 12 h −1 M . Symbol size and colour scale with mass and gas fraction, with lighter hues indicating gas paucity and dark hues indicating gas richness. The virial radius R 200 is represented by the solid circle, the radius at which the underlying cluster halo circular velocity profile reaches its maximum is represented by the dashed circle, while the shaded regions indicate the inner and outer boundaries of the temperature and radial velocity profile features. matching haloes identified in the full-physics run with those in the non-radiative and dark-matter-only runs, and confirm excellent agreement between the spatial and kinematic distributions. Differences that we note in the dark-matter-only run -that the counterpart to the outgoing halo is at a slightly larger cluster-centric distance than in the full-physics and non-radiative runs -can be readily understood as a consequence of the change in bound mass arising from the stripping of its gas content.
In Fig. 11, we make explicit how the infalling substructure, highlighted in Fig. 10, impacts the ICM. Recall that this infalling substructure is gas-rich [f g ∼ (0.7−0.9) b / m ], and it has a large radial velocity (∼−1700 km s −1 ) relative to the cluster; this velocity is significantly larger than the motions of cluster particles in the same region -for dark matter, v r ∼ −150 km s −1 with dispersions of ∼1000 km s −1 ; for gas v r ∼ −250 km s −1 , with dispersions of ∼600 km s −1 . We plot the projected gas distribution, weighted by temperature, in a region 2.1 h −1 Mpc wide and 1.1 h −1 Mpc thick, centred on and in the orbital plane of the substructure. Contours indicate the dark matter, gas, and stellar density associated with the substructure, while arrows show the direction to the cluster centre ( x C ) and the direction of motion of the substructure relative to the cluster centre ( v C ). Fig. 11 shows the infalling substructure driving a shock into the ICM -there is evidence for shocked gas in advance of the substructure and trailing gas in its wake in all of the runs. The temperature enhancement associated with the shock is more Figure 11. Infalling group shock in a subset of the full-physics runs; from the top right-hand side to bottom left-hand side, G3-MUSIC, G3-OWLS, G3-X, and AREPO. Here we zoom in on a region 2.1 h −1 Mpc wide and 1.1 h −1 Mpc thick, centred on and in the orbital plane of the galaxy group that generates the shock that gives rise to the feature in the temperature profile at 0.9R 200 . We show the mean logarithmic temperature distribution along the line of sight in this slab in greyscale, where each pixel is 20 kpc in size. We also show contour maps of the average dark matter density in green, the average gas density in orange, and the average stellar density in red. The galaxy group's direction of motion relative to the cluster's centre-of-mass is shown by a blue arrow, while the direction towards the cluster centre is shown by a magenta arrow. pronounced in the direction of the cluster centre, towards which the density gradient is higher, than in the direction of motion. 5 That a fast-moving, infalling, gas-rich substructure might generate a shock in the ICM is not surprising; but what is the influence of the gas-poor substructure moving outwards? At z = 0, it is moving with a radial velocity of ∼1500 km s −1 ; however, its position and velocity do not correspond to the inner feature seen in the temperature and radial velocity profiles -its radial velocity is too high, and it is at too great a cluster-centric radius. By tracking its orbit to earlier times, we can say that it had a mass on infall of ∼6 × 10 13 M and it was gas-rich. At its pericentric passage of ∼ 0.15R 200 , it deposited a gas mass of ∼8 × 10 12 M in ∼50 Myr as it moved at speeds of ∼2750 km s −1 , and led to its baryon fraction plummeting from 0.5 b / m to b / m . It is 5 Code-to-code variations are evident in the structure of the ICM -the classic SPH codes (G3-MUSIC,  show more small-scale variations than either the modern SPH codes (G3-X) or the moving-mesh code (AREPO), and the differences are particularly striking when one compares G3-MUSIC, which contains several cool, dense knots within this region (lower left-hand side), with G3-X, where such knots are absent. These small-scale variations are especially pronounced in the non-radiative runs, consistent with previous studies (see, e.g. fig. 12 in Power, Read & Hobbs 2014), but the impact of these variations on our conclusions is quantitative, not qualitative.
this stripping event and the gas associated with it that produces the shock at ∼(0.5−0.6)R 200 .
To illustrate these observations, in Fig. 12, we show the projected radial velocities and temperatures (upper and lower panels, respectively) of gas in the vicinity of the inner temperature enhancement (left-hand panels) and the gas stripped from the outward-moving substructure (right-hand panels). To show the generality of our results and the consistency between code behaviours, we consider two representative cases -G3-OWLS at the top and AREPO at the bottom. Inspection of the left-hand panels makes clear an arc of shock-heated gas, both trailing the outwardly moving substructure and colliding with infalling gas; the right-hand panels reveal that the shocked cluster gas is propelled outwards by the cold, dense, fast-moving gas stripped from the outwardly moving cluster. The fastest outwardly moving portion of this stripped gas is spatially coincident with the inner edge of the shock, pushing material in front of it, and it is this outwardly moving wake that shock heats when it encounters infalling gas.

The role of the cosmic web
A consistent feature of the results in Section 3.1 is the presence of an enhancement in gas density and temperature at ∼ (2-3)R 200 . In contrast to the transient enhancements with R 200 , this feature is Figure 12. Gas properties around the inner shock in the G3-OWLS and AREPO full-physics runs (upper and lower panels, respectively). We show, in projection, the mean radial velocity (upper panels) and temperature (lower panels) of gas in a shell around the shock feature (left-hand panels) and the gas stripped from the outward-going substructure (right-hand panels) in the AREPO full-physics run. We also show the position of the outward moving substructure and its direction of motion, with the point and arrow colour coded according to its radial velocity. The outer solid circle is at R 200 , while the inner/outer dashed circles are at 0.4/0.8R 200 and encompass the inner shock. For clarity, we do not include the gas interior to the shell in projection.
long-lived, and we have argued that it is likely associated with the gas accretion from the cosmic web. We have also demonstrated that the massive substructures whose passage through the ICM drives shocks have been accreted on to the cluster from a preferential direction, which we interpreted as a signature of the infall along filaments.
We use DISPERSE (cf. Sousbie 2011) to characterize the filamentary structure in the vicinity of the cluster. DISPERSE computes a discrete density field from a Delaunay tessellation of particles and identifies filaments as ridges connecting maxima through saddle points in the underlying density field. The result is a 'skeleton' of the density field, consisting of a fully connected network of segments, which can be smoothed (by pairwise averaging of neighbouring segments) to focus on a given spatial scale. 6 Here we compute the skeleton using both dark matter and gas in cubes of size 4 and 20 h −1 Mpc on a side, centred on the location of the maximum density of the cluster. Fig. 13 gives a visual impression of the cosmic web in which the cluster is embedded. Solid red lines show a 3 h −1 Mpc wide slice of the network of filaments with a persistence threshold c = 0.01 in a 20 h −1 Mpc cube smoothed on to a 128 3 mesh, consecutively trimmed of arcs with robustness smaller than the mean plus 1 standard deviation and smoothed over a 3 h −1 Mpc spatial scale. While these results are calculated for the gas, we note that we recover a very similar skeleton on large scales using the dark matter density. Dashed circles indicate the locations of the features in the radial velocity profiles, while the solid circles correspond to the gas-poor outwardly moving substructure (lower in x-y, left-hand side in y-z) and gas-rich infalling substructure, respectively.
The result is visually striking - Fig. 13 makes clear how the trajectories of both substructures can be traced to larger scale filaments, while also highlighting the connection between the enhancements in gas density and temperature, as well as the spherically averaged negative radial velocity, at ∼ (2-3)R 200 , to the merging of several smaller filaments along the direction of the most robust one. This explains why the location of this feature should be long-lived, albeit varying in magnitude, whereas the shocks associated with infalling substructures within R 200 should be short-lived and dynamic.

S U M M A RY
We have used a suite of cosmological hydrodynamical zoom simulations of a single galaxy cluster, run with a range of astrophysical codes and galaxy formation models as part of the NIFTY comparison project (see Sembolini et al. 2016a,b), to study how the thermodynamical properties of the ICM in a galaxy cluster's outskirts (0.3 ≤ R/R 200 3) at z = 0 are influenced by the accretion of substructure from the cosmic web. Using the radial profiles of gas density, temperature, and radial velocity, we highlighted the presence of enhancements in gas density and temperature, which we noted are associated with shocks (as measured by the Mach number) and which are commonly interpreted as signatures of the passage of substructure (see, e.g. Zhuravleva et al. 2013;Lau et al. 2017). We showed how such features within R 200 can be linked to the substructure evident in projected gas density and temperature maps, and can be connected to the passage of substructures through the ICM.
In particular, we looked in detail at two massive substructures -one gas-rich substructure infalling on a preferentially radial orbit with a relative speed of −1700 km s −1 , producing the clear signature of a shock in the direction of its motion, and the other 6 The skeleton can be smoothed on a given spatial scale, or it can be filtered on a persistence threshold. Persistence characterizes the robustness of the topology (i.e. number of holes, tunnels, threads, etc.) of the density field above a given excursion threshold, which can be interpreted as a measure of the significance of topological features, similar to the signal-to-noise ratio in observations. Note that persistence only characterizes the robustness of pairs of critical points (such as a maximum-saddle point height). We can, however, also derive a continuous robustness ratio that directly relates to the contrast of filaments (arcs between the critical points) relative to their background (Sousbie 2011). This allows us to distinguish between strong sustained filaments connecting high-contrast nodes (high persistence level, high robustness ratio) and fainter, possibly more transient, features (low persistence level, low robustness ratio). gas-poor substructure, moving outwards with a velocity of 1500 km s −1 , and gravitationally influencing the gas in its vicinity. By constructing the cosmic web in the vicinity of the cluster, we established that the cluster accretes from two dominant filaments, but that a network of filaments combine at ∼ (2-3)R 200 to give rise to density and temperature enhancements, associated with the infall of gas and dark matter along the cosmic web. Because the trajectories followed by the substructures that give rise to the shocks within R 200 are fixed by the filaments from which they were accreted, this suggests a way to reconstruct the infall time and direction of massive substructures from observational data (e.g. radio sychrotron measurements of termination shocks associated with the infalling galaxy substructures; Brown & Rudnick 2011).
We note general qualitative agreement between the selection of astrophysical codes and galaxy formation models (non-radiative and full physics) used in this paper in terms of their recovery of spherically averaged cluster gas profiles, at the level of ∼ 25 per cent or 0.1 dex between runs. Over the radial range 0.3 ≤ R R 200 , we find broad agreement between predictions for the locations of enhancements in gas density and temperature and when and where shocks occur, but there are differences between predictions for the sizes of enhancements and how strong shocks -as gauged by the rms Mach number, M rms , most notably between full-physics runs. This is as we would expect. Both gas and substructure accretion on to the cluster will be driven by the large-scale gravitational field, which dictates when and where, and we expect excellent agreement between the codes because of our requirement that they are aligned (cf. Sembolini et al. 2016a), whereas shock strength will be driven by a number of factors, including the resolution and choice of the hydro solver, and so we expect a greater variation between codes.
Although we have considered only a single galaxy cluster, albeit multiple realizations with different codes and galaxy formation models, our results imply that the dynamical effects of accreting substructure on ICM properties will need to be accounted for when making accurate predictions for next-generation surveys, such as, e.g. the Square Kilometre Array in the radio continuum, tracing synchrotron emission (e.g. Giovannini et al. 2015;Vazza et al. 2015), eROSITA (e.g. Merloni et al. 2012), X-Ray Imaging and Spectroscopy Mission (XRISM; Tashiro et al. 2018), and ATHENA (e.g. Nandra et al. 2013) in X-ray. Conversely, these surveys can open up a new view of the merger and accretion histories of clusters.
In summary, we have investigated in detail how the passage of massive substructures influences the ICM at a large cluster-centric radius and shown how it connects to features in radial ICM profiles, which have been previously quantified in a statistical sense (e.g. Zhuravleva et al. 2013;Nelson et al. 2014;Rasia et al. 2014;Lau et al. 2015;Ota et al. 2018;Zinger et al. 2018). Our case study of the two massive substructures shows that the manner in which such systems drive shocks into the ICM is not uniform; in one case, we have the classical picture of a gas-rich, infalling, substructure driving a shock, while in the other, it is stripped gas, which was once associated with a now gas-poor, outgoing, substructure, that is driving a shock. Finally, we have verified that the predictions of different astrophysical codes and galaxy formation models are in broad agreement, predicting spherically averaged radial properties of the ICM that agree at the 25 per cent (0.1 dex) level.