Turbulence and Vorticity in Galaxy Clusters Generated by Structure Formation

Turbulence is a key ingredient for the evolution of the intracluster medium, whose properties can be predicted with high resolution numerical simulations. We present initial results on the generation of solenoidal and compressive turbulence in the intracluster medium during the formation of a small-size cluster using highly resolved, non-radiative cosmological simulations, with a refined monitoring in time. In this first of a series of papers, we closely look at one simulated cluster whose formation was distinguished by a merger around $z \sim 0.3$. We separate laminar gas motions, turbulence and shocks with dedicated filtering strategies and distinguish the solenoidal and compressive components of the gas flows using Hodge-Helmholtz decomposition. Solenoidal turbulence dominates the dissipation of turbulent motions ($\sim 95\%$) in the central cluster volume at all epochs. The dissipation via compressive modes is found to be more important ($\sim 30 \%$ of the total) only at large radii ($\geq 0.5 ~r_{\rm vir}$) and close to merger events. We show that enstrophy (vorticity squared) is good proxy of solenoidal turbulence. All terms ruling the evolution of enstrophy (i.e. baroclinic, compressive, stretching and advective terms) are found to be significant, but in amounts that vary with time and location. Two important trends for the growth of enstrophy in our simulation are identified: first, enstrophy is continuously accreted into the cluster from the outside, and most of that accreted enstrophy is generated near the outer accretion shocks by baroclinic and compressive processes. Second, in the cluster interior vortex stretching is dominant, although the other terms also contribute substantially.

The resulting ICM driving motions on scales that range up to at least 100s of kpc will generally include weak-to-moderatelystrong shocks and hydrodynamic shear, both of which are expected to lead to turbulent motions that cascade downwards towards dissipation scales.
The solenoidal motions will stretch and fold structures, so are primarily responsible for amplifying and tangling the ICM magnetic field (e.g., Porter et al. 2015;Beresnyak & Miniati 2015). The compressive turbulence component will, itself, produce weak shocks that can, in turn, generate solenoidal motions (e.g., Porter et al. 2015). Both compressive and solenoidal turbulent components may accelerate cosmic rays through second-order Fermi processes (e.g., Fujita et al. 2003;Brunetti & Blasi 2005;Brunetti & Lazarian 2007, 2016. Several previous simulation efforts have measured the energy ratio between compressive and solenoidal motions in the ICM, finding a predominance of solenoidal motions (e.g., Ryu et al. 2008;Iapichino et al. 2011;Vazza et al. 2014, ). Interplay between the turbulence and shocks may be important in other respects, as well. For instance, turbulent amplification of magnetic fields by shocks and associated second-order Fermi acceleration leading to radio relic emission has been explored in several recent studies (Iapichino & Brüggen 2012;Donnert et al. 2016;Ji et al. 2016;Fujita et al. 2015;Fujita, Akamatsu, & Kimura 2016;Donnert et al. 2016). The relative contributions from solenoidal and compressive turbulent components will depend on the manner in which the turbulence is generated (Federrath et al. 2010;Porter et al. 2015) and its intensity (Vazquez-Semadeni 1994). In the ICM, each of these conditions is likely to vary significantly in both space and time.
The present work is motivated particular by the primary need to establish when, where, how and at what level the two turbulence components are produced and what is their relation to cluster formation dynamics. Here we focus on the turbulence itself, postponing its applications to subsequent works. We focus on turbulence generation, both solenoidal and compressive, and its connections to local ICM dynamical conditions. This complements previous simulation studies that have examined the global energetics of ICM turbulence evolved during cluster formation, including its association with major merger activity.
Particularly when issues such as magnetic field amplification and ICM dissipative processes, including cosmic-ray acceleration, are involved and when their dependences on local conditions are important (e.g. Subramanian et al. 2006;Brunetti & Jones 2014), it can be essential to separate solenoidal from compressive turbulent motions. For example, in recent work Miniati (2015) showed that the cluster-wide ICM compressive turbulence component is likely to have a steep (Burgers-law-like) spectrum, greatly reducing the power available for cosmic ray acceleration compared to a Kraichnan-like spectrum unless that power can cascade to very small scales, where it can more efficiently transfer energy to the cosmic rays.
In order to establish and evaluate the physical roles of turbulence it is essential to separate truly turbulent, uncorrelated flows from correlated, large scale bulk motions and shocks. Uncorrelated flows cascade energy and vorticity to small scales where they work to amplify magnetic fields and dissipate into heat and nonthermal particle energy. Coherent flows, on the other hand, carry signatures of global dynamical events, but are less directly connected to dissipation and magnetic field development.
Power spectra and structure functions constructed from simulation cluster-wide velocity fields typically suggest outer coherence scales ∼ 1 Mpc (e.g., Vazza et al. 2009b;Miniati 2014). While these scales correctly capture dominant, energy containing processes for the entire cluster, they do not, as emphasized above, necessarily discriminate against non-random, so, non-turbulent motions. They also span highly inhomogeneous, often stratified volumes whose motions on moderate to small scales are often too separated to be well connected causally when local driving conditions vary abruptly in response to nonspherical accretion or interactions (including mergers) with halos. So the ability of such global statistics to represent turbulent motions on the scales where they are most influential is limited. In that context, a more "local" approach seems better motivated. One strategy of this kind was suggested by Vazza et al. (2012) and Vazza et al. (2014). We will follow this strategy here in order to understand more clearly the generation, evolution and dissipation of the solenoidal and compressive turbulent motions produced during cluster formation. The following section outlines our simulations. Section 3 provides a summary of the several analysis tools we employ in this work, while Section 4 presents results of these analyses applied to a selected cluster simulation. Section 5 provides a brief summary and conclusion.

SIMULATIONS
We carried out, using the publicly available ENZO code (Bryan et al. 2014), multiple cosmological simulations designed to follow closely the formation of clusters selected to have a broad range of evolutionary histories.
The simulations applied the WMAP7 ΛCDM cosmology (Komatsu et al. 2011), with Ω0 = 1.0, ΩB = 0.0445, ΩDM = 0.2265, ΩΛ = 0.728, Hubble parameter h = 0.702, with a normalization for the primordial density power spectrum of σ8 = 0.8 and a primordial index of n = 0.961. All runs were non-radiative. No non-gravitational sources of heating were present except for an imposed temperature floor of T = 3 · 10 4 K in the redshift range 4 z 7, tuned to mimic the effects of reionization at low redshift.
We generated initial conditions (IC) separately at z = 30 for each individual simulation using 2 levels of nested volumes with comoving dimension, L0 = 44 Mpc/h = 62.7(≈ 63)Mpc. This technique is the same as introduced in Wise & Abel (2007).
First, low-resolution runs of several independent cosmological volumes were investigated in order to select the most massive objects in the volume. Second, new IC were generated by nesting two grids of 400 3 cells each and two levels of DM particles (400 3 each) with increasing mass resolution. The total volume was rotated in order to host the formation of the pre-selected cluster at the centre of the domain.
Inside the central (L0/10) 3 volume of each box (= 6.27 ≈ 6.3 Mpc) 3 (comoving), which is large enough to include the virial radii of most of our clusters, we further refined the grid by a factor 4. That increased our innermost spatial resolution to ∆x = 13.8/h kpc ≈ 20 kpc.
The generation of shocks and turbulence in simulated flows may be subject to spurious numerical effects, especially when adaptive mesh refinement is concerned (e.g., Miniati 2014;Schmidt et al. 2015). We wanted to avoid spurious effects caused by an uneven grid structure over time or space in the cluster formation region; the imposed fixed mesh refinement scheme puts 100 percent of the inner sub-volume uniformly sampled at 14/h kpc ≈ 20 kpc. The desired behavior was obtained in ENZO by "flagging" all cells within the volume of interest and using the AMR scheme to compute first the intermediate level of 40 kpc, and second the final level spanning the same sub-volume at 20 kpc. That procedure also ensured a conservative reconstruction of the fluxes across the coarse boundary of the 6.3 Mpc-sized, "zoom" region with the outer, lower-resolution volume, thus minimizing the noise in the refined reconstruction of infalling matter from the peripheral regions. This refinement procedure goes beyond previous sets of simulations by our group, where only a fraction (even if quite large) of the cluster volume was refined with tailored AMR (Vazza et al. 2009b).
The resulting full "Itasca simulated cluster" (ISC) 1 sample consists of 20 clusters with the above time and spatial information and has been designed to let us examine a variety of formation scenarios in detail. Each cluster run required about ∼ 12000 − 13000 cpu hours (about 1300 root-grid time steps and ∼ 10 5 top-level time steps in total). For analysis purposes we saved one data cube of hydrodynamical and N-body properties following sequences of 10 root grid time steps before z = 1, and then after every single root grid time step for the remaining z 1 evolution. This typically lead to ∼ 200 data dumps being retained for analysis. The dump time resolution after z = 1, while not constant, was generally ∼ 50 Myr.
In this paper we limit our analysis to one small-mass cluster, designated as "cluster it903", that underwent a major merger event ending around z ∼ 0.3. At the end of the simulation (z = 0), it903 had a total mass, Mtot ≈ 10 14 M ⊙ , with core and virial radii, rc ∼ 100 kpc, rvir ∼ 1 Mpc, respectively. The core temperature, Tc ∼ 2×10 7 K, corresponding to a sound speed, cs ≈ 660 km/sec. We defer the analysis of the full ISC cluster sample to future work. Figure 1 shows the evolution of the total enclosed gas mass and of the mean gas temperature for our "full" high resolution region (6.3 Mpc) 3 for cluster it903 and for an inner (1.44 Mpc) 3 , "clustercentred" region moving with that cluster. We use the inner, clustercentred volume extensively in the following sections for the study of turbulence and enstrophy since it reveals events involving the cluster more distinctly than the larger volume. The plot shows the relatively steady mass growth of it903 over cosmological times and the major merger event beginning after z ≈ 0.4 (followed by a decrease in the central gas mass, due to the outflow of the gas mass initially attached to a merger-involved subunit). This mass history is accompanied by a sharp peak in average gas temperature at z ≈ 0.35 in the cluster-centred, zoomed box and less dramatic, somewhat later (z ≈ 0.3) temperature spike in the full high resolution region 2 .

Visual impression of cluster it903
Three-dimensional volume rendering snapshots 3 of the (6.3 Mpc) 3 high resolution volume are shown in Fig. 2 at redshifts z = 1, z = 0.5 and z = 0.32, which outline the evolution of the cluster into the most significant merger events mentioned above. The left column of images highlights regions of high gas density, so provides a general sense of the mass merger history. The centre image column highlights the 3D shock distribution in this volume at the same times. The shocks are color coded by Mach number for 1 The "Itasca" label refers to the HPC cluster at the University of Minnesota used to compute the simulations. The website of the project is accessible at http://cosmosimfrazza.myfreesites.net/isc-project. 2 To avoid confusion later on we mention here that in some analyses we employ a larger, 5.76 Mpc 3 , cluster-centred box, since that size roughly matches the virial radius of the final cluster. For specificity we will refer to this as the "cluster virial volume". 3 These renderings assign color and opacity values to each voxel in a volume depending on the value of a rendered quantity (e.g., density), then construct perspective views using volume ray casting (http://www.lcse.umn.edu/hvr/hvr.html).

Figure 1.
Top: evolution of the integrated gas mass (black) and of the gas mass increment/decrement (red/blue) for cluster it903 from z=1 to z=0. Bottom: evolution of the volume-weighted mean temperature for it903. The solid lines give the evolution for the"innermost", cluster-centred (1.44 Mpc) 3 volume used for our primary turbulence analysis, while the dotted lines refer to the full (6.3 Mpc) 3 peak resolution (∆x ≈ 20 kpc) volume.

1.5
M 20 (red-yellow-white). The right column then displays the 3D distribution of enstrophy (ǫ =(1/2) vorticity 2 ) as an easy-to-compute and very useful proxy for the turbulence velocity distribution. It is clear from these images, as we discuss in detail below, that the enstrophy distribution is well-connected to the shock distribution, albeit very different in detail. Fig. 3 presents 2D slices at z = 0.32 through the centre of the same volume and along the same line of sight as Fig. 2. The left image shows the gas density, the middle image the gas temperature, and the right image the shock distribution, again, color coded by Mach number. M 20 red-yellow-white). Right panels render the enstrophy distribution (arbitrary units).

