© ESO 2022

1. 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⁷ − 10⁸ 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 2016). Those authors found a 1D velocity dispersion ∼ 160 km s⁻¹, 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 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. 2009, 2011a, 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, 2011a). As another alternative, a multi-scale iterative filtering approach has been implemented in, for example, Vazza et al. (2012), Angelinelli et al. (2020), Vallés-Pérez et al. (2021). 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:

(1)

where v_l is the characteristic turbulent velocity on a scale l, and N_BV is the Brunt-Väisälä frequency:

(2)

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 Sect. 2 we describe the simulated cluster sample and the numerical method used to estimate the relevant variables of a turbulent ICM; in Sect. 3 we discuss the outcomes of the analysis, looking also for a comparison with previous work; in Sect. 4 we summarise our findings, highlighting some observational implications.

2. Methods

2.1. 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, σ₈ = 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 non-gravitational 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., Schober et al. 2012; Donnert et al. 2018; Vazza et al. 2018b; Domínguez-Fernández et al. 2019). 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² − 10³. 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 (Vazza et al. 2018b). 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₀ = 110 h kpc⁻¹ ≈ 157 kpc and the dark matter (DM) mass resolution was m_DM = 8.96 ⋅ 10⁷ 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 55 h kpc⁻¹ ≈ 78.4 kpc and the DM mass resolution was m_DM = 1.12 ⋅ 10⁷ M_⊙. Finally, inside the central (L₀/10)³ ≈ 6.3 Mpc³ 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⁻¹ ≈ 20 kpc.

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¹⁴ M_⊙ < M₂₀₀ < 3.32 ⋅ 10¹⁴ M_⊙. For further information on the Itasca cluster sample, we refer to Vazza et al. (2017), Wittor et al. (2017), Angelinelli et al. (2020).

2.2. Filtering method and numerical analysis

In order to disentangle bulk from turbulent motions, we used a tailored filtering method. Vazza et al. (2012, 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. (2018a) 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:

(3)

(4)

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. Figure 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).

Fig. 1. Slice of the original density and velocity 3D distributions in one cluster (left panels). Right panels: turbulent density and velocity fluctuations maps, after the application of the fixed-scale filtering method.

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.

2.2.1. 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: and , 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. Figure 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, , and the velocity contrast δv/c_s. The PDFs are obtained from the fluctuations within the virial radius (defined as, R₂₀₀, i.e. the radius enclosing a 200 overdensity with respect to the cosmic critical density).

Fig. 2. Probability distribution function of the gas density contrast, logarithmic gas density fluctuations, and velocity fluctuations. The colours correspond to different clusters, at z = 0. All statistics refer to within the virial radius.

2.2.2. 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 Sect. 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 Sect. 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.

Fig. 3. Top panel: radial profile of the Richardson number in our clusters. Middle panel: radial profile of the Brunt-Väisala frequency for each cluster of the sample. Lower panel: radial profile of the Richardson number considering different filtering scales L, for a single cluster of our sample. The horizontal red line separates the turbulence-dominated (Ri < 1) and buoyancy-dominated (Ri > 1) regimes.

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 kinetic-energy 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 β,

(5)

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. (2018a) for a further analysis of the anisotropy parameter in the Itasca cluster sample.

2.2.3. 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. 2011a for a discussion):

(6)

where v_i(k) is the Fourier transform of each component

(7)

The gas density fluctuations, |ρ(k)²|, 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).

Fig. 4. Density and velocity power spectra in a simulated cluster from the Itasca sample, obtained at different radii. The vertical orange bars indicate the selected inertial subrange used for the analysis. The normalisation of density spectra has a radial trend, with larger fluctuations in the outer regions.

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. 2011b for an example). Considering that bremsstrahlung is the dominant radiative mechanism at the typical temperature of plasma in galaxy clusters (10⁷ − 10⁸K), we estimated the surface brightness by projecting S_X = n²T^1/2 (where n and T are the plasma density and temperature in each cell) along multiple lines of sight. Figure 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 fluctuations, 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):

(8)

Fig. 5. Projected X-ray surface brightness in a cluster of our sample before (top panel) and after (bottom panel) filtering out the clumps.

where P_2D(k) is the 2D power spectrum and k is the wavenumber in kpc⁻¹ 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 non-thermal 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. (2017) galaxy cluster sample, which is 0.9. We observe that our fluctuations are by a factor ∼2 larger.

Fig. 6. Distribution of the X-ray surface brightness fluctuations, , in the Itasca cluster sample.

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 X-ray photons. Real observations, however, are affected by the uneven effective spatial resolution, and by the fixed sensitivity to X-ray brightness fluctuations. Moreover, it is not straightforward to implement exactly the same procedure for the masking of point-like 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.

2.2.4. 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., Kritsuk et al. 2007). Figure 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 ∼ 800 kpc. We refer the reader to Vazza et al. (2011a, 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₃(l)∝l, within the separation range 100–800 kpc, similarly to the velocity power spectra. As in previous work (Vazza et al. 2011a), 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₃(l) with the density cube.

Fig. 7. Structure function of the tangential (upper panel), transverse (middle panel) and numerical density (lower panel). 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.

Recent studies investigated the structure function in galaxy clusters (Li et al. 2020; Wang et al. 2021; Mohapatra et al. 2021). 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.

2.3. 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. (2020, 2022) used in their studies to our cluster sample. For instance, in Mohapatra et al. (2020) and Shi & Zhang (2019), 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 Sect. 2.2). Therefore, we applied our analysis to data from the simulations by Shi & Zhang (2019). 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 (). 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. (2020).

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³ 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.

