{'fraction_non_alphanumeric': 0.03393098129940235, 'fraction_numerical': 0.020397146712936187, 'mean_word_length': 4.628974498100923, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exploring the relation between turbulent velocity and density fluctuations in the stratified intracluster medium

The dynamics of the intracluster medium (ICM) is affected by turbulence driven by several processes, such as mergers, accretion and feedback from active galactic nuclei. X-ray surface brightness fluctuations have been used to constrain turbulence in galaxy clusters. Here, we use simulations to further investigate the relation between gas density and turbulent velocity fluctuations, with a focus on the effect of the stratification of the ICM. In this work, we studied the turbulence driven by hierarchical accretion by analysing a sample of galaxy clusters simulated with the cosmological code ENZO. We used a fixed scale filtering approach to disentangle laminar from turbulent flows. In dynamically perturbed galaxy clusters, we found a relation between the root mean square of density and velocity fluctuations, albeit with a different slope than previously reported. The Richardson number is a parameter that represents the ratio between turbulence and buoyancy, and we found that this variable has a strong dependence on the filtering scale. However, we could not detect any strong relation between the Richardson number and the logarithmic density fluctuations, in contrast to results by recent and more idealised simulations. In particular, we find a strong effect from radial accretion, which appears to be the main driver for the gas fluctuations. The ubiquitous radial bias in the dynamics of the ICM suggests that homogeneity and isotropy are not always valid assumptions, even if the turbulent spectra follow Kolmogorov's scaling. Finally, we find that the slope of the velocity and density spectra are independent of cluster-centric radii.


Introduction
In the current paradigm of hierarchical structure formation, accretion of smaller structures and mergers play a key role in the evolution of galaxies and galaxy clusters. Gas that falls into dark matter halos is heated, mostly via shocks or turbulent dissipation, to the virial temperature of a cluster (10 7 − 10 8 K). This results in a fully ionised X-ray emitting intracluster medium (ICM). During the process of cluster formation, mergers and accretion generate hydrodynamic instabilities that create turbulence. Observations and simulations have investigated the properties of turbulent motions in the ICM. In this work, we extend previous theoretical efforts by examining the competing roles of stratification and large-scale infall motions.
There is observational evidence that supports the idea of a turbulent ICM. Examples of observed features used to measure ICM turbulence include X-ray surface brightness fluctuations, interpreted as density fluctuations induced by turbulence, or pressure fluctuations obtained from X-ray maps (Schuecker et al. 2004, Sanders & Fabian 2012. Another method to infer information about turbulence is the study of emission lines in the X-ray band. Sanders et al. (2010) placed limits on the turbulent broadening of the emission lines in cool-core clusters. One of the main findings so far for the internal kinematics of the hot gas was obtained observing the cool-core Perseus cluster with the Hitomi satellite (Hitomi Collaboration et al. 2016). Those authors found a 1D velocity dispersion ∼ 160 km s −1 , on a scale of L ∼ 50 kpc in the core of the Perseus cluster. Similar subsonic turbulent velocities are also suggested by studies of a relatively cold ICM in cool-core clusters (e.g. Olivares et al. 2019;Russell et al. 2019). Indeed, if the cold gas originates by hot ICM cooling, the two phases share the same kinematics (Gaspari et al. 2018). Turbulence requires a large Reynolds number, which implies a low viscosity. This is also indicated by the evidence of sharp features connected to Kelvin-Helmholtz instabilities and cold fronts in the ICM (Markevitch & Vikhlinin 2007), which otherwise would be suppressed due to thermal conduction and viscosity (Roediger et al. 2013, Wang & Markevitch 2018. Finally, it is possible to exploit the thermal Sunyaev-Zeldovich effect to infer turbulent pressure maps (Khatri & Gaspari 2016).
Turbulent motions can indirectly affect other observables of galaxy clusters because they affect non-thermal components, such as magnetic fields and cosmic rays. Combining X-ray and radio observations, Eckert et al. (2017) and Bonafede et al. (2018) found that the brightness fluctuations correlate with diffuse radio emission, suggesting that turbulence might be linked to the acceleration of non-thermal particles. Rebusco et al. (2006) investigated the effect of turbulent diffusion on the metal abundance profile. Moreover, turbulence might play a role in the measurement of cluster masses, as it creates additional kinetic pressure that affects the hydrostatic equilibrium typically Article number, page 1 of 14 arXiv:2112.00774v1 [astro-ph.CO] 1 Dec 2021 assumed to apply to the ICM. In fact, observational and numerical work has confirmed that the ratio between turbulent and thermal energy typically varies between ∼ 10 − 30% (e.g. Angelinelli et al. 2020;Eckert et al. 2019;Ettori et al. 2019;Lau et al. 2009). Since precision cosmology depends on very accurate total cluster mass, turbulent pressure profiles need to be well constrained. Finally, turbulence has been proposed as a source of heat to prevent development of catastrophic cooling flows (Zhuravleva et al. 2014a). The most promising scenario there involves turbulent feedback induced by active galactic nuclei. The central gaseous medium is heated by the interaction with buoyantly rising bubbles that are created by jets launched from the central black hole. While such feedback flows can be energetically sufficient, the process by which the mechanical energy is transferred to the ICM and thermalised is still debated (Brüggen & Kaiser 2002;De Young 2010;Gaspari et al. 2013a,b;Randall et al. 2015;Fabian et al. 2017).
Fully cosmological, or more idealised simulations have been used to investigate turbulence in galaxy clusters (Dolag et al. 2005;Kim & Ryu 2005;Wittor et al. 2017;Vazza et al. 2009Vazza et al. , 2011bVazza et al. , 2017. Zhuravleva et al. (2014b) and Gaspari et al. (2014) studied the relation between density and velocity perturbations, Vazza et al. (2012) investigated the turbulent diffusion coefficient and Angelinelli et al. (2020), Nelson et al. (2014), Lau et al. (2009) constrained the non-thermal pressure contribution coming from turbulent motions. One of the difficulties in determining ICM turbulent quantities is the disentanglement of bulk and turbulent motions. Several algorithms have been developed for this purpose. For example, Lau et al. (2009) estimated ICM turbulent velocities as the residual with respect to the velocity field averaged over spherical shells. Alternatively, one can estimate the turbulent motions by interpolating the original 3-D velocity field in order to map the local mean field and to identify turbulent velocity fluctuations on scales smaller than the interpolation scale (Vazza et al. 2009(Vazza et al. , 2011b. As another alternative, a multi-scale iterative filtering approach has been implemented in, for example, Vazza et al. (2012), Angelinelli et al. (2020), ?. In the present work, we apply a simpler approach, that is by using a fixed scale filtering method, based on previous results by our group.
Describing the turbulent nature of the ICM is not trivial. Astrophysical turbulent flows are often described in terms of Kolmogorov's theory (Kolmogorov 1941), which assumes homogeneity and isotropy. However, the ICM is a stratified plasma close to hydrostatic equilibrium, where buoyancy may well change the character of the turbulence. In particular, the relation between density and velocity fluctuations, often used to estimate the turbulence strength in the ICM (e.g. Zhuravleva et al. 2014b) might vary in the presence of a varying stratification (and forcing mechanism). The dimensionless Richardson number, given by the square of the ratio of the turbulent eddy turn-over timescale to the buoyancy oscillation timescale (or, equivalently, by the ratio of the buoyancy to inertial forces), is a measure of the importance of the stratification in the turbulence dynamics. It is defined as: where v l is the characteristic turbulent velocity on a scale l, and N BV is the Brunt-Väisälä frequency: Here, g(r) is the gravitational acceleration, P(r) and ρ(r) are the pressure and density, respectively, and γ is the adiabatic index of the ICM. Thus, Ri 1, implies short turbulent eddy turn-over times compared to buoyancy times, suggesting homogeneous, isotropic turbulence unaffected by density stratification. On the other hand, when Ri > 1 the buoyancy times are shorter than the turbulent eddy times. When the buoyancy force becomes dynamically important, it suppresses radial motions and it leads to preferentially azimuthal turbulence (Brethouwer et al. 2007).
There are previous works investigating turbulence in a stratified ICM-like medium (Zhuravleva et al. 2014b;Gaspari et al. 2014;Valdarnini 2019;Shi & Zhang 2019;Mohapatra et al. 2020). However, they focused on idealised plasma flows, or were limited to the innermost regions of galaxy clusters (≤ 500 kpc). Previous observational and numerical work investigated the possibility of inferring the statistics of turbulence in the ICM from gas density and temperature perturbations, which can be recovered from existing X-ray observations (Arévalo et al. 2016;Churazov et al. 2012;Zhuravleva et al. 2018). However, due to the finite sensitivity of X-ray telescopes, this technique is still limited by the large number of photons required to sample large enough cluster scales, with sufficient spatial detail. In order to further compare the results with observational data, past analyses of simulations were restricted to the virial region of the cluster. Moreover, idealised simulations are also bound to miss substructures and filaments, which can bias the analysis compared to the statistics derived from more realistic, dynamically formed clusters.
In this paper we explore the character of the turbulence in realistic, stratified ICM generated in high-resolution, cosmological simulations of galaxy clusters. We mainly focus on the relation between density and velocity fluctuations, their dependence on the Richardson number (see Mohapatra et al. 2020), and the anisotropy of the turbulent velocity field. This paper is structured as follows: in Section 2 we describe the simulated cluster sample and the numerical method used to estimate the relevant variables of a turbulent ICM; in Section 3 we discuss the outcomes of the analysis, looking also for a comparison with previous work; in Section 4 we summarise our findings, highlighting some observational implications.

The Itasca Simulated Cluster sample
We analysed a sample of galaxy clusters from the 'Itasca Clusters' set of simulations. Each simulation was carried out at high spatial resolution using ENZO, a (parabolic) cosmological numerical code for magneto-hydrodynamics (Bryan et al. 2014). A key feature of the code is its Adaptive Mesh Refinement (AMR) capability, which enables a large spatial and temporal dynamical range.
Here we applied a pre-determined set of nested grids in order to apply a constant, comoving grid resolution. Starting from the root grid, which covers the entire volume with a coarse, uniform grid, we placed finer grids as soon as interesting regions started to evolve. In the Itasca suite, we just forced the code to refine the 100% of the innermost zoom-in region, up to the highest allowed AMR level. If the root grid spacing was ∆x, then the spacing of a refined patch at level l was ∆x/r l , where r is the integer refinement factor. The simulations presented here employed a purely hydrodynamical fluid solver based on the Godunov Piecewise Parabolic Method (PPM) scheme. (Colella & Woodward 1984).
The cosmological Itasca simulations assume a WMAP7 ΛCDM cosmology (Komatsu et al. 2011), with Ω B = 0.0445, Ω DM = 0.2265, Ω Λ = 0.728, Hubble parameter h = 0.702, σ 8 = 0.8 and a primordial index of n = 0.961. All runs were non-radiative and did not include feedback from star-forming regions or active galactic nuclei. Even though the effects of nongravitational heating are small on the 100 kpc scale compared to the heat caused by mergers (Valdarnini 2019), the combination of cooling and feedback leads to a higher number of substructures with different temperatures and densities could create small-scale turbulent motions. In this sense, such structures are absent in our simulations and their presence in the observed galaxy clusters might alter gas density and temperature estimates. Our simulations do not include magnetic fields even though the interplay between turbulence and dynamo amplification is well-established in galaxy clusters (e.g. ????). However, the magnetic pressure is small, considering that the typical plasma β parameter (i.e. the ratio between thermal and magnetic pressure, β) is very large, β ∼ 10 2 − 10 3 . This means that we do not expect any relevant modification in the density or velocity structure of the ICM by the (neglected) presence of magnetic fields, except on very small scales (≤ 30 kpc) where the tension of amplified magnetic field lines is significant (?). Nevertheless, the gas dynamics on these small scales is not our primary concern and we proceed with a purely hydrodynamical approach.
In brief, the Itasca simulations were run according to the following procedure. Firstly, independent cosmological boxes were simulated in order to select the most massive objects (clusters) in the volume. Initial conditions were generated separately for each simulation at redshift z = 30. The spatial resolution, on the coarse grid, was L 0 = 110 h kpc −1 ≈ 157 kpc and the dark matter (DM) mass resolution was m DM = 8.96 · 10 7 M . A second more refined grid was created, centred on the selected cluster formation region, which covered the innermost ≈ 31 Mpc. In this case, the comoving spatial resolution reached 55h kpc −1 ≈ 78.4 kpc and the DM mass resolution was m DM = 1.12 · 10 7 M . Finally, inside the central (L 0 /10) 3 ≈ 6.3 Mpc 3 volume of each box, a fixed further refinement was enforced, meaning that 100% of cells were refined. This volume was large enough to include the virial radii of most of the clusters, while the spatial resolution was increased to 13.8 h kpc −1 ≈ 20kpc.
The full Itasca cluster sample consists of 20 galaxy clusters and has been used to investigate the dynamic and thermal properties of cluster ICMs. In this work, we used 8 Itasca clusters at z = 0 in the mass range 0.54 · 10 14 M < M 200 < 3.32 · 10 14 M . For further information on the Itasca cluster sample, we refer to Vazza et al. (2017), Wittor et al. (2017), Angelinelli et al. (2020).

Filtering method and numerical analysis
In order to disentangle bulk from turbulent motions, we used a tailored filtering method. Vazza et al. (2012), Vazza et al. (2017) have shown that galaxy clusters might present a distribution of turbulent injection scales, since ICM turbulence can be driven by many processes (e.g AGN feedback or accreted clumps). However, Vazza et al. (2018) utilising the Itasca clusters found that typical scales from turbulence injected by accretion, are in the range ∼ 200 -400 kpc. For this reason, we choose L = 300 kpc as a filtering scale in the following analysis. Nevertheless, we tested the robustness of our analysis by using other filtering scales.
Our filtering method estimates the turbulent fluctuation (δv, δρ) as follows: where v and ρ are the local velocity and density and ρ L and v L are the means estimated on cubic boxes of (linear) size L, centred on the cell to filter. In the case of the velocity, we use a density-weighted mean. Fig. 1 shows maps of density and velocity fields, before and after we performed this algorithm in a cluster of our sample. We then estimate the density contrast, δρ/ ρ L and the normalised velocity fluctuation, δv/c s , where c s is the sound speed at a given radius, obtained from the (azimuthally averaged) temperature profile. We investigate the relation between these two quantities to directly compare to the work by Zhuravleva et al. (2014b).
In the following, we calculate the statistical quantities (for instance the variance of the density contrast) in cubic boxes of linear size 600 kpc, located at different radii and different directions. In this way we map the whole simulation box to assess the spatial variation of the turbulence properties. We also performed a similar analysis by dividing the main box in spherical shells, with similar but somewhat less clear results. This can be understood by the azimuthally asymmetric ICM dynamics, generated by accretion events along the large scale cosmological filaments. Because of this we only discuss the sub-boxes analysis in the rest of the paper.

Density and velocity fluctuations
We adopt as a measure for the typical density fluctuation the standard deviation of the probability distribution function (PDF) of the density contrast δρ/ ρ , which we indicate in the following with σ ρ , computed in every box. In order to limit the contamination from shock compression and substructures, we eliminate the cells corresponding to the top 5% of the density contrast distribution in every box, following earlier work concerned with X-ray modelling (e.g. Roncarelli et al. 2013). For testing purposes, we also experimented with other thresholds, obtaining: σ 2 ρ,99 = 0.26, σ 2 ρ,95 = 0.22 and σ 2 ρ,90 = 0.20, which refer to thresholds of 99%, 95%, 90%, respectively.In a similar way, the mean normalised velocity fluctuation, σ v , is estimated as the standard deviation of δv/c s , defined in the previous section, using the same cells used to calculate σ ρ .
Finally, to compare our results on stratified turbulence with previous work, we also consider the logarithmic density fluctuations, s = ln(ρ/ ρ L ) and calculate the standard deviation in every box, σ s . Fig. 2 shows, for each cluster of the sample, the probability distribution function of the variables we use in this work: the gas density contrast, δρ/ρ, the logarithmic gas density fluctuations, s = ln(ρ/ρ), and the velocity contrast δv/c s . The PDFs are obtained from the fluctuations within the virial radius (defined as, R 200 , i.e. the radius enclosing a 200 overdensity with respect to the cosmic critical density).

Richardson number
The Richardson number, Ri, is central to our study, but, in contrast to some previous work in which the authors considered Ri as an input parameter of the simulation (Mohapatra et al. 2020, Shi & Zhang 2019, we had to compute it a-posteriori following Eq. 1. We obtained the radial profile of the Richardson number in 600 kpc wide boxes, starting from the gravitational acceleration and gas pressure and density radial profile in the box that we used to estimate the Brunt-Väisala frequency , N BV (Eq. 2). We also evaluated the turbulent velocity profile in the box, computing a root mean square after the excision of the densest density fluctuations (see Sec. 2.2.1). For the sake of simplicity, we fixed the spatial scale in Eq. 1 as the size L of the filtering scale (see Sec. 2.2). In order to characterise each box with its Richardson number, we averaged the Ri radial profile.
Furthermore, we obtained the radial profile of the Richardson number over the whole volume of each galaxy cluster. The top panel of Fig. 3 shows the Richardson number profile in our cluster sample for a filtering scale L = 300 kpc (top panel), together with the Brunt-Väisala frequency profile (middle panel). In all clusters, except one, buoyancy is always found to be dominant (Ri 1) within the virial radius, while the turbulent cascade overcomes gravity only in the outskirts. Moreover, we find that the Richardson number has a strong dependence on the filtering scale L (Fig. 3, lower panel). It nearly follows the Kolmogorov scaling Ri ∝ L 4/3 . Consequently, a local estimate of Ri must always be referred to a turbulent scales.
As noted above, in the remainder of this work we restrict ourselves to the fixed filtering scale of L = 300 kpc; this is motivated by our previous works on the Itasca dataset. For example, in Vazza et al. (2017) and Angelinelli et al. (2020) we measured, with an iterative multi-scale filtering algorithm, that the kineticenergy weighted distribution of cells in the volume of our clusters have a local turbulent outer scale ≤ 400 kpc, with little dependence on the dynamical state of the host cluster. Based on these earlier studies, we can conclude that L ≥ 500 kpc appears to be too large based on the actual kinematic analysis of these systems, while L ≤ 200 kpc scales are too small to capture the total turbulent kinetic energy budget developed in the ICM of these simulated systems.
In order to quantify the degree of anisotropy, possibly caused by the buoyancy, we compute the parameter β, where v t and v r are the tangential and radial components of the turbulent velocity distribution (that is, the r.m.s. values computed in 600 kpc wide boxes, with the filtering scale again set to L = 300 kpc), respectively. The anisotropy parameter can span the range −∞ < β ≤ 1. Positive values mean that the dynamics is dominated by the radial motions, whereas, when β < 0 the dominant components of the velocity field are those perpendicular to the radial direction. We refer the reader to Vazza et al. (2018) for a further analysis of the anisotropy parameter in the Itasca cluster sample.

Power spectra
In this section, we compare the density and velocity power spectra at various cluster radii. We measure the power spectra of the ICM velocity field in large boxes, ≈ 2.5 Mpc (linear size) wide, located at different radii using a Fast Fourier Transform, with a periodic domain for simplicity (e.g. see Vazza et al. 2011b, for a discussion): where The gas density fluctuations, |ρ(k)| 2 , are obtained in a similar way, that is by computing the Fourier transform of the real space density distribution. In this case we did not perform any filtering, since filtering with a fixed spatial scale would just suppress the power on k ≤ 2π/L scales, but we select the inertial subrange highlighted with the orange bars ( Fig. 4): 150 kpc l 1 Mpc. The small scale limit is being fixed to avoid measuring spectral shapes in a regime already affected by the numerical dissipation of the PPM hydro scheme. The normalisation of the power spectra in Fig. 4 has been rescaled to compare the slopes at various radii. In general, the normalisation of the density spectra increases outside the virial radius, except for the power spectrum computed in the outermost region of the cluster which shows a drop. This increasing trend is likely due to the larger clumping factors in the outskirts (Angelinelli et al. 2021).
We evaluated the slope of each k-bin within the inertial subrange and we computed an average value for a more robust characterisation of the power spectrum. For the velocity power spectra we evaluated the slope in each bin of each velocity component and then we averaged over the all slope values. Finally, we averaged the slopes of the power spectra performed on boxes at the same radius, computing an inverse-variance weighted mean, to investigate the radial trend of the slope of the power spectrum, compared to the ∝ k −11/3 trend expected for the Kolmogorov model in three dimensions.
As a preliminary check on the reliability of our cluster simulations for this study, we computed the statistics of gas fluctuations. As previously explained, some spurious features that can affect our numerical analysis might be present in our galaxy clusters (e.g. shocks, clumps or filaments), as well as to some degree also in real observations. In order to excise most of such features not associated with true turbulence-driven fluctuations, for each radial shell we computed the probability distribution function of density and excluded the 5% densest cells. In order to have uniformly filled maps, we replaced the original density values of these cells with the median density within the entire radial shell. This method has been implemented in previous works (see Vazza et al. 2011a, for an example). Considering that bremsstrahlung is the dominant radiative mechanism at the typical temperature of plasma in galaxy clusters (10 7 − 10 8 K), we estimated the surface brightness by projecting S X = n 2 T 1/2 (where n and T are the plasma density and temperature in each cell) along multiple lines of sight. Fig. 5 shows the projected X-ray surface brightness in one of our simulated clusters before and after removing the clumps. We can thus estimate the surface brightness fluctu-ations, dividing the X-ray images by their radial average value within each radial bin. This method is commonly used in observational studies, even though it does not consider any azimuth asymmetries of the ICM. Thus, we obtain the 2D power spectra of X-ray fluctuations on 1 Mpc wide boxes centred on the gas density peak of the cluster. The characteristic amplitude is computed as follows Churazov et al. (2012): where P 2D (k) is the 2D power spectrum and k is the wavenumber in kpc −1 units. The distribution of surface brightness fluctuations in our cluster sample is shown in Fig. 6. Eckert et al. (2017) computed the distribution of the amplitude in surface brightness fluctuations for 51 galaxy clusters, and looked for the possible correlation between larger fluctuations and the presence of nonthermal radio emission. The results can be used to investigate the connection between turbulence and the re-acceleration of relativistic electrons via Fermi II mechanism (e.g. Brunetti & Jones 2014, and references therein). The red line in Fig. 6 represents the mean value of the X-ray surface brightness fluctuations in Eckert et al. galaxy cluster sample, which is 0.9. We observe that our fluctuations are by a factor ∼ 2 larger. We should mention that our statistics are not affected by observational limitations. The surface-brightness fluctuations can be reconstructed for the entire cluster volume with equal fidelity and ignoring such limiting factors as the number of emitted Xray photons. Real observations, however, are affected by the uneven effective spatial resolution, and by the fixed sensitivity to Xray brightness fluctuations. Moreover, it is not straightforward to implement exactly the same procedure for the masking of pointlike X-ray sources in observations, which can bias the observed amplitude low by removing a part of gas density enhancement associated with the surroundings of excised galaxies and AGN.

Structure function
For the sake of completeness, we also computed the third-order structure function of the velocity and density for the same reference galaxy cluster used above for the power spectra analysis, as usually this is more commonly done in study of other turbulent astrophysical environments (e.g. ?). Fig. 7 shows the third-order structure function of the velocity component tangential (upper panel) and longitudinal (central panel) to the separation, l, and of the particle density distribution (lower panel). The analysis was carried out on 800 kpc wide boxes, located at multiple distances from the centre of the cluster. We do not use any previously applied filtering method because the influence of the laminar flow is dominant mainly at the largest scales ( 1 Mpc). The maximum separation distance explored here is l ∼ 800kpc. We refer the reader to Vazza et al. (2011b), (2017) for further details on how the analysis was performed. Our analysis shows that, overall, both the gas density and velocity third-order structure functions have a Kolmogorov-like scaling, S 3 (l) ∝ l, within the separation range 100-800 kpc, similarly to the velocity power spectra. As in previous work (Vazza et al. 2011b), the steepening at the smallest separations is mostly driven by the increasing dissipation of the PPM scheme, which dampens fluctuations faster than in the Kolmogorov model. The normalisation of the velocity structure function increases from the outskirts to the centre of the cluster, with some variance related to the presence of substructures. The normalisation of the density structure function, on the other hand, decreases in a more steady way towards the  periphery, following the drop of the cluster gas density profile (see Fig. 1), which contributes to S 3 (l) with the density cube.
Recent studies investigated the structure function in galaxy clusters (???). Overall, the structure functions have a Kolmogorov slope (∝ l 1/3 for first-order functions), or steeper. In those studies the turbulence is driven by the black hole energy output and the observed scales range between 1 and 100 kpc, which cannot be probed by our simulations. On the other hand, our simulations are more suitable to capture the development of turbulence injected by mass accretion on ∼ 0.1 − 1 Mpc scales, and the complementary view of the structure function confirms that the dynamics of gas motions are reasonably close to a Kolmogorov model.

Testing our method in idealised simulations
Next, we validate our method against more idealised simulations. As previously discussed, due to the different nature of the simulations, it is not feasible to apply the very same algorithm that Mohapatra et al. 2020Mohapatra et al. , 2021 in their studies to our cluster sample. For instance, in Mohapatra et al. and Shi & Zhang, the Mach number and the Richardson number were parameters of the simulation. However, in cosmological clusters these quantities cannot be estimated a priori. Furthermore, the algorithm to evaluate the density fluctuations must be revisited as well, due to the presence of self-gravitating structures inherent in forming clusters, including the Itasca clusters (see Sec. 2.2). Therefore, we applied our analysis to data from the simulations by Shi & Zhang. In their work, they used the FLASH code (Fryxell et al. 2000) to simulate 3D local boxes with L x = L y = 200 kpc, L z = 250 kpc and a resolution of 256 × 256 × 320 cells. The aim of their study was to reproduce local ICM conditions in various regions for both relaxed and disturbed galaxy clusters. For this reason, they performed four simulation runs with increasing levels of stratification corresponding to different Froude numbers (Fr = √ 1/Ri). They used the Froude number to quantify the role of the buoyancy in each simulation. The initial velocity field was generated with the purpose to mimic a 3D Kolmogorov turbulence, with an injection scale L = 43.6 kpc. They also removed the compressive modes to generate merely solenoidal turbulence. Given the initial conditions, the evolution of the velocity distribution is fully dictated by the Mach and Froude number. In this sense, these simulations are quite similar to those of Mohapatra et al..
For comparison, we analysed the most stratified cluster data cube at 2 Gyr and estimated all the relevant quantities, such as the Richardson number and the logarithmic density fluctuations. In table 1 we compare our results with the parameters of the snapshot (S&Z method), both analysing the whole cube (full data cube) and boxes with 80 3 cells located in various spot of the cube itself. To estimate the Froude number we used the injection scale of their simulations, which is L = 43.6 kpc. The slope is consistent with Kolmorov's theory in the range between 100 -800 kpc, while it has a steeper behaviour at the smallest scales (l < 100 kpc). The normalisation increases with cluster-centric distance, especially for the density structure function due to the stratification of the ICM.
After filtering for density fluctuations, we obtained the logarithmic perturbations both performing a fixed-scale filtering Fig. 8. A comparison between the logarithmic density fluctuations estimated with a fixed scale filtering method and calculating the density mean in slices at a given z coordinate. The red line is the bisector, which highlights the one-to-one relation. method, with a scale L ≈ 43.6 kpc and computing the mean density on a slice at a given z coordinate, which is the method that Mohapatra et al. 2020Mohapatra et al. , 2021 used in their studies. We plot the results in Fig. 8 for a comparison between the two filtering algorithms. The black points are the results of the analysis performed in a box with 80 3 cells. Finally, we estimated the mean logarithmic density fluctuation considering the full data cube, obtaining σ 2 s = 0.048. This result is in close agreement with Fig. 8 of Mohapatra et al. (2021), which predicts σ 2 s ≈ 0.04 − 0.05 for a Froude number Fr ≈ 0.3 − 0.4. In summary, our results match reasonably well those in the literature, which instills confidence in our filtering methods, both for idealised as well as cosmological simulations.

Results
First, we focus on the relation between the density and velocity fluctuations in our cluster sample. Then, we investigate how the stratification can affect this relation and test the dependence of the density fluctuation on the Richardson number. A large part of the analysis was performed within the virial radius of the galaxy clusters in order to compare to observations.

The σ ρ − σ v relation
In observations, turbulence in galaxy clusters has been investigated using the relation between the gas density and velocity fluctuations. This has been followed up by numerical work in, both, cosmological and idealised simulations (Zhuravleva et al. 2014b, Gaspari et al. 2014). The goal of those studies was to derive the amplitude of gas turbulent motions solely from the analysis of surface brightness fluctuations. Those studies found a one-to-one scaling relation suggesting that such a relation is valid in, both, perturbed and relaxed clusters. Fig. 9 shows the δρ/ρ − δv/c s relation in our sample for, both, relaxed and disturbed clusters. The classification of the dynamical state is based on the centroid shift parameter, w, which measures the offset between the centre of mass and the position of the gas density peak of the cluster (Rasia et al. 2013). The centroid shift of each cluster, w i , was estimated averaging the parameter values of multiple lines of sight. Galaxy clusters are then classified comparing w i with the median of the centroid shift of the sample, which is w c = 0.0067 (which is half of the value used in Cassano et al. 2011), with relaxed clusters having w ≤ w c . The relation was obtained within the virial volume by applying a fixed scale filtering method with scale L = 300 kpc and averaging the density and velocity fluctuations over 600 kpc wide boxes (see Sec. 2.2). The outcome is a linear relation that depends on the dynamical state. Relaxed clusters have a slope m = 1.15 ± 0.06, where σ v ∝ mσ ρ , consistent with Zhuravleva et al. (2014b). Contrary to expectations, disturbed clusters show a flatter relation, with m = 0.40 ± 0.05.
In order to test the goodness of the fit, we used the Pearson's correlation coefficient wherex andȳ are the means of the variables used for the fit. R is a measure of linear correlation between two sets of data, but it ignores other types of relationship. The correlation coefficient ranges between -1 and 1. These limits represent a perfect (anti)correlation, while a value of 0 implies that the variables are uncorrelated. For this fit, we obtained R = 0.43 and R = 0.82 for perturbed and relaxed clusters, respectively. Moreover, in Fig. 11 we show the distribution of the ratio of the rms values of density and velocity fluctuations which depends depends slightly on the cluster-centric radius, regardless of the dynamical state of the cluster. Blue and black lines (or points) are related to relaxed and disturbed clusters, respectively, while the red line represents the relation found by Zhuravleva et al. (2014b). The Itasca simulation are non-radiative. This implies that some features of cool-core clusters, such as the peak (drop) in the density (temperature) radial profile in the centre of galaxy cluster, are missing. For this reason, it is not easy to make a proper classification, dividing our clusters in relaxed or disturbed, in this case. Due to this uncertainty, we performed the fit of the σ ρ − σ v relation also considering the total cluster sample (Fig. 10). We found a flatter relation, compared to a one-to-one relation, with a slope m = 0.63 ± 0.04 and a Pearson's correlation coefficient R = 0.60. Finally, we noticed that some fluctuations that are not driven by turbulence might affect the σ ρ − σ v relation, even after filtering. As a result, some of the σ v and σ ρ values in Fig. 9 might be caused by such contamination. Angelinelli et al. (2021) recently studied the radial profile of the clumping factor and the number density of clump in the Itasca sample. Here, we performed a similar analysis. We divided the inner region (within 3R 500 ) into six bins. For each of them we calculated the number of clumps as a function of the respective volume. The mass threshold for the clumps is 10 8 M and we set their maximum size to be 300 kpc. Fig. 12 shows the radial profile of the median clump number density. The shadowed regions represents the 16th and 84th clump number density distribution percentile boundaries. The number of clumps in the innermost region is higher in the perturbed clusters, making it less straightforward to study the relation between density and velocity fluctuations. However, the presence of selfgravitating structures is a problem for X-ray observations, too. In observations the identification of clumps is even more complicated, due to their low surface brightness. A solution might be to perform the analysis using a smaller filtering scale, in order to Fig. 9. σ ρ − σ v relation in our sample, which shows a dependence on the dynamical state of the galaxy cluster. The classification is based on the centroid shift parameter. The threshold value is w c = 0.0067. The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b). filter out also the smallest clumps. Nevertheless, the driver of turbulence in this simulation is the radial accretion. Thus, a smaller filtering scale would not be physically motivated since the injection scales generated by such event are likely larger than 200 kpc (Vazza et al. 2018;Angelinelli et al. 2020).

The
Ri − σ 2 s correlation In the previous section we discussed the relation between density and velocity perturbations (Fig. 9). However, the stratification Article number, page 9 of 14 A&A proofs: manuscript no. main Fig. 10. σ ρ − σ v relation for the total cluster sample.The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b). and, consequently, the buoyancy make the interpretation of gas fluctuations less straightforward. Hence, we computed the mean of the Richardson number and the logarithmic density fluctuations over 600 kpc wide boxes using a filtering method with a fixed scale of L = 300 kpc (see Sec. 2.2). In Mohapatra et al. (2020), the authors found a tight relation between density fluctuations and the Richardson number (highlighted by the orange line in Fig. 13), which we cannot detect in our sample (black points in Fig. 13). In particular, our simulated clusters show higher density fluctuations, up to three orders of magnitude, for Ri < 1, and a lower limit on log 10 (σ 2 s ) ≈ −1.5, regardless of Ri. We have already presented in Sec. 4 our tests on the reliability of the density fluctuations computed in the Itasca cluster. Even though our distribution does not perfectly match that of Eckert et al. (2017), the deviation in the surface brightness fluctuations is not sufficient to obtain density fluctuations similar to those of Mohapatra et al. (2020). Moreover, we have also shown in Sec. 2.3 that when applied to a similar (idealised) simulation setup, our pipeline correctly recovers the expected density fluctuations and Froude or Richardson number. We are therefore led to the conclusion that the setup in Mohapatra et al. 2020 and 2021 is much more idealised than more realistic ICM as produced by the Itasca sample, both concerning the solenoidal forcing of turbulence and the absence of self-gravitating gas substructures and bulk motions associated with accretion. In our cluster sample, the accretion of clumps is ubiquitous and it produces both solenoidal and compressive turbulence. This results in a distribution of flow velocities with typical averaged (on 600 kpc wide boxes) Mach numbers of ≈ 0.4 − 0.5 and likely locally supersonic motions (see Fig. 14). Another crucial difference is the distribution of the logarithmic density fluctuations, which is found to be significantly different from a log-normal distribution in our work. In summary, in our simulations we find that density fluctuations do not depend on Ri. Therefore, estimating the density perturbations as the variance of the probability distribution function, and relating it to the physical properties of the the distribution of the ratio between density and velocity rms. The black points (line) correspond to relaxed clusters and blue points (line) to the perturbed ones. The classification into relaxed and perturbed clusters is based on the centroid shift parameter, whose discriminating value is w c = 0.0067. underlying turbulent flow, appears to be challenging in realistic situations, even without taking into account observational effects and sampling issues.

The Ri − β correlation
The stratification and buoyancy might also affect the velocity distribution of the plasma, suppressing the radial motion. As a result, in the buoyancy-dominated regime (Ri > 1), the gas dynamics should be dominated by tangential motions (see Brethouwer et al. 2007, Mohapatra et al. 2020. Therefore, we look for a relation between the Richardson number and the anisotropy parameter. Since for Ri < 1 the hypotheses of homogeneity and isotropy are still valid, we should expect a mean β ∼ 0. On the other hand, a mean β < 0 should arise when we move to buoyancy-dominated regimes. However, looking at Fig. 15, we see that the dynamics of the ICM is biased by radial motions despite the stratification. As already noticed in previous work on this sample (e.g. Vazza et al. 2018), this ubiquitous dynamical feature of the ICM is most likely a consequence of radial accretion flows that are important drivers of ICM turbulence. Finally, we have explored the sensitivity of our results on the choice of the filtering scale and box size, by considering larger and smaller filtering scales. We find that qualitatively the results still hold (see Appendix 1 for more details).
The important conclusion is that the Richardson number is not able to characterise the effect of ICM stratification for cosmologically evolving clusters. Since accretion occurs largely via the infall of gas clumps whose motion is predominantly radial (e.g. Angelinelli et al. 2021), a predominantly radial (turbulent) Fig. 12. Radial profile of the number density of clumps for both relaxed (blue) and perturbed (red) galaxy clusters. The shadowed region represent the 16th and 84th percentiles of the bin's number density distribution. velocity is found at most cluster-centric distances in all clusters of our sample. This radial turbulence is unlike the stratified turbulence expected for a Ri ≥ 1 flow.

Radial profile of the power spectra slope
The velocity power spectrum quantifies the strength of turbulent fluctuations on different scales, and the slope of the power spectrum is customarily compared to Kolmogorov's theory.
We compute the radial profile of slope of the density and velocity power spectra in the Itasca cluster sample, which is shown in Fig. 16. The error-bars represent the 1σ error, which were estimated as the standard error of the weighted mean σ = , where σ i are the 1σ errors of each power spectrum at a given distance from the centre of the cluster. The slope of   15. Ri−β relation in our cluster sample. The radial bias is an imprint that the dynamics is dominated by radial motions, driven by the accretion. In this sense, the Richardson number does not have a predictive power about the real dynamics of the ICM. The orange line is used to highlight the turbulence (left) and buoyancy (right) dominated regime, while the blue line distinguishes a radially (β > 0) from a tangentially (β < 0) bias. the velocity spectrum shows a constant radial profile and agrees with the Kolmogorov's model, indicated by the blue line. On the other hand, the density power spectra have steeper slopes, albeit with a bigger scatter. The slope of the velocity spectrum is flatter compared to Vazza et al. (2017), who used the same cluster sample. The reason for this difference is explained by the varying normalisation of the power spectrum, which decreases as we move to larger radii (see Fig. 4). While in previous work we measured a unique power spectrum across a large volume (hence mixing local spectra with a slightly decreasing normalisation going to the cluster outskirts), here we computed spectra in smaller boxes. On scales of ≤ 600 kpc the turbulence appears even more Kolmogorov-like, despite the presence of radial biases and a non-stationary forcing. At the same time, the larger scatter in the slope of the density power spectrum at all radii questions whether any simple relation can be robustly derived between the two quantities.

Discussion and conclusion
We have investigated whether the turbulent gas velocity fluctuations in the ICM can be constrained using projected X-ray surface brightness fluctuations. Previous work (Gaspari et al. 2014;Zhuravleva et al. 2014b;Schuecker et al. 2004;Sanders & Fabian 2012) have already shown that the pattern of surface brightness fluctuations can, under certain limiting assumptions, Fig. 16. Radial profile of the slope of density and velocity power spectra in our cluster sample. The velocity slope is in agreement with Kolmogorov's theory and shows a constant profile as well as the slope of the density spectra, which is steeper. The blue line traces Kolmogorov's prediction for the velocity slope of a 3D spectrum. be translated into density fluctuations. However, in a stratified ICM, the δρ − δv relation is affected by buoyancy effects. Mohapatra et al. (2020), (2021) further derived an accurate relation between the two, parameterised by the Richardson number and the Mach number of the flow, in the very simplified case of turbulence-in-a-box simulations. So far, such analyses have been limited to idealised simulations or to the innermost cluster regions. Here, we used the Itasca cluster sample, which is a set of galaxy clusters simulated with the cosmological code ENZO with a uniform comoving spatial resolution of about ≈ 20 kpc. In order to compare the results to observations, we limited our analysis of the ICM to within the virial radius. Fig. 4 shows that while the density fluctuations increase as we move to larger radii, the velocity power spectra remain fairly constant. The increasing clumping factor in the outer region of galaxy clusters is the likely driver for this trend. As a consequence, a more accurate algorithm to mask clumps is required in observational and numerical analyses. This will be particularly challenging for observations, where the masking small-scale and faint clumps will remain limited. We can summarise our results as follows: • Consistent with previous work (Zhuravleva et al. 2014b, Gaspari et al. 2014, we find that a linear relation exists between the standard deviation of the density contrast, σ ρ and the normalised velocity fluctuations, σ v , albeit with a significant scatter across the sample. However, the slope of this correlation does not generally agree with the findings of the aforementioned studies. For the whole Itasca sample (Fig. 10) it results σ v ≈ 0.6σ ρ , with smaller velocity fluctuations than expected by applying the Gaspari et al. (2014) relation.
When we separate the sample in relaxed and perturbed clusters, we find that relaxed objects show a steeper slope, consistent with 1, and again a large scatter. Perturbed clusters follow instead a flatter, weaker relation, σ v ≈ 0.4σ ρ , with a low correlation coefficient. Fig. 12 shows that this relation may be contaminated by unfiltered clumps. This demonstrates the problem in using this relation to infer the turbulence in realistic clusters since this could lead to incorrect estimates, for example, of the turbulent heating rate (Zhuravleva et al. 2014a) and the kinetic energy available to explain radio halos (P H ∝ ρσ 3 v , e.g. Eckert et al. 2017).
• Although buoyancy is the dominant process in the largest part of the cluster volume (Fig. 3), we did not find a strong relation between the logarithmic density fluctuations and the Richardson number. The results in Sec. 3.2 show that our fluctuations are much higher compared to those of Mohapatra et al. (2020), (2021), who performed an analogous analysis in more idealised simulations, and are independent of the Richardson number. We attribute the difference to the unavoidable presence of compressive turbulence, consistent with previous work (e.g. Miniati 2015;Vazza et al. 2017;?) . Furthermore, the turbulent anisotropy parameter shows the effects of substantial radial accretion (see Fig. 15), which appears to be the dominant driver for the gas fluctuations in our simulations (e.g. Angelinelli et al. 2021). • Turbulence in the ICM is usually assumed to be homogeneous and isotropic. We show that this may not be valid, even if the spectra are found to follow Kolmogorov's scaling. This is primarily caused by radial accretion that dominates the dynamics of the ICM (see Fig. 15). In cluster cores, jets from central radio galaxies may produce a similar effect (Brüggen & Vazza 2015).
Finally, we point out that our cosmological simulations are optimised for a uniform resolution across the cluster volume and cannot attain the high spatial resolution of observations in cluster cores. Therefore, our simulations do not extend to scales below ≤ 10 kpc probed by some X-ray observations of nearby clusters of galaxies.
5. Appendix 1: On the Ri − σ 2 s relation at different filtering scales In Sec. 2.2.2 we showed that the Richardson number is very dependent on the filtering scale. In order to validate the robustness of our results, we performed tests with a filtering scale of 100 and 800 kpc. Due to the resolution of our simulation, an analysis with a much smaller scale cannot be accomplished. Fig. A1 and A2 show the Ri − σ 2 s and Ri − β relation for L = 100 kpc and L = 800 kpc, respectively. In essence, the results are almost identical to those of Sec. 3.2, as we did not detect any strong relation between the variables. It looks there is a Ri − σ 2 s relation for L = 800 kpc. However, there is a lot a scatter and the fit would definitely be different from that of Mohapatra et al. since the density fluctuations are much higher. Furthermore, a filtering algorithm with such a large scale, it might not be a properly choice, since we are essentially removing the bulk motions on scale of typical size of the cluster virial radius, leaving out the laminar flows on smaller scales. Besides, the bottom panel of Fig. A2 shows that a radial bias is still present, confirming that the Richardson number is not the ideal proxy to quantify the effect of stratification on turbulence in galaxy clusters.  s and Ri − β plot for a filtering scale L = 800 kpc. The orange line is used to highlight the turbulence (left) and buoyancy (right) dominated regime, while the blue line distinguishes a radially (β > 0) from a tangentially (β < 0) bias.