Shock finder
We detect shock waves in post-processing analysis using the algorithm presented in Vazza et al. (2009). The scheme is based on an analysis of 1D velocity jumps across cells. The minimum of the 3D divergence of the velocity, ∇ · v, is used to identify the centre of the shock region (see also Ryu et al. 2003;Skillman et al. 2008). Typically, shock transitions span about 2-3 cells along the shock normal. The one-dimensional Mach number for flagged transition is constructed from the 1D velocity jumps along each scan axis using the Rankine-Hugoniot relations. The final Mach number is constructed from a combination of the three 1D solutions. Shock surfaces are then approximated as the ensemble of the face areas of cells tagged as shocked by the scheme. This method has been extensively tested against similar methods used in grid codes (e.g., Vazza et al. 2011b), and has proven to be an efficient and accurate measure of shock waves in cosmological runs. The kinetic power across each flagged cell shock surface, which provides a useful metric for energy available to dissipation in shocks, is given by where ρu is the co-moving up-stream density, v sh = Mcs is the co-moving speed of the shock, M is the inflowing Mach number and ∆x is the cell size. Total kinetic energy flux through shocks is then a sum from Eq. 1 over flagged shock cells. However, we point out that the identification and characterization of shocks, and especially shocks with M 1.5 and oblique on the grid is made uncertain by the relatively larger numerical errors associated with very small jumps in velocity, by numerical smearing of the shock transition (Skillman et al. 2008) as well as by the presence of significant temperature and velocity fluctuations in the ICM, which add uncertainties to estimates of pre-shock values (Vazza et al. 2009). In order to bracket the role of, mostly inconsequential, weak shocks in the following turbulent analysis, we will there, also, present complementary results obtained by masking out regions close to identified shocks with M thr 1.2. We apply this lower M thr bound to simplify procedures.