After filtering for density fluctuations, we obtained the logarithmic perturbations both performing a fixed-scale filtering 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. (2020, 2022) 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³ cells. Finally, we estimated the mean logarithmic density fluctuation considering the full data cube, obtaining . This result is in close agreement with Fig. 8 of Mohapatra et al. (2022), which predicts 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.

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.

3. 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.

3.1. 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.

Figure 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 Sect. 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.

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 wc = 0.0067. The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b).

In order to test the goodness of the fit, we used the Pearson’s correlation coefficient

(9)

where 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.

Fig. 10. σρ − σv relation for the total cluster sample. The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b).

Fig. 11. Top panel: ratio of the amplitudes of density and velocity fluctuations as function of the distance from the cluster centre. Bottom panel: 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 wc = 0.0067.

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₅₀₀) 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⁸ M_⊙ and we set their maximum size to be 300 kpc. Figure 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 self-gravitating 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 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. 2018a; Angelinelli et al. 2020).

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.

3.2. The Ri −  correlation

In the previous section we discussed the relation between density and velocity perturbations (Fig. 9). However, the stratification 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 Sect. 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 , regardless of Ri.

Fig. 13. plot. Density perturbations are independent of the Richardson number and much higher than those found by Mohapatra et al. (2020) (orange line). The analysis was performed on the total cluster sample.

We have already presented in Sect. 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 Sect. 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, 2022) 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 underlying turbulent flow, appears to be challenging in realistic situations, even without taking into account observational effects and sampling issues.

Fig. 14. Distribution of the Mach numbers in our cluster sample, using averaging over 600 kpc. The tail of the distribution at high Mach numbers is caused by local, supersonic turbulence.

3.3. 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, 2022). 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. 2018a), 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 A for more details).

Fig. 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 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) 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.

3.4. 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 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.

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.

4. 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, be translated into density fluctuations. However, in a stratified ICM, the δρ − δv relation is affected by buoyancy effects. Mohapatra et al. (2020, 2022) 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.

Figure 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.

Figure 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 (, 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 Sect. 3.2 show that our fluctuations are much higher compared to those of Mohapatra et al. (2020, 2022), 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; Vallés-Pérez et al. 2021). 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.

Acknowledgments

The authors thank Xun Shi for providing data for Sect. 2.3. The cosmological simulations described in this work were performed using the ENZO code (http://enzo-project.org), which is the product of a collaborative effort of scientists at many universities and national laboratories. F. V. acknowledges financial support from the European Union’s Horizon 2020 program under the ERC Starting Grant ‘MAGCOW’, no. 714196. Work by T.W.J. was supported by the US National Science Foundation grant AST1714205. MB acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306, Germany’s Excellence Strategy – EXC-2094 “Origins” – 390783311.

References

Appendix A: On the 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. A.1 and A.2 show the 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 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. A.2 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.

Fig. A.1. and Ri − β plot for a filtering scale L = 100 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.

Fig. A.2. 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.

All Tables

All Figures

	Fig. 1. Slice of the original density and velocity 3D distributions in one cluster (left panels). Right panels: turbulent density and velocity fluctuations maps, after the application of the fixed-scale filtering method.
In the text

	Fig. 2. Probability distribution function of the gas density contrast, logarithmic gas density fluctuations, and velocity fluctuations. The colours correspond to different clusters, at z = 0. All statistics refer to within the virial radius.
In the text

Fig. 3. Top panel: radial profile of the Richardson number in our clusters. Middle panel: radial profile of the Brunt-Väisala frequency for each cluster of the sample. Lower panel: radial profile of the Richardson number considering different filtering scales L, for a single cluster of our sample. The horizontal red line separates the turbulence-dominated (Ri < 1) and buoyancy-dominated (Ri > 1) regimes.

In the text

	Fig. 4. Density and velocity power spectra in a simulated cluster from the Itasca sample, obtained at different radii. The vertical orange bars indicate the selected inertial subrange used for the analysis. The normalisation of density spectra has a radial trend, with larger fluctuations in the outer regions.
In the text

	Fig. 5. Projected X-ray surface brightness in a cluster of our sample before (top panel) and after (bottom panel) filtering out the clumps.
In the text

	Fig. 6. Distribution of the X-ray surface brightness fluctuations, , in the Itasca cluster sample.
In the text

Fig. 7. Structure function of the tangential (upper panel), transverse (middle panel) and numerical density (lower panel). 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.

In the text

	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.
In the text

	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 wc = 0.0067. The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b).
In the text

	Fig. 10. σρ − σv relation for the total cluster sample. The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b).
In the text

Fig. 11. Top panel: ratio of the amplitudes of density and velocity fluctuations as function of the distance from the cluster centre. Bottom panel: 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 wc = 0.0067.

In the text

	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.
In the text

	Fig. 13. plot. Density perturbations are independent of the Richardson number and much higher than those found by Mohapatra et al. (2020) (orange line). The analysis was performed on the total cluster sample.
In the text

	Fig. 14. Distribution of the Mach numbers in our cluster sample, using averaging over 600 kpc. The tail of the distribution at high Mach numbers is caused by local, supersonic turbulence.
In the text

Fig. 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.

In the text

	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.
In the text

	Fig. A.1. and Ri − β plot for a filtering scale L = 100 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.
In the text

	Fig. A.2. 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.
In the text