Filtering of turbulent motions
The extraction of turbulent motions within the cluster 3D velocity field requires a proper filtering of often comparable coherent velocity components (characteristically larger-scale) from uncorrelated, turbulent velocity components that cascade to small scales. As noted above, the roles of solenoidal (rotational, incompressive) and compressive (irrotational) turbulent motions (defined respectively by ∇ · v sol = 0 and ∇ × vcomp = 0) each have important roles in the ICM. So, as an additional step we also separated the velocity field (both filtered and unfiltered) into those elements. To accomplish these objectives we combined several steps that we previously proposed and tested individually in Vazza et al. (2012) and Vazza et al. (2014).
As an initial step to reduce numerical noise and finite difference cross-talk between divergence and curl operations, we applied a first order velocity smoothing filter to the initial full velocity field of the simulation. This has a benign influence on extracted turbulent velocity fields 4 .
For our primary turbulence filter we applied the iterative, multi-scale velocity filtering techniques based on Vazza et al. (2012) to the (6.3 Mpc) 3 subvolume of each ENZO snapshot. This filter reconstructs the local mean velocity field around each position, r, by iteratively computing the local mean velocity field in the "n th " iteration as where the sum is over cells within a domain of radius, Ln, and where wi is a weighting function. In this work we simply set wi = 1 and use a volume-weighting, while in other applications in more stratified media wi = ρ (i.e. density-weighting) is a more appropriated choice. However, given the rather small filtering scales reconstructed by our algorithm in the innermost cluster volume considered in this work (Sec.4), the differences between the wi = 1 and the wi = ρ are very small, as discussed below. The local small-scale, fluctuating velocity field within the radius, Ln( r), relative to position r, is then computed as δv(Ln( r)) = v − v(Ln) for increasing values of Ln. Iterations are continued until the change in δ v between two iterations in Ln falls below a given tolerance parameter, which, based on our tests, we set to 10%. The resulting | δv(Ln)| provides our best estimate for the turbulent velocity magnitude for an eddy-size L eddy ≈ 2 · Ln.
We observe that, while in Vazza et al. (2012) we used the local skewness of the velocity field as a fast proxy to tag shocks, in the present work we can access this information in a more accurate way through the (obviously more computationally intensive) shock finding procedure outlined above (Sec.3.1). Therefore, we excluded shocks by simply stopping the iterations whenever a shocked cell entered the domain. That is, the length, Ln then represents the distance to the nearest "influential" shock. On the other hand, our procedure is not designed to explicitly filter out the contribution from velocity shears, e.g. at the contact discontinuity generated by sloshing cold fronts (e.g. Zuhone & Roediger 2016). While in principle the presence of such discontinuities might introduce a small spurious contribution to our measured turbulent budget, this spurious signal is small compared to the turbulence induced by mergers (e.g. Vazza et al. 2012). In particular, the cluster studied in this first paper is a highly perturbed one, where the formation of sloshing cold fronts is highly unlikely (e.g. Zuhone & Roediger 2016).
Our results here, as well as previous cluster simulations are roughly consistent with the behaviour of solenoidal turbulence fol-lowing the classic, Kolmogorov picture in which |δ v| ∝ L 1/3 n ((e.g., Ryu et al. 2008;Xu et al. 2011;Vazza et al. 2011a;Miniati 2014). Consequently, while the r.m.s turbulent velocity or the total turbulent pressure depend on Ln, the solenoidal kinetic energy cascade flux, defined as: is insensitive to the specific value of Ln as long as it is measured within the inertial range of the turbulence. Obviously, departures from this behavior can appear if the turbulent behaviour is very different from Kolmogorov and other, distinctive flow patterns are important (e.g., coherent shock waves, see Sec.3.3.1). In the ideal case, our procedure constrains f KE,turb at the outer scale of turbulence, and L turb = 2L, as shown in the Appendix (Sec. A1.1). However, in practice the iterations are stopped before reaching this scale, and f KE,turb is computed within the turbulent cascade. For this reason, in the following we will regard L as a filtering rather than a turbulent scale. Reconstruction of the second is difficult for multiple reasons, including nonequilibrium and highly inhomogeneous flows on large scales. Therefore, in general L turb 2L. As a comparison to the iterative filter we also present below some results in which a simple, fixed scale (L f ∼ 1 Mpc) was used.
As an additional test, we have verified that the usage of a density-weighting in Eq.2 leaves our results basically unchanged. In particular, the kinetic energy flux measured by Eq.3 is increased only by a factor ∼ 2 at most, when using the density-weighting within the central (1.44 Mpc) 3 volume studied in the following (Sec.4).

Solenoidal & Compressive Motion Decomposition
In order to characterize the dynamical properties of cluster turbulence we decomposed, both, the filtered and un-filtered 3D velocity fields of the simulations into solenoidal and compressive elements using the Hodge-Helmholtz projection in Fourier space (e.g., Kritsuk et al. 2011). Here we outline our fiducial method to carry out the decomposition, while in the Appendix (Sec. A1.2) we compare alternative approaches and control tests on our procedure.
Our fiducial decomposition algorithm first constructed the Fourier space velocity vector field, The compressive component in Fourier space, k ×˜ V comp( k) = 0, was found as the residual, , then produced the associated physical solenoidal and compressive velocity distributions, v sol ( r) and vcomp( r), where, again, r represents a point in the spatial domain. This procedure was performed both on the full, primitive 3D velocity data, yielding v sol ( r) and vcomp( r) and on the small-scale filtered field, yielding δv sol ( r) and δvcomp( r). One of our tasks in the present analysis is to compare properties of the two solution sets.
It is useful to point out here that the products of this analysis offer a useful way to estimate the local turbulent energy flux across scales. Specifically, for uncorrelated velocities, δvL, filtered on scale, L( r), the turbulent energy flux per unit volume was estimated using Eq. 3, for either the compressive, δvL,comp, or solenoidal component, δv L,sol . In steady, Kolmogorov turbulence (δvL ∝ L 1/3 ) these fluxes would be scale-independent, so, would provide robust estimates for the local turbulent energy dissipation rate per unit volume, ρ ǫ d . In the section 3.4 we outline an alternate, complementary approach to estimation of the solenoidal kinetic energy dissipation not requiring the above turbulence scale filtering.

The influence of shocks in Turbulent Energy Flux Budgets
The presence of cluster formation shocks is problematic to the turbulence component analysis. First of all, the numerically smoothed profiles of shocks contaminate to some degree the solenoidal Fourier field, V i,sol ( k) for larger k and thus v sol and δv sol on small scales. Fortunately, this issue is significantly mitigated by the velocity smoothing mentioned above, as discussed in Porter et al. (2015), for example. More importantly, while some shocks contribute appropriately to the compressive velocity element, vcomp, and sometimes to δvcomp, not all shocks, and in particular, structure formation shocks driven by coherent flows, are not elements of uncorrelated, compressive turbulence, δvcomp. The difficult issue, then, is one of judging which shock compressions to count as part of the compressive turbulent motions, δvcomp. In practice, some weak shocks are integral to the turbulence, including those generated by colliding turbulent motions (even solenoidal motions) (e.g., Porter et al. 2015), while other shocks, especially stronger ones and those whose extents exceed the cluster core scales, are more properly associated with the generation of (random) turbulence (e.g., Federrath et al. 2010;Porter et al. 2015), but are not elements of the turbulence per se. We are not aware of any simple, clean and robust way to establish this dichotomy. To explore the significance of this complication we carried out a series of numerical experiments in which we masked out patches of cells around shocked cells flagged using methods outlined in Sec. 3.1. That specific algorithm is identified below as shock limiting, since, as mentioned above, the length Ln is then limited by the separation scale of shocks with M thr . The resulting kinetic energy statistics were then compared to those obtained without "shock limiting". In each case we also examined the velocity fields as extracted from our interative turbulence filter (section 3.2) and for turbulence motions obtained using a fixedscale filter. To be conservative with respect to the numerical smearing of shocks, when we applied the masking procedure we removed shock centres, as well as the adjacent ±2 cells along the shock normal to ensure that the numerical shock profile ( 3 cells) is fully contained by the masking region.
While this procedure obviously excludes some kinetic energy on scales of a few cells as well as larger scale flows, we found that the net turbulent, kinetic energy fluxes were not sharply reduced by the shock limiting algorithm. As illustrated in the bottom panel of Fig. 4, the dominant influence of the shock limiter appears to be the reduction in the filtering length, Ln. Specifically, the shock-limited distribution for Ln is offset to smaller Ln by roughly a factor 3 from the non-shock-limited Ln distribution. On the other hand, the top panel, which presents radial profiles of the kinetic energy flux, FKE (Eq. 3) demonstrates that the difference between energy fluxes with or without shock limiting is generally less than ∼ 20 − 40%. The application of a cluster-scaled, fixed length turbulence filter, L f = 1 Mpc (comparable to the rvir of the system), however, led Solid lines show results with shocks masked out, while for dotted lines no masking was used. The dashed lines show the results obtained using a fixed L = 1 Mpc filtering scale without shock masking. Bottom: fractional number of cells with a given iterated filtering scales, 2 L, for the it903 cluster with (dotted curve) or without (solid curve) masking of shocks. The contribution of each cell has been weighted by its turbulent kinetic energy.
to seriously reduced energy fluxes, typically by factors ∼ 5 − 10. The fact that the fixed filtering scale L f = 1 Mpc systematically underestimated the kinetic energy flux reconstructed in the other approaches suggests that this scale is already larger than the true outer injection scale of turbulence in the domain. This follows from the Kolmogorov picture of turbulent cascades in the ICM, because on scales larger than the injection scale the kinetic energy flux (Eq. 3) is not conserved. Therefore, the values of r.m.s. velocities measured on these large scales are not truly representative of ICM turbulence.

Enstrophy as a Metric for Solenoidal Turbulence
Previous studies (Ryu et al. 2008;Vazza et al. 2014;Miniati 2014;Schmidt et al. 2014) and our results below suggest that ICM tur-bulent motions are predominantly solenoidal in character. The distinguishing property of solenoidal turbulence is, of course, that the motions are rotational; that is, they have non-vanishing local vorticity, ω = ∇ × v = 0. The vector vorticity, ω tends to average towards zero, so the vorticity magnitude is more a more useful tool. It is useful in this regard to recall that the eddy turn-over rate on a scale ℓ, 1/τ eddy ∼ δv ℓ,sol /ℓ ∼ |ω ℓ |, where the subscript on |ω| identifies this as representing circulation on the specific scale, ℓ. That is, vorticity is a measure of the rate at which eddies turn over. We will use this concept below to normalize vorticity measures in convenient units; i.e.,ω = ω · τ0, where τ0 is a representative timescale. Then,ω represents a characteristic number of eddy turnovers in the chosen interval, τ0.
As an additional perspective, we note that the square of the vorticity, or more directly the enstrophy, ǫ = (1/2)ω 2 , can be related in a turbulent flow to the kinetic energy content per unit mass, (1/2) v 2 sol , and dissipation of the solenoidal turbulence. This measure can then be matched to the solenoidal turbulent energy extracted through the filtering algorithm discussed in Sections 3.2 and 3.3. But, since no such filtering is involved in finding the enstrophy, the methods are complementary.
Formally, in terms of the turbulence one-dimensional velocity power spectrum, Es(k), the mass weighted solenoidal kinetic energy can also be written as: wherek 2 is the spectral-weighted mean of k 2 . The enstrophy can be obtained directly from the simulation data by application of the numerical, finite difference, curl operation on the primitive flow fields. A potential issue is that finite difference gradient operations on a compressible flow can pick up unphysical, numerical noise that obscures the signals of interest. Previously we referred to this as finite difference "cross-talk". However, Porter et al. (2015) demonstrated that these effects can be significantly ameliorated by using a simple smoothing operation on the velocity fields ("Favre filtering"), without significantly reducing the desired signal. Consequently, we employ the same approach in our enstrophy analysis here, employing a simple 3 3 cell-average smoothing. For Kolmogorov solenoidal turbulence the power spectrum can be written as Es(k) = C0η 2/3 d k −5/3 (e.g., Gotoh et al. 2002), where η d is the rate of solenoidal turbulent energy dissipation per unit mass, and C0 ∼ 3/2 − 2, is the so-called Kolmogorov constant. Applying the Kolmogorov form over a range of wavenumbers [k0, k1 = a k k0], with a k > 1 Eq. 5 leads to the relationships, The information in Eq. 6 can also be used to express the turbulent energy dissipation rate, or energy flux rate, in terms of either the solenoidal velocity or the enstrophy. In the limit a k ≫ 1 these become, where in the final forms we have set C0 = 1.8, v 2 sol ) 3/2 ∼ δv 3 L with k0 = 2π/L to match the energy flux relation in equation 3 and set a k = L/ℓ1. We tested the relations Eq. 7 and Eq. 8 using steady, driven, homogeneous turbulence simulation data (of known dissipation rate) presented in Porter et al. (2015) and found good agreement to within ∼ 15 % for both the velocity-based and enstrophy-based predictions. The turbulent energy dissipation rate per unit volume can be expressed as f turb = ρ · η d . We will apply these relations to our cluster simulation data in Sec. 4.2. In developed hydrodynamical turbulence the total rate of that dissipation is independent of the microphysical details, although it obviously provides upper bounds to energy input rates. One of the keys to a useful understanding of ICM turbulence is an understanding of how, where and when it is generated. Enstrophy tracking provides an effective and practical tool to study this for the dominant, solenoidal component. To this end we can derive an equation for enstrophy evolution using the curl of the Navier-Stokes equation (e.g., Porter et al. 2015); namely, where enstrophy source, sink and flux terms, F , (all called "source terms" below) are defined as, The enstrophy advection term, F adv , in Eq. 10 is conservative, so that its integral over a closed system must vanish. We will see over cluster volumes, however, that the integral of this term does not vanish. The F stretch , F baroc and Fcomp terms account for vortex stretching, enstrophy production in baroclinic flows and in compressions, respectively. Note that the fluid compression rate, −∇ · v, enters into both the F adv and Fcomp terms. However, whereas F adv always integrates to zero in a closed system, Fcomp does not if there is a net alignment of the velocity with the enstrophy gradient field. This alignment is usually present in shocks, so that Fcomp dV = v · ∇ǫ dV > 0 there, but is mostly small elsewhere in the absence of systematic compression (Porter et al. 2015). In driven turbulence "in a box" simulations this term was found to provide a good, overall measure of enstrophy growth by way of shocks (Porter et al. 2015). During cluster mergers there will be systematic compressions and rarefactions, and the behavior of the Fcomp term will reflect that, as well. True vorticity source terms such as vorticity creation in curved or intersecting shocks, or F baroc 6 , for example, while necessary to seed enstrophy in an irrotational flow, are mostly sub-dominant in homogeneous turbulence simulations once any vorticity exists in the flow. We will examine the varied roles of each of these terms in our cluster simulation data in Sec. 4.2.
For completeness we include in Eq. 9 and Eq. 10 the explicit viscous dissipation term, F diss , where G = (1/ρ)∇ρ· S, with S the traceless strain tensor (e.g. Porter et al. 2015). However, our simulations are based on Euler-limit hydrodynamics, where there is no 6 Curved or intersecting shocks can create vorticity even in isothermal flows, whereas F baroc cannot. So, these sources represent distinct physics. explicit viscosity, ν. Thus it is not possible to evaluate F diss explicitly in our simulations. On the other hand, effective net turbulence dissipation rates can be estimated using the turbulence relations in Eq. 7 and Eq. 8, which we do in Sec. 4.2. See Zhu et al. (2010); Schmidt et al. (2014) for previous analogous turbulence dissipation analyses in clusters. We will examine this issue more broadly for the ISC simulations in a subsequent paper.

RESULTS
As mentioned in Sec. 2, we focus this paper on the one cluster designated it903, while we defer the study of the complete ISC sample to future work.

Preliminary turbulence analysis.
We start our analysis of turbulence in it903 by studying the spatial distribution of the gas velocity field, filtered according to the methods presented in Sec. 3.2-3.3. Fig. 5 illustrates a 2D slice of the velocity field at z = 0.32, processed in several ways to reveal its turbulence properties: unfiltered total velocity along with its compressive and solenoidal components (top row) or small-scale filtered velocity field (as in Sec. 3.2) and its components (bottom row). The distribution of filtering scales used to remove large-scale motions has been shown in the bottom panel of Fig. 4.
The general evolution of this cluster, until the last major merger close to z ≈ 0.3, can be seen in the sequences of 3D volume-rendered images in Fig. 2. The merger developed along the upper-left, lower-right diagonal of this view. While the volume distribution of shocks visible in Fig. 2 is quite complex, the 2D slice in Fig. 3  Before entering through accretion shocks, the accreted gas reaches typical infall velocities of ∼ 500−700 km/s at this epoch. This flow is predominantly compressive, yet significant solenoidal velocity components are present even before crossing accretion shocks, where filaments break the spherical geometry of accretion shocks (e.g., in the top right corner of the image). Shock interactions also significantly enhance the solenoidal motions (see Sec. 4.3).
Inside the cluster the coherence of infall motions gets broken by irregular, converging flows and resulting shear motions. The maxima in the velocity field, associated with the density peaks of substructures, can reach ∼ 500 − 800 km/s. Even the unfiltered velocity field gives a clear visual impression of a predominance of solenoidal motions within the cluster. The principal exceptions are associated with strong shock waves sweeping through the volume at all times, but especially during the major merger.
The actual predominance by solenoidal turbulent motions is clearly revealed when large-scale laminar motions and shocks are filtered out (lower right panel of Fig.5). The total filtered velocity (lower left panel) is dominated by solenoidal motions everywhere except near shocks, where small-scale velocity structures (∼ 100 km/s) are impossible to distinguish between real smallscale turbulence and simpler shock jumps. Everywhere at distances 100 kpc away from shocks the small-scale velocities are almost entirely solenoidal. The absolute maxima of the small-scale solenoidal velocity field are found at the interface of filamentary accretions within the cluster volume, and also downstream of shocks. We will analyse the generation/amplification of vorticity at shocks in the next section.
We next analyse the turbulent nature of the flow using power spectra of the (unfiltered) velocity field shown in Fig. 6. For that we computed the 3-D power spectra within increasing volumes around the centre of it903 at z = 0.32, assuming, for that exercise, periodic boundaries in application of FFTs. The spectra for the total velocity field show the typical power-law behaviour of clusters simulated in this way (e.g., Vazza et al. 2011a), up to nearly two orders of magnitude in scale when the cluster virial volume(5.76 Mpc) 3 is considered. The spectrum flattens at low k-values, but the powerlaw is close to E(k) ∝ k −2 (i.e. steeper than Kolmogorov turbulence) for most scales. Fig. 7 shows the complementary view of the second-order longitudinal structure function for the same volumes, obtained by randomly extracting ≈ 5 · 10 4 paris of cells in the domain. The trend of the structure functions is similar to the results of Miniati (2014), with hints of a flattening at ∼ Mpc scales and of a steeper behaviour of the compressive component at small scales. In both panels we also show in colors the spectra/structure functions of the solenoidal and compressive components. Again, this shows how the solenoidal component is larger at most scales. However, in the largest box the difference between the two modes is reduced, and in this case the smallest scales are dominated by the compressive component, suggesting the relevant contribution of shocks forming in cluster outskirts. When larger volume are included, the compressive structure functions steepen at small scales, strengthening the view that shocked regions become increasingly more relevant on large scales. In both views the recovered trends are consistent overall with the picture of a turbulent ICM, mixed with large-scale regular velocity components for scales 0.1 − 1 Mpc and punctuated by small-scale velocity perturbations due to shocks.
We notice that our analysis detects a significantly steeper slope (both in the power spectra and in the structure function) in the solenoidal component, compared to the compressive component. This is at variance with some other recent numerical studies of turbulence in the ICM (e.g. Porter et al. 2015;Miniati 2015). However, understanding the origin of this difference is not trivial. Most of the difference is seen at small scales, when the impact of numerical dissipation is larger (e.g. Kritsuk et al. 2011). Moreover, the mass/dynamical state of it903 is different from the one analysed  in Miniati (2015), and also our method for the mode decomposition of turbulent modes is different (Sec. 3.3, see also Appendix A.1.2). Constraining the slope of turbulent modes at small scales in the ICM is relevant to estimate of (re)acceleration of radio emitting particles (e.g. Brunetti & Jones 2014;Miniati 2015;Brunetti 2016), but we defer a more extensive exploration of this issue to future work with our ISC sample.
In Fig. 8 we show the evolution of the power spectrum in the inner, cluster-centred (1.44 Mpc) 3 region of the two components at different epochs (top panel), and the ratio between the compressive and the total velocity power spectrum for the same epochs (lower panel). While the solenoidal component only shows variations within a factor 2 for k 2 at all epochs, the compressive component varies more significantly over time and at all scales. Consequently, the ratio between the compressive and the total power spectrum also shows significant variations: in the investigated interval it ranges from ∼ 5% to ∼ 50% at the largest scales, and from ∼ 30% to ∼ 50% at the smallest scales. However, it is worth noticing that that the actual difference between solenoidal and compressive turbulence is not properly captured by this simple ratio. In reality large-scale velocity fields introduce regular components at all spatial scale, which are best removed only by our multi-scale filter. Likewise, shocks very significantly bias the estimate of the real compressive turbulence at the smallest scales. It is not simple to evaluate the proper role of such shocks. Some, especially weak shocks with large curvature, are truly components of the compressive turbulence (i.e., uncorrelated flows), while stronger shocks with small curvature are not. In either case said shocks can become sources of turbulence and our approach is to try to allow for a range of possibilities. The average (volume-weighted) radial velocity profiles for different epochs of it903 are given in Fig. 9, and can be compared with the small-scale filtered profiles (bottom panel). We also show the average volume-weighted profiles of the compressive and solenoidal components. In the early stages of cluster formation the total velocity was large in the centre, ∼ 500 km/s. At later stages, after the major merger, it flattened at all radii and decreased to ∼ 200 km/s in the central regions, slightly increasing outwards. The unfiltered compressive velocity field is found to be larger than the solenodial field only in the centre of it903 at high redshift (z = 0.84), following supersonic bulk motions associated with fast infalling gas substructures (Fig. 2, top panel), while it is always smaller later on. The measured small-scale velocities (referred to within a ∼ 200 kpc scale in the shock limiter case, or ∼ 400 kpc in the case without masking of shocks, as in Fig. 4) using our filtering approach are of the order of ∼ 100 − 200 km/s at most epochs, with a very flat profile outside the cluster core. The small-scale compressive velocity component is found to be significant (but still smaller than the solenoidal one) close to the major merger event at z = 0.32. In particular, outside rvir we measure a very significant jump in the compressive small-scale velocity, ×2 − 3 larger than the increase in the solenoidal component at the same radius.
An important point to stress here is that neither of these two velocity profiles characterizes the full turbulent (uncorrelated) velocity field of the ICM. For the first, unfiltered case the contribution from laminar infall motions (clearly visible in Fig.5) biases the velocities high compared to random components, while in the second case our filtering procedure computes the r.m.s. random velocities only within the Ln reached before the filtering algorithm stops when finding a shock or converges on an average velocity within the scale, Ln. The scales Ln will generally underestimate the true outer scale for uncorrelated motions. Moreover, the inhomogeneity of the cluster limits the meaning of these scales, as r.m.s. velocity extracted using different Ln are not easily compared.
As we already commented in Sec. 3.2, the kinetic energy flux represents a more robust tool than velocity magnitudes to measure the consequences of random velocity components, because it is a relatively scale-independent measure across the turbulent Figure 11. Radial profile of the ratio between the kinetic energy flux of the small-scale compressive components and the total small-scale kinetic energy flux ((with or without our shock masking procedure in both cases). The profile are drawn for the three epochs of Fig.4. cascade. 7 In addition, the kinetic energy flux itself has important physical meaning. It bounds the dissipation rate of kinetic energy of gas motions into thermal energy (e.g., ZuHone et al. 2016; Zhuravleva et al. 2016) and into cosmic ray energy (e.g., Brunetti & Jones 2014;Miniati 2015). That energy flux also feeds the amplification of magnetic fields via small-scale dynamo action (e.g., Porter et al. 2015), although on smaller scales than we simulate here (e.g., Beresnyak & Miniati 2015). Therefore, in the remainder of this paper we will use the turbulent energy flux as our primary turbulence metric, rather than the velocity field or the kinetic energy to describe the turbulence in it903.
The energy flux across shocks provides an additional, specific and important ICM dynamical metric, since some fraction of this energy is dissipated into heat, while some it also feeds the generation of turbulence, as outlined above. In Fig. 10 (top left) we show the kinetic energy flux though shocks (f KE,shock , Eq. 1) and, for comparison, the kinetic energy flux of solenoidal (bottom right) and compressive (bottom left) filtered velocity fields (f KE,turb , Eq. 3) for the same slice as in Fig. 5. To better highlight the role played by shocks, in the same Figure (top right panel) we also show the kinetic energy flux of the total filtered velocity field, after masking the region tagged as shocked (Sec. 3.3.1).
The kinetic energy flux in the cluster is dominated by central, merger shocks, which process ∼ 10 40 − 10 41 erg/s per cell. However, in the innermost (1 Mpc) 3 cluster volume, downstream of the expanding merger shocks, the kinetic energy flux in the solenoidal component displays many large patches of high dissipation rate, with values of order ∼ 10 39 − 10 40 erg/s.
The radial profiles of the ratio of compressive to total kinetic energy flux, f KE,turb.comp /f KE,turb , are given in Fig. 11, both for the small-scale filtering procedure without masking shocks, and for the filtering procedure including shock masking when M 1.2 (Sec. 3.3.1). At all epochs, the flux ratio displays a marked increase with radius. With shocks included, the flux ratio ratio is ∼ 2 − 3 times larger (i.e. the relative energy flux of the compressive components is increased). When shocks are not included, the flux ratio is only ∼ a few percent in the central Mpc 3 volume ( 0.5 rvir) at all investigated epochs. This further justifies our use of enstrophy as a trustworthy proxy of turbulence in the following Section. Interestingly, close to the major merger epoch (z = 0.32 in the Figure) and at larger radii, the flux ratio jumps to ∼ 30% (∼ 40% if shocks are included), highlightling the significant generation of compressive turbulence triggered by the merger. We remark that, in order to better generalise this results, the analysis of a more extended set of clusters is necessary.

Enstrophy Analysis
In Sec. 3.4 we outlined properties of fluid enstrophy, ǫ = (1/2)ω 2 , that can be used effectively and efficiently to probe properties of solenoidal turbulence. Here we apply those tools to the simulation data for the it903 cluster. Fig. 12 illustrates the evolution since z = 1 of volumeaveraged enstrophy within the full (6.3 Mpc) 3 high-resolution volume (dot-dashed line) and within the smaller, cluster centered, (1.44 Mpc) 3 volume (solid line). It is obvious that signatures of the cluster dynamical history are much more evident if one focuses on a relatively small volume more closer to the virial size of the cluster. The z ∼ 0.3 major merger is quite obvious in the smaller volume, but not evident at all in the larger one. The enstrophy is expressed again as the normalized quantity,ǫ = ǫ · τ 2 0 , where the characteristic time, τ0, with τ0 = 1 Gyr. Then, √ǫ represents a representative number of turbulent eddy turnovers per Gyr. For com-parison, a characteristic eddy velocity, δv sol ∼ 100 km/sec and a characteristic coherence scale L ∼ 100 kpc, lead to 1/ √ ǫ ≈ 1 Gyr.
Here in the smaller volume the characteristic ǫ ∼ 10, implying eddy turn-over times ∼ 300 Myr. In addition, the mean enstrophy in the smaller, cluster-centred volume is several times larger than in the bigger volume. Although the mean enstrophy evolution in the larger volume does not reveal distinct events, it does show a slow, monotonic increase over time by roughly a factor of two, thus reflecting a gradual increase in turbulent energy per unit mass over time. The absence of clear signals for discrete events in this larger volume is due to the strong cluster concentration of the enstrophy evident in Fig. 2, or, analogously, concentration of the turbulent solenoidal velocity field shown in Fig. 5 or Fig. 10.

Comparison with solenoidal turbulent velocity field
One of our objectives in this discussion is to establish the degree of concordance in our simulations between enstrophy as outlined in section 3.4 and solenoidal turbulent velocity metrics as determined using methods outlined in Sec. 3.3. The ǫ values in Fig.  12 provide one simple test. In Sec. 4.1 we found characteristic turbulent solenoidal velocities v sol ∼ 80 km/sec (see Fig. 10). Those values lead toǫ ∼ 10 (with τ0 = 1 Gyr) provided length scales, ℓ = 2π/ √k 2 ∼ 100 kpc, which is quite consistent with the coherence scale analysis in Sec. 4.1.
A second valuable example of cross-comparison between the enstrophy and velocity analysis of turbulence comes through evaluation of the energy dissipation rate of the solenoidal turbulence, which we expressed in terms of the solenoidal turbulence velocity in Eq. 7 and in terms of the enstrophy in Eq. 8. In Fig. 13 we illustrate the spatial distribution in a 2D slice of the turbulence dissipation rate per volume at z = 0.32 from, on the left, the solenoidal turbulence (filtered) velocity field itself (Eq. 7 times ρ), and, on the right, the enstrophy field (Eq. 8 times ρ). There are minor difference in the details, but, on the whole the agreement is remarkably good.
In Fig. 14 we provide a comparison of the volume-integrated turbulence dissipation rate over time inside the (1.44 Mpc) 3 "cluster-centred" volume computed using the two formulations, again with δvL and L derived from the filtering analysis in Sec. 4.1 and with ℓ1 = ∆x. For the most part the two dissipation rate estimates agree to much better than a factor of two. The good agreement between these two independent estimates of the turbulent dissipation rate (also applying information on very different spatial scales) stresses once more that we are capturing reasonably well a turbulent-like cascade in its inertial range.
The only significant exception to the match between the two energy dissipation rates occurs early on, around z ∼ 1. At this epoch our zoom volume is mostly transected by large-scale converging motions on the proto-cluster. These motions have coherence scales of the order of the box size, which makes it impossible for our multi-scale filter to correctly disentangle bulk and turbulent component. Thus the filter identifies as turbulence even largescale shear motions outside of the proto-cluster, which did not have enough time to cascade down to the scale where enstrophy is measured. However, this problem quickly disappears as the cluster volume grows and the "zoom" region is mostly filled by the virial cluster volume.
In addition to the major merger event around z ∼ 0.3, there are other, recognizable turbulence evolution features visible in Figs. Figure 13. 2D slice at z = 0.32 (same as in Fig. 3) showing (left) the turbulence dissipation rate per unit volume based on the solenoidal turbulence velocity and Eq. 7 (or Eq. 3) and (right) the equivalent turbulence dissipation rate based on the enstrophy and Eq. 8.

Figure 14.
Turbulence energy dissipation rate integrated over the (1.44 Mpc) 3 comoving cluster-centred volume from the (filtered) solenoidal turbulence velocity field (red dashed curve) and from the enstrophy distribution (black solid curve). 12 and 14. These include a broad enstrophy peak in Fig. 12 between z ∼ 0.8 and z ∼ 0.6 with a maximum around z ∼ 0.7. This peak breaks into a pair of peaks near z ≈ 0.85 and z ≈ 0.6 in the turbulent energy dissipation rate. The enstrophy evolution plot Figure 15. Evolution of different energy fluxes for it903 (averaged within the central 1.44 3 Mpc 3 region: a) total kinetic energy flux over shock kinetic power (black); b) compressive kinetic energy flux over total kinetic energy flux (blue); c) compressive kinetic energy flux over shock kinetic power (red). The spread in the ratios at each redshift indicates ranges related to the inclusion or exclusion of shocked regions to compute the kinetic energy flux (the upper bounds being the estimates including shocked regions). exhibits shoulders at those times, but the peak is clearly offset in time.
Similarly, a close comparison of the enstrophy peak associated with the major merger shows that the turbulent energy dissipa-tion peaks around z ≈ 0.32, whereas the enstrophy itself peaks ∼ 1/2 Gyr later. Both behaviors are associated with the major merger event, but represent somewhat different dynamics. In particular, the turbulence dissipation rate actually peaks just before the closest approach of the two subcluster cores, when the enstrophy is most highly concentrated into the regions of highest gas density. The turbulent energy and its dissipation are also then focused into these regions. Dissipation rates outside the core regions remain relatively smaller. On the other hand, the sharp decrease in the mean enstrophy after core passage near z ≈ 0.32 is actually not so much a consequence of dissipation as of the strong outflows following the merger shocks generated during the event. In fact, as we point out in the section below, there is a net outflux of enstrophy from this central volume. There is, in addition, systematic decompression of the gas, which as equation 9 emphasizes, leads to enstrophy reduction.
The two earlier spikes in turbulent energy dissipation evident in Fig. 14 also correspond to merging activity, although minor mergers only. In each case there are brief intervals when turbulent motions are concentrated into the cluster core, which leads to sharp rises in the turbulent dissipation rate. The immediate impact on the total enstrophy budget is less significant in these cases.
Related information is illustrated in Fig. 15. It shows how ratios of various kinetic energy fluxes evolve with time in the inner, cluster-centred volume. For each redshift we show, both, the values obtained by removing the contribution from shocked regions (lower bound of each color) or by including M 1.2 shocks in the turbulent kinetic flux (upper bound). While the specific value of each ratio can change up to a factor ∼ 2 for most of the evolution, most of the time features are seen in both cases, and are in phase with the spikes in turbulent dissipation already noticed in Fig.14. In particular, each of the large spikes (z ≈ 0.85, ≈ 0.6 and ≈ 0.38) is associated also with the increase of the compressive kinetic flux, which reaches ∼ 10 − 15% of the total flux. Away from these spikes, the dissipation of turbulent motions is generally contributed by solenoidal motions at the ∼ 95% level. In this central volume the total kinetic energy flux is smaller than the kinetic power of shocks at most redshifts, with the exception of z 0.1 when the two becomes comparable in the absence of significant shock waves crossing the domain.

How, where and when is solenoidal turbulence generated in the ICM?
The previous analysis suggests that, while the budget of purely compressive turbulence is subject to uncertainties related to the presence of shocks, enstrophy provides consistent measures for the local and global solenoidal turbulence. Now we look at the processes that generate enstrophy described in Eq. 9 and 10, which allows a deeper understanding of the sources and amplification of turbulence in the cluster over time. Fig. 16 shows the time evolution of the enstrophy source terms defined in the above equations. Analogous to the enstrophy plots in Fig. 12, we apply a normalization factor τ 3 0 . The various ratios, ǫτ 2 0 /Fxτ 3 0 , then provide measures of the growth (or damping) timescale due to a given source term, Fx, measured in time units, τ0 = 1 Gyr. Indeed, by comparing Figs. 12 and 16 we can confirm that the timescales for the various source terms in this volume are τ ∼ 1 − 3 Gyr, consistent with the apparent evolution timescale of the enstrophy in Fig. 12.
As a reminder, the F adv term, analogous to ∇ · vρ in the mass conservation equation, measures the net enstrophy influx rate, while F stretch relates to the net rate at which vortex tubes are lengthening, F stretch > 0, or shortening, F stretch < 0. The Fcomp term identifies regions where enstrophy concentration correlates with ongoing gas compression. That can be reversible or not, depending on whether the gas compression is reversible or not (i.e., in shocks). The F baroc enstrophy source term identifies where nonvanishing cross products of density and entropy gradients align with the local vorticity (which may be expected downstream of nonplanar shock structures, for example).
We note three obvious properties of the individual source terms as revealed in Fig. 16. The first is that all the source terms, F adv , Fcomp , F stretch and F baroc averaged over this (1.44 Mpc 3 ) volume are roughly comparable. During merger events around z ∼ 0.7 and z ∼ 0.3, when the mean enstrophy is most rapidly increasing (Fig. 12), the F stretch term dominates, but only at most by a factor ∼ 2. During the major merger event around z ∼ 0.3, there is also a sharp peak in the compressive source term, Fcomp , in this volume.
The second notable outcome revealed in Fig. 16 is that, despite the fact that F adv is a conservative quantity, so that F adv dV = 0 over a closed volume, during most of the cluster history F adv > 0 in this volume. This highlights the fact that as a part of the accretion building the cluster, substantial enstrophy is also added into the central regions from outside. In fact, this enstrophy accretion is, most of the time, competitive with locally generated enstrophy growth (the other three terms). Enstrophy accretion into the cluster central region is identified below as at least partly a consequence of vorticity generated near the accretion shocks outside the cluster. We note, finally in reference to Fig. 16 that there is a brief period around z ∼ 0.25, following the major merger when both F adv and Fcomp become slightly negative in this volume. That behavior reflects an expansion and outflow of gas from the cluster during this same interval that is visible also in Fig. 1.
Key insights into the origins of the enstrophy within the cluster can be gained by examining the spatial distributions of Figure 18. Zoomed, 1.7 Mpc × 1.7 Mpc, central portion of Fig. 17, focusing on the "double relic" looking shock structure that has formed as a result of the ongoing merger. The units are as in Fig. 17. Identified shocks are again indicated by yellow contours. the enstrophy source terms. For instance, Fig. 17 shows in the same 6.3 Mpc × 6.3 Mpc, 2D slice used before at z = 0.32 (so during the major merger) the distribution of the four enstrophy source terms. The distribution of the enstrophy itself in this slice (not shown) is qualitatively similar to the distribution of filtered solenoidal velocity in the lower panel of Fig. 5. The advective source term, F adv , and especially the stretching source term, F stretch , have very roughly similar distributions to the enstrophy, as we might expect, since both are proportional to ǫ. In detail, however, all four distributions are quite distinct, as we should also expect, since each depends on unique dynamical behaviors.
It is also useful to compare the source term distributions to the distribution of shocks in the same region (yellow contours). There is, in this context, an anticipated association between shocks and the F baroc and Fcomp terms. Indeed, Porter et al. (2015), for example, demonstrated in compressible turbulence simulations that where the existing enstrophy is intense, and a strong shock exists, the compressive source term tends to be large. In this case, however, the relationship is complex.
This last point and its explanation are made clearer by examining two zoomed ∼ (1.5 Mpc 2 ) sub-sections of Fig. 17. In particular, we focus on the core cluster region between two merger shocks waves (Fig. 18), and on the the lower right accretion shock (M ∼ 10 2 ), well outside the cluster core (Fig. 19). Once again recall that these slices come from z = 0.32 during the major merger and when, according to Fig. 16, the enstrophy growth in the cluster centred volume peaks. In between the two internal merger shocks in Fig. 18 the F adv and F stretch terms show strong, distributed patterns of enstrophy development. In comparison, the F baroc and Fcomp terms are significant only very close to the merger shocks. Both show strong positive and negative peaks. In contrast both the F baroc and Fcomp terms associated with the accretion shock in Fig. 19 are predominantly positive, and their net contributions dominate those of F adv and F stretch in the area shown.
The picture that develops from this analysis is that the vorticity (enstrophy) generation within the cluster is a two step process: • Initial generation takes place in the cluster periphery in association with accretion shocks via baroclinic influences and compression. The compression contributions involve both external, accreting vortical motions and, as discussed in Porter et al. (2015) creation of enstrophy through shear stresses generated within complex shock structures. Some enstrophy also is accreted, especially along filaments.
• Subsequently, this enstrophy (solenoidal turbulence) is accreted into the cluster (F adv ), and in the innermost cluster regions the cluster enstrophy evolves especially by stretching, F stretch , but also by compressing those vortex structures initially generated in the cluster outskirts. Since those latter source terms are largest during merger activity, vorticity develops most strongly in those periods. Cluster-core-scale flows and shear are dominant contributors to those drivers.
To summarize this section, we surmise, using a combination of analysis tools, that the solenoidal turbulence found in the innermost regions of clusters is predominantly the result of injection and enhancement of accreting vortical flows at accretion shocks, followed by significant amplification by central advection, stretching and compression, particularly during merger events.The compressive and advective contributions can be at least partially reversible, whereas the stretching and baroclinic contributions generally are not.

Comparison to previous work
The  Gaspari et al. 2014;Schmidt et al. 2014) or structure functions (Vazza et al. 2011a;Miniati 2014Miniati , 2015. The character of the ICM gas flows and the presence of turbulence in our simulated cluster is consistent with other highresolution simulation studies in the literature (e.g., Dolag et al. 2005;Vazza et al. 2011a;Miniati 2014;Schmidt et al. 2014), with quantitative differences related to the specific analysis techniques. In the limited context of high-resolution grid-based simulations of non-radiative clusters, the pressure support from turbulence in the cluster core has, for example, been estimated in the range ∼ 5 − 15% of the gas pressure using the total, global gas velocity dispersion (e.g., Lau et al. 2009;Nelson et al. 2014), or, say ∼ 5 − 20% from large-scale, incoherent, but fixed scale velocity distributions, Vazza et al. (2011a). On the other hand, significantly smaller turbulent pressures tend to result from multi-scale filtering algorithms (e.g., ∼ 1 − 5%, Vazza et al. (2012) These numbers just highlight the underlying difficulty in disentangling large-scale from small-scale motions and assessing what is the best way to extract the purely turbulent component of complex and stratified 3D flows. At that level the issue is partly one of the intended meaning of the term "turbulence". Is the aim, for instance, to identify an energy reservoir for heat and/or cosmic rays, or is it simply to characterize deviations from hydrostatic equilibrium? On top of that issue is the inherent uncertainty in the meaning of the velocity variations due to inhomogeneity and stratification of the cluster, the usually unsteady cluster dynamical state (e.g., Lau et al. 2009;Vazza et al. 2011a), and the impact of additional physics such as radiative cooling and feedback (Vazza et al. 2013;Nagai et al. 2013;Battaglia et al. 2015) (neglected in this study). Non physical influences, including simulation numerical resolution and algorithms (Dolag et al. 2005;Battaglia et al. 2015) also are certain to influence these outcomes at some level. It is important to isolate as much as possible the issues, so to address them as cleanly as possible.
The closest analogy to the present study is the "Matryoshka" simulation by Miniati (2014), who, using a simulation strategy similar to ours, followed the evolution of a ∼ 10 15 M⊙ cluster using ≈ 10 kpc resolution inside the virial radius (which roughly spans a ∼ 6 times larger dynamical range) . The "Matryoshka" cluster gas velocity field was examined using structure functions and Hodge-Helmholtz decomposition to follow global evolution of solenoidal and compressive turbulence properties separately within roughly the cluster core and inside the virial radius (Miniati 2014(Miniati , 2015. Even though the "Matryoshka" cluster was ∼ 10 times more massive than it903, the general history of the cluster and the character of evolving turbulence are roughly similar. Miniati (2015) found compressive turbulence to be in the range ∼ 0.2 − 0.4 of the total turbulent energy in the central (Mpc) 3 region.
Our results are roughly similar, but suggest a significantly smaller contribution from compressive turbulence because of our filtering of shocked regions. We found in it903 that the compressive dissipation rate is in general ∼ 5% of the total dissipation rate for most of the cluster lifetime, with a spike of ∼ 15% during the major merger event (the jump in the compressive dissipation is larger in cluster outskirts, where it reaches ∼ 30% close to the merger). This translates roughly into a 15 − 30% energy fraction in compressive turbulence, assuming both components have identical outer scales. Additionally, we measure a steeper spectral distribution of energy in the solenoidal velocity component than in the compressive component.
We surmise that the reduced compressive turbulence role we found is mostly due to our procedure of extracting turbulent motions (both solenoidal and compressive) from the filtered, smallscale uncorrelated velocity field (Sec. 3.2), rather than from the total velocity field, but also from our removal of contributions from stronger shocks to the turbulent motions on the grounds that most of those strong shocks are not directly involved in the uncorrelated motions.
Even if the procedure of imposing a fixed filtering scale for the entire volume, as in Miniati (2014), mostly removes the largescale laminar velocity component (see also Vazza et al. 2011a), that procedure cannot fully account for the variation in the turbulent coherence length in the stratified cluster atmosphere on intermediate scales. Our previous tests in Vazza et al. (2012) showed, indeed, that the turbulent energy budget in the cluster centre is usually reduced by a factor ∼ 2 when a spatially varying filtering scales is adopted.

CONCLUSIONS
Understanding the dissipation of turbulent energy is key to understand the heating of the plasma, the acceleration of cosmic rays in the ICM, as well as the growth of intracluster magnetic fields (e.g. Subramanian et al. 2006;Brunetti & Jones 2014;Miniati 2015).
While current X-ray line spectroscopy can provide only upper limits on the chaotic motion velocities in relatively bright cluster cores (e.g., Sanders et al. 2011;Pinto et al. 2015) 8 , future X-ray satellites with superior spectral resolution (e.g. ATHENA) should be able eventually to detect directly the driving-scale turbulent motions in the ICMs of multiple clusters (Nagai et al. 2013;Ettori et al. 2013;Zhuravleva et al. 2013;). In the meantime, turbulent motions in the ICM induce moderate pressure fluctuations that may be detected in X-rays (e.g., Schuecker et al. 2004;Sanders & Fabian 2012;Churazov et al. 2012;Gaspari et al. 2014;Zhuravleva et al. 2014), or through the S-Z effect (e.g., Khatri & Gaspari 2016). 8 The Hitomi satellite in its short life did successfully measure velocity profiles for the Perseus cluster (The Hitomi Collaboration 2016) Numerical simulations of the ICM are fundamental to assessing the real impact of ICM turbulence on all the above. In this work, we focused on the analysis of the connection between accretiondriven shock waves and turbulent motions in the ICM. In particular, we explored both the local and the statistical causal connections between shocks and the emergence of solenoidal and compressive turbulent motions during a simulated cluster lifetime.
Our main conclusions from this study can be summarized as follows: • Gas flows in the ICM are characterised by a turbulent behaviour across a wide range of scales, roughly consistent with a Kolmogorov-like model. However, these flows are mixed with larger-scale regular (correlated) velocity components for scales 0.1 − 1 Mpc and are punctuated by small-scale velocity perturbations due to shocks, which makes it difficult to isolate accurately uncorrelated turbulent fluctuations of the flow at most scales (Sec. 4.1).
• Using Hodge-Helmoltz decomposition within domains distributed across our simulated cluster, we measure dominant solenoidal velocity fields everywhere within the cluster and at most epochs (with the exception of high redshift epochs, when the cluster is still forming and is far from a virialised state). The solenoidal component makes 50−80% of the amplitude of the total velocity field at most epochs and scales (Sec. 4.1).
• The kinetic energy dissipation rate of the small-scale velocity field is a powerful tool to measure the ratio of compressive and solenoidal motions in a nearly scale-independent way. The dissipation in compressive modes only accounts for a few percent of the total turbulent dissipation rate in the central ∼ Mpc 3 volume. This can increase to about ∼ 15% in the central Mpc 3 during major merger events, and to ∼ 30% in cluster outskirts (Sec. 4.1).
• Vorticity and enstrophy are trustworthy proxies of the dominant solenoidal turbulent component. In particular, the volumeintegrated dissipation rate of solenoidal turbulence and of enstrophy are very well correlated in the 0 z 1 redshift range, and they show remarkably similar spatial patterns (Sec. 4.2).
• For the first time, we apply the Navier-Stokes formalism to analyse in detail how enstrophy evolves in the simulated ICM, by decomposing its growth rate into advective, stretching, compressive and baroclinic terms (Sec. 4.2).
• At accretion shocks baroclinic generation of enstrophy along with enstrophy enhancement during flow compression are the most important source terms of enstrophy. In cluster interiors vortex stretching dominates the growth of enstrophy, although advective concentration of enstrophy and, especially during mergers, enstrophy enhancement through compression can be comparable. Merger shocks largely seed the enstrophy enhanced by vortex stretching and advective concentration in the cluster interior (Sec. 4.2).
The study of this first cluster of the ISC sample showed how rich is the complexity of simulated ICM turbulence, even in this rather restricted physical setup. Our analysis suggests that a careful combination of filtering techniques is mandatory to identify all major components of the turbulent energy budget reliably, and to give them a physical meaning as a function of scale. Through the extensive analysis of our full ISC sample in planned follow-up work it will be possible to generalise the results obtained for this first cluster in a more robust statistical way.

A1.1 Spatial filtering for turbulence
Similar to Vazza et al. (2012), we ran several tests of our iterative filtering procedure for turbulence (Sec.3.2) over control turbulent boxes with predefined outer scales and slopes for the power spectrum. In the example given here, we analysed a purely solenoidal velocity field on a 100 3 grid, first generating a vector potential and then computing its curl. The vector potential was drawn from 3 k 50 and had a power-law spectrum with a slope Ev(k) ∝ k −5/3 . We then ran our algorithm with the same set of parameters as the analysis in the paper and checked whether the input velocity field was correctly recovered. The panels in Figure 1 present the performance of this test, comparing the power spectrum and the 2 nd order structure function of both input and reconstructed velocity field using our filtering procedure 3.2.
The results are representative of the performances of the filter in this kind of test, which also applies by construction when regular fields on large scales are imposed on the setup (see, e.g. Vazza et al. 2012, for similar tests)). Both the outer injection scale of the flow and the spectral behaviour of the fluctuating field are accurately recovered by our filtering procedure. The slope of the power spectrum and of the structure function are recovered within ∼ 5 − 10% accuracy. In these tests, typically only ∼ 10% percent of the kinetic energy in the fluctuating field is not recovered after the filtering procedure. That is due to a misidentification with a large-scale smooth flow, which can occur when the outer correlation scales of the flow are of the order or larger than the computing domain. (Hence, they mimic a large scale regular component.) On the other hand, the kinetic energy flux within the cascade is captured extremely well. For the main goals of our ICM analysis this small difference plays only a minor role, and for the purpose of our analysis in the main paper this filtering technique is suitable for capturing the most important small-scale turbulent features of the flow.

A1.2 Hodge-Helmoltz modes decomposition
We tested the ability of the Helmholtz-Hodge ("HH") procedures employed in Sec.3.3 to decompose blended solenoidal/compressive velocity fields correctly and in comparison to a simple, "straw man" alternative. The basis for the HH procedure is that a vector field with suitable asymptotic/boundary behaviors can be expressed in a form v = −∇φ + ∇ × a = vc + vs (e.g., Arfken & Webber 1995). In particular, assuming the vector field, v, has a Fourier transform, V k , one can obtain the two components, vc and vs. The algorithm applied in 3.3 (and which we call here Method 1) assumed explicitly that V k = i kΦ + i k × A = V k,c + V k,s (where capital letters in this discussion generally refer to Fourier transforms of lower case spatial fields); that is, that the Fourier transform, V k,c = k( k · V k )/|k| 2 , can be found by projection of V k onto k and that V k,s = − k × ( k × V k )/|k| 2 can be found normal to k. The alternate approach (which we call Method 2 here), instead applied fourth order finite difference methods to estimate d = −∇ 2 φ = ∇ · v, then obtained Φ = D/k 2 . The inverse Fourier transform yielded φ, while vc = −∇φ was obtained from φ using fourth order finite differences. We note that use of high Figure 1. Top: velocity field power spectra for the input velocity of our test (black) and for the output velocity field as reconstructed by our iterative filter (red). Bottom: 2nd order longitudinal structure functions for the same test field, with identical meaning of colors. order differences for spatial derivatives is especially important in attempts to extract compressive component variations close to the Nyquist frequency, since those difference algorithms can construct exact derivatives of polynomials up to the order of the differencing. Equivalent procedures to those employed to find vc could be used to find vs directly, of course. In these tests, however, we found vc explicitly, then obtained vs as vs = v − vc.
To carry out the HH tests we first constructed idealized 3D velocity fields in a cubic box of arbitrary length, L, with both compressive, vc, and solenoidal, vs, velocity components using discrete Fourier sums; v( x) = V k cos( k · x + ψ k ), with 0 ψ k 2π selected from random deviates. For simplicity we aligned all the component wave vectors, k along the x axis; that is k = kxx, with |k| = 2πn/L and |n| > 0. Then, of course, the velocity fields were Figure 2. Helmholtz decomposition into compressive (left) and solenoidal (right) component magnitudes of a 3D velocity field plotted along x using Method 1(red) and Method 2 (green), as explained in the text. The constructed, "input" components are shown in black. The different rows correspond to different sized portions of the box: full 256 3 zone, periodic box; middle: 100 3 zone sub-volume; top: 50 3 zone sub-volume. periodic on length, L, while each compressive (solenoidal) Fourier component satisfied k × V k = 0 ( k · V k = 0). The compressive and solenoidal fields were given distinct, non-vanishing power spectra over distinct ranges in kx, vc and vs can be obtained from Φ and A, respectively, as outlined in the previous paragraph. Note that there is no power in n = 0 modes.
The constructed velocity field was evaluated on 256 3 uniform spatial grid points spanning L. Although the constructed vector fields were periodic over the full, 256 3 domain, velocity fields in our cluster simulations generally are not periodic in domains of interest. Therefore, we tested the accuracy of both HH algorithms to recover correct velocity information from nonperiodic, 100 3 and 50 3 cell subvolumes.
An example test velocity field is illustrated by the solid, black lines in Fig.2. The compressive (solenoidal) component is on the left (right). The bottom panels show the velocity distributions across the full domain, while the middle (top) row show them in a 100 3 (50 3 ) subvolume. The case shown in Fig.2 is for a compressive component with modes spanning 2 kx 128 with a E(k) ∝ k −2 spectrum, and a solenoidal component with a ∝ k −5/3 spectrum in modes spanning 4 kx 32. Fig. 2 shows the outcomes of the Method 1 and Method 2 HH decompositions of this velocity field with Method 1 shown in red and Method 2 in green. Both on the full, periodic volume and the non-periodic subvolumes the Method 1 reconstructions are almost identical to the input velocity field. That is, this method proved to be quite accurate. In contrast Method 2 clearly misses substantial power in the compressional component, despite the use of high order difference algorithms to estimate spatial derivatives. In conclusion, then, the Method 1 HH decomposition that depends entirely on FFTs is a reliable approach to obtain the compressive and solenoidal velocity components, even when the domain of interest is not periodic. In contrast, the use of finite differences to avoid FFT extraction of non-periodic spatial gradient information was not successful.