Merger-driven multi-scale ICM density perturbations: testing cosmological simulations and constraining plasma physics

The hot intracluster medium (ICM) provides a unique laboratory to test multi-scale physics in numerical simulations and probe plasma physics. Utilizing archival Chandra observations, we measure density fluctuations in the ICM in a sample of 80 nearby (z<1) galaxy clusters and infer scale-dependent velocities within regions affected by mergers (r<R2500c), excluding cool-cores. Systematic uncertainties (e.g., substructures, cluster asymmetries) are carefully explored to ensure robust measurements within the bulk ICM. We find typical velocities ~220 (300) km/s in relaxed (unrelaxed) clusters, which translate to non-thermal pressure fractions ~4 (8) per cent, and clumping factors ~1.03 (1.06). We show that density fluctuation amplitudes could distinguish relaxed from unrelaxed clusters in these regions. Comparison with density fluctuations in cosmological simulations shows good agreement in merging clusters. Simulations underpredict the amplitude of fluctuations in relaxed clusters on length scales<0.75 R2500c, suggesting these systems are most sensitive to missing physics in the simulations. In clusters hosting radio halos, we examine correlations between gas velocities, turbulent dissipation rate, and radio emission strength/efficiency to test turbulent re-acceleration of cosmic ray electrons. We measure a weak correlation, driven by a few outlier clusters, in contrast to some previous studies. Finally, we present upper limits on effective viscosity in the bulk ICM of 16 clusters, showing it is systematically suppressed by at least a factor of 8, and the suppression is a general property of the ICM. Confirmation of our results with direct velocity measurements will be possible soon with XRISM.


INTRODUCTION
Galaxy clusters are the most massive gravitationally-bound systems in the universe and continue to grow via hierarchical mergers and continuous accretion of matter from the intergalactic medium.These accretion processes are the dominant drivers of gas dynamics within the hot ( ∼ 10 7.5 K) intracluster medium (ICM) outside the innermost regions (e.g., Walker et al. 2019;Heinrich et al. 2021;ZuHone & Su 2022) that are often affected by feedback from central supermassive black holes (see, e.g., Hlavacek-Larrondo et al. 2022, for a recent review).Probing turbulent and bulk motions within the ICM is crucial for understanding energy partition during large-scale structure evolution, astrophysical processes that drive the evolution of galaxies within the ICM, and plasma physics, yet direct measurements of these motions remain challenging (for exceptions, see, e.g., Gatuzz et al. 2023;Sanders et al. 2020; Hitomi Collaboration et al. ★ E-mail: amheinrich@uchicago.edu 2018; Sanders & Fabian 2013).Directly resolving bulk and turbulent velocities of the hot gas in galaxy clusters, groups, and massive galaxies will soon be possible with the recently-launched XRISM satellite (XRISM Science Team 2020).
While waiting for data from X-ray microcalorimeters, such as XRISM, NewAthena (Barret et al. 2020) and LEM (Kraft et al. 2022), to become available, indirect techniques have been used to constrain ICM velocities, including resonant scattering (e.g., Hitomi Collaboration et al. 2017;Ogorzalek et al. 2017), inference from metal abundance distributions (Rebusco et al. 2005), shock widths (Nulsen et al. 2013), and X-ray surface brightness or Sunyaev-Zel'dovich fluctuations (e.g., Schuecker et al. 2004;Churazov et al. 2012;Zhuravleva et al. 2015;Walker et al. 2015;Romero et al. 2023).This last method is particularly attractive as it provides characteristic velocity amplitudes as a function of spatial scale (or wavenumber) -information that is difficult to extract, even from direct velocity measurements through Doppler shifts.
The X-ray surface brightness method utilizes a relation between the scale-dependent amplitude (related to a power spectrum) of gas density fluctuations, (/)  , and the one-dimensional characteristic velocity,  1D,k .This method has previously been calibrated using relaxed clusters from non-radiative cosmological simulations (Zhuravleva et al. 2014b) and idealized Coma-like cluster simulations (Gaspari et al. 2014).It was first applied to a sample of the brightest cores of relaxed clusters to measure velocity power spectra (e.g., Zhuravleva et al. 2014aZhuravleva et al. , 2018)).Recently, this relationship was verified for unrelaxed objects in cosmological simulations (Simonte et al. 2022;Zhuravleva et al. 2023), showing that the correlation is still valid, though the scatter is larger by a factor of ∼ 2 compared to relaxed systems (Zhuravleva et al. 2023).High-resolution idealized simulations of stratified turbulence (Mohapatra et al. 2020(Mohapatra et al. , 2021) ) showed that the relationship is more complicated and depends on the level of stratification, turbulence driving, and pressure/entropy scale height of the system, which could explain the scatter in the simplified relation.
Using mostly imaging data from XMM-Newton, the technique has been applied to unrelaxed objects and regions outside of cool-cores (Eckert et al. 2017;Zhang et al. 2022), where gas dynamics are driven mainly by cluster mergers.The relatively large PSF of that telescope, however, limits measurements to rather large scales.Given highresolution images of galaxy clusters obtained with Chandra, this method can provide velocity measurements on scales small enough to test plasma effects and small-scale physics in simulations.Luckily, after nearly 25 years of operations, the Chandra archive contains enough data to measure these velocities in a large sample of clusters.
Velocity measurements outside of cool-cores open up a variety of phenomena that can be studied.For example, turbulence in the ICM is theorized to be linked to diffuse synchrotron emission from giant radio halos, mostly found in dynamically disturbed clusters (e.g., Cassano et al. 2010).The GeV cosmic-ray electrons (CRe) powering the radio emission in halos are likely associated with 2ndorder Fermi acceleration (e.g., Brunetti & Lazarian 2007;Petrosian & Bykov 2008;van Weeren et al. 2019) of a seed CRe population.Measurements of turbulent velocities and their correlation with the radio emission strength could, therefore, provide valuable insights into the CRe acceleration mechanisms.
Additionally, the ICM is threaded by a (likely) tangled magnetic field with  ∼ a few .Although the magnetic pressure is small compared to the thermal pressure, the field can significantly alter the ICM microphysics, including the ICM's transport properties (e.g., Schekochihin & Cowley 2006;Kunz et al. 2022;Arzamasskiy et al. 2023).Multi-scale velocity measurements to scales close to the Coulomb mean-free-path of particle interactions ( mfp ) are crucial for understanding the ICM as a weakly-collisional high- plasma.Since  mfp typically increases with the radius in the ICM, probing microphysical scales with observations is most promising in regions outside dense cluster cores.
Finally, the partition of energy during matter accretion onto virialized halos and cluster mergers also depends on the poorly known microphysics of the hot plasma in the ICM.Measuring gas motions outside cluster cores, one can directly probe what fraction of gravitational energy converts into thermal and kinetic energies.The latter contribution also biases cluster mass measurements based on the assumption of hydrostatic equilibrium (e.g., Lau et al. 2009;Nelson et al. 2014a;Vazza et al. 2018).Measuring the non-thermal pressure and demonstrating that it is responsible for the mass bias are important for a few percent-level precision cluster cosmology with X-ray and SZ techniques.
In this work, we will, for the first time, measure ICM gas velocities across a range of scales using the density fluctuations technique in a large sample of galaxy clusters within a cosmologically-interesting region of radius of  2500c using archival Chandra data.The high spatial resolution of Chandra allows us to measure velocities on smaller scales (and a broader range of scales) than studied previously.We will probe gas dynamics driven by mergers and large-scale structure evolution, minimizing the effects of AGN feedback and gas cooling.We assume a flat ΛCDM cosmology with ℎ = 0.7 and Ω m = 0.3. c is defined as the cluster-centric radius within which the average density is  times the critical density of the universe at the cluster's redshift.
In Section 2, we describe our sample of clusters, image processing techniques, and power spectral analysis methods.In Section 3, we present individual and averaged measurements of (/)  , the characteristic ICM gas velocity  1D,k , and place limits on the ICM clumping factor and scale-independent characteristic velocity.Section 4.1.1discusses these constraints in relation to cosmological simulations.Finally, Section 4.2 focuses on relevant plasma physics, namely, testing the turbulent re-acceleration mechanism of CRe and constraining transport properties in the bulk ICM.

SAMPLE SELECTION AND DATA ANALYSIS
We aim to measure the power spectra of hot gas density fluctuations down to the smallest possible scales in a large sample of galaxy clusters.The Chandra X-ray Observatory is best-suited for our goal as it offers the highest angular resolution.We select systems where archival Chandra observations cover (nearly) the entirety of the cluster within a projected radius of  2500c from the cluster center, as determined from X-ray surface brightness contours at the same radius.The clusters are also brighter than the background level by at least a factor of two in this region.To have a sufficient number of photons in the X-ray images, we further require the total exposure time of each cluster to be greater than 100 ks.However, we also examine a few bright clusters below this threshold.Our final sample totals 80 clusters, with Chandra ObsIDs and cleaned exposure times presented in Table B1.We collected measured total masses ( 2500c ) and characteristic radii ( 2500c ) from the literature1 , summarizing their values and corresponding references in Table B2. Figure 1 summarizes our sample.Most clusters in our sample have masses ( 2500c ) between ∼ 8 × 10 13 and ∼ 10 15  ⊙ , and have redshifts between ∼ 0.05 and ∼ 0.5.

Initial data processing
The data for each cluster are reprocessed using the standard methods described in Vikhlinin et al. (2005), implementing scripts from ciao version 4.14.Event lists are reprocessed using CalDB version 4.9.7.Background flares are removed by filtering lightcurves using the lc_sigma_clip script.Non-X-Ray background files are created using combined products of the blanksky and readout_bkg routines.Event lists are restricted to the 0.5-3.5 keV band, where X-ray emissivity is almost independent of temperature for typical temperatures of galaxy clusters (Forman et al. 2007, see their Figure 1).The resultant event lists are combined using the merge_obs script, creating mosaic images of X-ray counts, background, a weighted instrument exposure map, and an exposure-weighted point-spreadfunction (PSF) map.A selection of exposure-corrected images of galaxy clusters from our sample is presented in Figure 2.

Classification of galaxy clusters
Based on these X-ray images and a literature review, we classify clusters in our sample into "relaxed", "intermediate", and "unrelaxed" systems.Relaxed clusters have cool-cores 2 , exhibit regular spherical or elliptical morphologies in X-rays, and show no signs of recent or ongoing mergers.Most intermediate clusters maintained their coolcores and showed evidence of minor mergers (e.g., gas sloshing).Clusters that are undergoing a merger but show regular morphology (e.g., A1689, which has spherical X-ray morphology and a weak cool-core, but hosts a radio halo (Vacca et al. 2011) and is likely undergoing a merger along the line-of-sight based on the galaxy redshift distribution (Girardi et al. 1997)) are also included into the "intermediate" category.Unrelaxed clusters are generally non-cool-core and show clear signs of ongoing major mergers, including shock fronts and "bullet-like" morphologies.Representative examples of these three categories are shown in Figure 2.While the divisions between relaxed and intermediate clusters, and between intermediate and unrelaxed clusters, are somewhat uncertain, the relaxed and unrelaxed subsamples are likely well separated in terms of merger effects.The "relaxed", "intermediate", and "unrelaxed" subsamples include 24, 30, and 26 clusters, respectively.

Deprojection analysis
To obtain radial profiles of gas density, temperature, and abundance of heavy elements for a subsample of galaxy clusters discussed in 2 Here, cool-core regions are defined as central regions with cooling times less than the Hubble time.
Section 4.2.2, we extract projected spectra in a set of radial annuli.These spectra are deprojected using a flavor of the onion-peeling technique described by Churazov et al. (2003).The resulting deprojected spectra are fit in the broad 0.6 − 8 keV band using XSPEC code (Arnaud 1996) and APEC plasma model (Smith et al. 2001) based on AtomDB (Foster et al. 2012) version 3.0.9.The abundance of heavy elements is fixed at 0.5 relative to Solar abundances of heavy elements taken from Anders & Grevesse (1989) in the fitting.However, we verify that treating the abundance as a free parameter does not affect the measured density and temperature profiles.

Residual images
Point sources in the X-ray images are identified via the wavdetect tool (accounting for the PSF), visually confirmed, and masked.The exposure-corrected and background-subtracted image is then radially binned into logarithmically-spaced annuli to create a surface brightness profile, (), as a function of the projected radius , for each cluster.We check that X-ray emission from each cluster is observable above the background within the region of interest.We then perform a least-squares fit of this profile to a spherically symmetric double -model: where the free parameters are the surface brightness normalization   , core radius  c, , and the   parameter, and  is the number of components (one or two −models).This fitted -model is used as the unperturbed surface brightness model for each cluster.For unrelaxed clusters that do not have the excess central emission associated with a cooling flow, the second -model component,  2 , is fixed to zero.To obtain the images of X-ray surface brightness fluctuations (residual images), we divide the initial cluster images by their best-fitting models.
Many clusters exhibit an elliptical X-ray morphology reflecting the underlying gravitational potential of the dark matter halo (e.g., Umetsu et al. 2018;Gonzalez et al. 2021).We correct our fitted -model to better reflect this ellipticity by performing a Levenberg-Marquardt two-dimensional fit of the exposure-corrected image to an elliptical -model, where the projected radius  in equation ( 1) has been replaced with the elliptical radius  e , with additional free parameters  (ellipticity) and  (orientation).The ratio of the bestfitting elliptical and spherical models is then used to correct the spherically symmetric -model described above.
Additionally, the gas density distribution in dynamically disturbed clusters is often asymmetric.Assuming a spherically-or even elliptically-symmetric underlying model in such systems can significantly bias the measurement of density fluctuations on large scales.To model an asymmetric cluster atmosphere, we "patch" the -model on large spatial scales following Zhuravleva et al. (2015).For details, see Appendix A. In short, the residual image is smoothed by a Gaussian with a large width to obtain a correction for the large-scale deviations.The unperturbed model (normalized by the image mask smoothed on the same scale) is then multiplied by this correction to get a patched model.This "patched" model effectively removes the large-scale asymmetries in residual images.We caution that the choice of patching scale is non-obvious.In our work, we calculate the power spectrum assuming several patching scales to determine which scale visually removes the asymmetry from the residual image.Our analysis only includes measurements of (/)  at scales that are unaffected by the choice of patching scale.images are shown, lightly smoothed with a Gaussian for display purposes.Solid white circles represent  2500c , while dashed circles mark the cool-core radius if the cluster has one.A1835 is a relaxed, cool-core cluster that does not show any signs of mergers (Schmidt et al. 2001;Govoni et al. 2009).A401 is a cool-core cluster, however, it is undergoing a merger as evidenced by the sloshing cold front to the south (white arrow), the proximity of the cluster to A399, and the presence of a giant radio halo (Murgia et al. 2010).For these reasons, we place it in the intermediate category.A2034 is a clear example of an unrelaxed, merging cluster as it lacks a cool-core, shows irregular morphology, contains a merger shock north of the cluster center (white arrow Owers et al. 2013), and hosts a radio halo (Botteon et al. 2022).

Power spectral analysis
To measure the power spectra of density fluctuations from the observed residual images, we follow the analysis steps described by Churazov et al. (2012).We employ the Δ-variance method (Arevalo et al. 2012) to measure the power spectrum of surface brightness fluctuations,  2D ()3 , in the residual image.The method is tuned for the non-periodic data with gaps and is best suited for power spectra that are smooth functions.In brief, given a masked image  of the background-subtracted X-ray counts and the exposure map  (product of the Chandra exposure map, unperturbed density model, and mask ), the power spectrum ( 2D ()) of the residual image at a wavenumber  is given by where with Gaussian kernels   1 and   2 of widths 0.225 , where  ≪ 1.We can deproject the 2D power spectrum of surface brightness fluctuations,  2D (), to a 3D power spectrum of X-ray emissivity fluctuations,  3D (), and obtain the amplitude of density fluctuations at each wavenumber, i.e., (/)  = √︁ 4 3D () 3 .Assuming that density fluctuations are a homogeneous and isotropic random field,  2D is related to  3D as where  ( z ) is the Fourier transform of the normalized emissivity distribution along the line of sight.The emissivity, and therefore the deprojection factor (integral in equation ( 4)), depends on the projected radius .We calculate this deprojection factor at each line of sight and include it in the exposure map from equation (3).We measure (/)  at each wavenumber  and relate it to the characteristic velocity of subsonic gas motions through the statistical relationship where  1D,k is a scale-dependent characteristic one-dimensional velocity of gas motions in the ICM,  s is the sound speed and  is a proportionality coefficient.Recently updated analysis of gas fluctuations in galaxy clusters from cosmological numerical simulations (Zhuravleva et al. 2023, hereafter Z23) confirmed a strong correlation between velocities and density fluctuations, finding that  varies slightly depending on the cluster dynamical state.Note that this result differs from the conclusions of Simonte et al. (2022), who filter the velocity field in order to separate bulk and turbulent motions.Given the difficulty in unambiguously defining the filtering scale (both in observations and simulations), we conservatively measure the whole velocity field.Following Z23, who average  between scales of 60 to 300 kpc, our adopted values for  are summarized in Table 1.
In the near future, this relationship will be tested using direct measurements of the turbulent velocity dispersion with XRISM.The first such measurement performed by Hitomi was broadly consistent with  ≈ 1 in the Perseus cluster core (Hitomi Collaboration et al. 2018).
As our goal is to probe merger-driven gas motions and density fluctuations, we derive the density fluctuation power spectrum within a cluster-centric projected radius of  2500c , where the effect of the non-X-ray-background is minimized.We omit (mask) cool-cores where appropriate, as these are often perturbed by AGN activity (e.g., Heinrich et al. 2021;Walker et al. 2018) and strong cooling (Fabian 1994;Hudson et al. 2010).We define  cool as the radius within which the gas cooling time  cool is less than the age of the universe  age ().

Bright substructures
Galaxy clusters continuously grow through mergers and often exhibit bright substructures that can dominate the signal if not properly excluded.Our primary goal is to characterize volume filling (bulk) gas perturbations.As such, for each cluster with prominent substructures, we produce two power spectra: one with substructures included and one without.Each cluster with identifiable substructures is tested to determine the region size that sufficiently masks each feature.Figure 3 provides an example of this test on the Bullet cluster.This cluster is undergoing a major merger in the plane of the sky that is responsible for the bright cold front (the eponymous "bullet"), the remnant of the infalling system's cool-core, and the sharp surface brightness discontinuity to the west of the "bullet", a bow shock (Markevitch et al. 2002).We first calculate (/)  for the entire cluster within  2500c , shown in Figure 3 as the orange spectrum.The blue and green spectra are then calculated after masking their respective regions shown on the inset.One can see the spectrum significantly decreases when the blue region is masked, but shows little change in the case of the larger green region.Therefore, we conclude the blue region is sufficient to remove the effect of these substructures on the spectrum as a whole.The slight difference between the blue and green spectra at large spatial scales is due to the lower sampling of these scales when the region size is increased.In our final analysis, we only include scales that are largely unaffected by the choice of masked region.

Limits on observable scales
Each cluster must be carefully examined to determine which length scales are reliably probed within  2500c , given the available data.On large scales, our measurement is limited by the region size.When 1/ approaches this radius, the number of large-scale modes decreases significantly, causing the spectrum to flatten and decrease around  < 1/ 2500 .This problem of "sample variance" is treated in detail by Dupourqué et al. (2023).In this work, we conservatively limit our measurements to scales below this peak, and caution that this minimum wavenumber does not necessarily represent the turbulent injection scale.
On small scales, the observations are limited by Poisson noise, unresolved point sources, and suppression of structures due to the Chandra PSF.Note that though the PSF is nominally quite small, Gaussian smoothing in the Δ−variance method makes the PSF effect noticeable at larger scales than the nominal PSF size.
The Poisson uncertainty associated with the number of X-ray counts per pixel contributes a scale-independent component to the measured  3D of order Σ raw /Σ 2 , summing over the raw X-ray counts image and the square of the exposure map.We can accurately model this contribution by creating many (in practice,  pois = 60) realizations of random images scaled by the Poisson uncertainty, the square root of the X-ray counts/pixel.The average power spectrum of these randomized images is subtracted from the nominally measured  2D or  3D and used to inform the 1 uncertainties.
While we detect and remove visible point sources from our data, there is likely a contribution of unresolved point sources to the power spectrum on small scales.Following Churazov et al. (2012), we attempt to quantify the power spectral contribution of these unresolved point sources.We first measure the power spectra of bright, identifiable point sources,  bright , by calculating the difference between power spectra when they are excised and when they are included.This  bright can then be rescaled to estimate the contribution of undetected point sources based on their flux distribution.In each cluster, we measure the background-subtracted flux of each point source, determining the number of sources per flux bin /.Between the brightest point source in the field ( max ) and the dimmest ( min ), this distribution follows a power law with respect to F. It is then reasonable to conclude that  bright is proportional to the total flux of these point sources, namely, By definition, unresolved point sources occupy fluxes between  = 0 and  =  min .Assuming that the slope of the flux distribution does not change in this region, we can find the relative contribution of unresolved point sources by numerically integrating The "flux ratio" between these two integrals is calculated for each cluster and multiplied by  bright to find  faint .We then limit our measurements to scales where  3D ≫  faint .Finally, we consider the smoothing of small-scale fluctuations by the Chandra PSF.As the Δ-variance method calculates a lowresolution, smoothed power spectrum, it is not sufficient to limit our measurements to scales greater than ∼ a few arcseconds.Instead, we create mock observations of many randomly placed point sources in the region of interest and measure two power spectra: one of the raw image where each point source occupies one pixel, and one of the point sources after they have been convolved with the exposureweighted Chandra PSF.The ratio of these,  PSF =  conv / raw , represents the fractional effect of the point spread function on our measured  3D .We calculate this correction several (10) times, each time generating a new map of randomly-distributed point sources, and take the average values of  PSF at each wavenumber.We divide the Poisson-subtracted  3D by this correction to obtain the suppressed, high- power, and limit our measurement to scales where 1 −  PSF < 25 per cent.
After all the corrections and checks of systematic uncertainties, we find a range of scales for each cluster where our power spectra measurements are most reliable.These scales are listed in Table B3 and shown in Figure 4.The median cluster in our final sample is observable between ∼ 120 and ∼ 360 kpc, spanning a range of  max / min ≈ 3. Most clusters ( = 68) are observable between scales of 200 to 300 kpc.In our sample, 24 clusters are measurable at scales less than 100 kpc.The smallest scales in the sample (1/ ≈ 60 kpc) are achieved in A3667 and A1367.

Density fluctuations and clumping factors
We measure the scale-dependent amplitude of density fluctuations (/)  and velocity power spectrum  1D,k in a total of 80 galaxy clusters across length scales ranging from 50 kpc ≲ 1/ ≲ 750 kpc.Measurements from the bulk ICM (substructures excluded) are shown in Figure 4, where the sample is divided into relaxed, intermediate, and unrelaxed clusters.In this figure, we plot (/)  and  1D,k against wavenumber/spatial scale in units of  2500c (bottom/top axis), where line color corresponds to the relative uncertainty  ( /)  /(/)  with 1 uncertainties derived from Poisson noise as described in Section 2.7.Weighted sample averages for each morphological category are additionally shown in grey, with the shaded region corresponding to the average ±1 standard deviation.The top panels show the number of clusters with measured (/)  at a given scale.Figure B2 shows equivalent measurements taken from the entire cluster, with substructures included in the calculation.We stress that the latter measurements do not represent fluctuations in the bulk ICM.
Typical amplitudes of density fluctuations in the bulk ICM are ∼ 12 ± 2 per cent for relaxed clusters, increasing to ∼ 19 ± 5 per cent for unrelaxed clusters when measured at 1/ = 0.5  2500c (Figure 4).If substructures are included in the measurement, the level of fluctuations in relaxed clusters remains about the same, while the typical amplitude increases to ∼ 40 per cent in unrelaxed clusters B2).One would expect to see increasing density fluctuations with merger activity, and while this is the case in the entire region (Figure B2), the effect is smaller once prominent substructures are removed.
The difference in amplitudes in the bulk ICM translates to a factor of ≳ 2.5 more ICM kinetic energy at length scales smaller than 0.5 2500c in unrelaxed clusters compared to relaxed ones.This suggests that fluctuations in the bulk ICM are not strongly linked to the current merger state of the cluster.This may be because there has not been sufficient time for fluctuations in clusters with ongoing mergers to cascade to these relatively small scales.Another possibility is that fluctuations in the bulk ICM of unrelaxed clusters dissipate quickly.The opposite may also be true, that those fluctuations present in relaxed clusters are "left over" from mergers in the cluster's more distant history.
It is additionally interesting to assess whether the density fluctuation amplitudes in relaxed vs. unrelaxed clusters belong to different populations.A Kolmogorov-Smirnov test performed at each wavenumber shows that it is difficult to distinguish between relaxed and intermediate clusters.However, density fluctuations measured in unrelaxed clusters (including substructures) are distinct from relaxed and intermediate clusters at the  < 0.05 level between scales of 1/ ≈ 0.2 to 0.8  2500c .Therefore, at these wavenumbers, the amplitude of density fluctuations could be used as an indicator of cluster dynamic state, complementing the commonly used indicators based on, e.g., the asymmetry of the X-ray surface brightness, centroid shifts, power ratios, and symmetry-peakiness alignment (e.g., Cerini et al. 2022;Yuan et al. 2022;Cui et al. 2017;Mantz et al. 2015b).Although there is some overlap between the intermediate and unrelaxed populations, (/)  > 20 per cent at 1/ = 0.5 2500c is typically indicative of a major merger.
The scale-dependent amplitude of density fluctuations (/)  can additionally be used to provide constraints on the gas clumping factor, , in the ICM through a simple relation (Zhuravleva et al. 2015): This expression should be integrated over the entire range of wavenumbers, however, our measurements are limited to a smaller range of scales described in Section 2.7 and listed in Table B3.Therefore, we integrate across this range of measurable scales, interpreting this result as a lower limit on , which we plot in pink in Figure 5.
To estimate the upper limit on , we fit (/)  to a power-law via an MCMC chain (Hogg & Foreman-Mackey 2018) using the emcee package (Foreman-Mackey et al. 2013), and integrate this power-law from  = 1/ 2500c to  = ∞.An example of a fitted power spectrum is shown in Figure B1.The choice of  = 1/ 2500c is rather arbitrary, however, reasonable given the size of the considered region.These upper limits are plotted as unfilled blue markers in Figure 5, with uncertainties derived from the 84th and 16th percentile of the MCMC chain.Both lower and upper limits on the clumping factor, as well as the best-fit parameters of this power-law, i.e., the normalization of (/)  at 1/ = 2500c and the slope , are presented in Table B3.
The weighted-average median upper limit on the clumping factor in our sample is  ≈ 1.05±0.04.In relaxed and intermediate clusters, the average is 1.04 ± 0.02, compared to 1.09 ± 0.04 in unrelaxed clusters.We find reasonably good constraints on  in 23 clusters, where the median lower limit on  − 1 is within a factor of 3 of the upper limit.Namely, these clusters are Cygnus A, A1795, A2597, the Bullet cluster, A2219, A3667, A2319, A2142, A2034, A1664, A1650, MACS0717, A401, A1835, A1689, A2256, RXC1504, A3112, A2626, A2204, 4C+37.11,A2029, and A2063.Of these, most of the median upper limits are within the 1.02 ≲  U ≲ 1.06 range.On the high end are the Bullet cluster ( U ≈ 1.12) and MACS0717 ( U ≈ 1.17), objects with (likely) the highest clumping factors in our sample.A1914 and RBS1316 also have high  L but have significantly larger uncertainties.We find the lowest median upper limit on the clumping factor ( U ≲ 1.02) in A1650, A2199, and A2029, despite the presence of sloshing in A1650 and A2029.In comparison with other observations, our measurements are consistent with those of Eckert et al. (2015), who found  ≈ 1.05 ± 0.04 within 0.5 500c in 31 clusters.We only confirm a single cluster with  outside of this range (MACS0717).

Velocity power spectra and non-thermal pressure
We convert our measured (/)  into velocity power spectra using equation ( 5).We estimate the median sound speed in the region of interest from either our measured temperature profiles or those found in the literature, shown in Table B2.The proportionality coefficient  has been recently updated by Z23 (see their Figure 7) using a sample of 78 clusters from non-radiative cosmological simulations.In that work, the amplitude of density fluctuations and the Mach number of gas motions were calculated with an equivalent methodology to the one used here within a similar region.They find average values of  = (/)  /M 1D,k between 1/ = 60 and 300 kpc, where the simulation resolution is not expected to affect the measurement.The  values we use are summarized in Table 1.Resulting velocity power spectra ( 1D,k ) are presented in Fig 4 .Additionally, velocities in each morphological category are averaged in the same fashion as /, plotted in grey regions with uncertainties as in the middle row.We see that, on average,  1D,k is correlated to cluster dynamical state.On a length scale of 0.5  2500c , the typical  1D is ∼ 150 km s −1 for relaxed clusters, ∼ 200 km s −1 for intermediate clusters, and ∼ 250 km s −1 for unrelaxed clusters.The velocity amplitude decreases with scale, and the slope is consistent with the Kolmogorov model within the uncertainties.
To estimate the characteristic (i.e., scale-independent) velocity associated with these gas motions, we use the integrated / and calculate the one-dimensional Mach number, M 1D = (/)/.This measurement can additionally be used to calculate the kinetic pressure fraction: where  kin is the kinetic gas pressure,  thm is the thermal pressure,  tot is the total pressure, and  1D =   M 1D is the characteristic one-dimensional velocity, assuming an adiabatic index of  = 5/3.As stated previously, we do not attempt to distinguish between bulk ( bulk ) and turbulent ( turb ) gas motions.To calculate  kin / tot , we assume that High-resolution X-ray spectroscopy from XRISM, NewAthena, and LEM (XRISM Science Team 2020; Barret et al. 2020;Kraft et al. 2022) will be able to distinguish these components to confirm our measurements and this relationship.
Figure 6 shows the derived values for M 1D and  kin / tot for all clusters in our sample (see also Table B3) with equivalent notation  8) from  min to  max , providing a lower limit on each value.We additionally fit the observed spectrum to a power-law via MCMC and integrate from  = 1/ 2500c to  = ∞.The median values, along with upper and lower uncertainties, are presented in blue with unfilled markers.These values could be interpreted as upper limits.Clumping factors measured in the Omega500 cosmological simulations (Nagai & Lau 2011;Zhuravleva et al. 2013) are plotted in shaded regions in their corresponding morphological categories.Grey regions show simulations where radiative processes (e.g., cooling and star formation/feedback) are included, while yellow regions show results from non-radiative simulations.
to Figure 5.We find an average best-fit Mach number in the range of 0.18 ± 0.05 across the entire sample, consistent with subsonic gas motions.The highest characteristic velocities we can confirm are found in MACS0717 ( 1D ≳ 575 km s −1 , M 1D ≳ 0.27) and RBS1316 ( 1D ≳ 435 km s −1 , M 1D ≳ 0.22).In all respects (lower limits and best-fit clumping factor/characteristic velocity), this makes RBS1316 (a.k.a RX J1347.5-1145Johnson et al. 2012) the most disturbed cool-core cluster, indicating that significant merger effects are taking place.
The lowest velocities are found in A2597, A1650, A2626, A2199, and A2029, with upper limits on  1D < 200 km s −1 , all of which are cool-core clusters.Of these, A1650, A2626, and A2029 show some evidence of mergers in the form of sloshing, while the rest are significantly relaxed.For non-cool-core clusters, the lowest velocities we confirm are A3667 and A2256, with upper limits below  1D = 300 km s −1 .Both clusters have significant substructures in the form of cold fronts (e.g., de Gasperin et al. 2022;Rajpurohit et al. 2023) and bow shocks in the case of A3667.Confirmation of the characteristic velocities in the bulk ICM with direct observation should carefully exclude emission from these prominent substructures.
Measuring non-thermal pressure support is important for cosmology based on hydrostatic cluster mass measurements (e.g., Allen et al. 2008;Vikhlinin et al. 2009;Mantz et al. 2016a;Clerc & Finoguenov 2023).We place upper limits (denoted with the lower index ) on the fraction of kinetic pressure in the bulk ICM at ( kin / tot )  < 10 per cent in a total of 22 clusters within the  2500c region, with ( kin / tot )  < 5 per cent in A2219, A1650, PKS0745, A1835, A2199, and A2029.A2029 has the lowest non-thermal pressure contribution of ∼ 1.9 per cent.Most clusters among these 22 are classified as relaxed or intermediate.Three clusters in this subsample, however, are unrelaxed clusters: A3667, A2034, and A2256, each with ( kin / tot )  < 8 per cent in the bulk ICM despite significant ongoing merger processes.For these clusters, it is important to note that the measured kinetic pressure in the bulk ICM is unlikely to reflect the total kinetic pressure in the entire cluster, contributions from merger-driven substructures could be significant.In our sample, the weighted-average best-fit upper limit for the kinetic pressure fraction (⟨ kin / tot ⟩  ) are approximately 6 ± 2 per cent for relaxed clusters, 4 ± 2 per cent for intermediate clusters, and 7 ± 4 per cent for unrelaxed clusters.Note that intermediate and unrelaxed clusters will likely have significant kinetic pressure associated with discrete substructures, which are removed from these measurements.
Recent constraints have been placed on the kinetic pressure fraction in a similar region of the Perseus (de Vries et al. 2023) and Ophiuchus (Gatuzz et al. 2023) clusters, finding 6 ≲  kin / tot ≲ 13 per cent and  kin / tot ≲ 18 per cent, respectively.All of our clusters agree with the latter measurement, and all but 12, with upper limits on  kin / tot < 6 per cent, agree with the former.work, the authors split the simulated clusters into relaxed, intermediate, and unrelaxed systems and considered fluctuations in the bulk ICM (i.e., removing prominent substructures) within a similar region ( ≤ 0.5 500 ).The authors considered fluctuations in simulations with non-radiative physics, which apply to the regions considered in our work (i.e., beyond cool-cores and least affected by cooling and feedback processes).
Figure 7 shows the average amplitude of density fluctuations, with one standard deviation as uncertainties, in simulated and observed clusters in pink and blue, respectively.In clusters with evidence of mergers (the intermediate and unrelaxed subsamples), observed and simulated amplitudes of density fluctuations are consistent across all scales.In relaxed clusters, however, the amplitudes diverge on small scales, reaching a factor ≳ 2 difference between the mean values at 1/ ≲ 0.3 2500c (typically ∼ 120 kpc).Performing a Kolmogorov-Smirnov test at each length scale, we confirm that the distribution of density fluctuations observed in our relaxed cluster sample are distinct from the simulations between 1/ = 0.75 and 0.3  2500c (∼ 300 − 120 kpc, the dark blue bracket in the top-left panel), with the probability of the two samples being derived from the same distribution being < 0.05 at these scales.On scales smaller than 1/ ∼ 0.3 2500c , we observe a similar tension, however, the number of clusters where we could probe these scales is small, less than 10.If the observational and simulated measurements are compared with wavenumber in kpc −1 (instead of re-scaling with  2500c ), we observe this tension in relaxed clusters between ∼120 and 300 kpc.
We checked whether observational biases could cause this tension.Fluctuations in five clusters in the relaxed sample (A2537, A3847, A478, MACS1621 and MACS2140, see Figure 4) are measured on large scales (1/ ≳ 0.6 2500c ) with (/)  < 10 per cent.Extrapolating these amplitudes to smaller scales gives very small amplitudes of fluctuations, which are difficult to measure reliably given Poisson noise, background, and contribution from unresolved point sources.However, if these spectra were observed on small scales, they could reduce the average amplitude of observed density fluctuations, bringing it closer to the simulation average.To examine this potential bias, we extrapolate the power spectra in these five clusters across smaller scales, fitting the observed spectra with a power-law as described in Section 3.1, and assuming the slope of (/)  extends across the entire range of length scales.With these extended power spectra, we reanalyze the tension in relaxed clusters.While the distributions of (/)  in the observed and simulated subsamples become compatible at small scales (where we initially observed less than 10 clusters), the tension remains strong between 0.3 ≲ 1/ ≲ 0.75  2500c .This conclusion remains if we assume a fixed Kolmogorov slope of −1/3, or if the extrapolation across all scales is applied to every cluster in the sample.
Our measured power spectra could also be affected by unresolved point sources or, e.g., additional variations on the detector, especially on small scales.To check the potential effect of these, we take a more conservative value for the maximum wavenumber measured reliably in each cluster by only considering wavenumbers where we do not observe any flattening of the (/)  slope.We find that even with reduced limits, the divergence between observation and simulations remains across most of the range.We additionally compare the observed spectrum when no correction for the Chandra PSF (see Section 2.7) is applied, and find no change in the length scales where the divergence is observed.
One could expect that density fluctuations also correlate positively with redshift since the merger rate increases with redshift according to hierarchical structure formation.Our sample covers a range of redshifts up to  ∼ 1, while the Omega500 clusters were analyzed at  = 0. To check whether the observed tension could be caused by this redshift evolution of /, we recalculate the Kolmogorov-Smirnov test excluding high-redshift systems, and conclude that the redshift evolution did not significantly affect the result.
Finally, we check if our morphological classification of galaxy clusters has any effect on the observed tension between the observations and simulations.We identify clusters that may be on the boundary between morphological categories and recalculate the averaged spectra and K-S test when the morphological classification is randomly varied between categories for each boundary object.We find that the tension is strong even if the classification is varied or if these boundary objects are removed from the sample entirely.To conclude, the tension appears real and is stable to observational systematics.Simply put, we do not observe clusters with the lowamplitudes fluctuations that are common in simulations at scales of 1/ ≲ 0.75 2500c .
Note, the displayed scales from simulations are not expected to be affected by the resolution of simulations, which becomes evident on length scales below ∼ 60 kpc (≲ 0.15 2500c , Z23, see their Section 4.2).Moreover, the divergence in averaged (/)  is only observed in relaxed clusters and not in the intermediate or unrelaxed subsamples.Therefore, this divergence in relaxed clusters is unlikely to be driven primarily by numerical viscosity.
As the factor of ∼ 1.5 to 2.5 difference between the amplitude of density fluctuations on length scales smaller than 0.75 2500c in relaxed clusters cannot be explained by known observational/instrumental effects and is unlikely to be caused by the resolution of cosmological simulations, it could be attributed to missing physics in the simulations we consider.Radiative processes (e.g., gas cooling), AGN feedback, galactic outflows, and magnetic fields could all increase (/)  at these scales.Fluctuations in relaxed clusters appear to be more sensitive to these physics compared to perturbed systems, in which ongoing mergers continuously "stir" the ICM and contribute significantly to the observed amplitude of density fluctuations.While it is straightforward to extend the analysis of density fluctuations to other simulations and simulation runs that include these additional physics, it is beyond the scope of this work.The observed tension also means that relaxed clusters should be selected for testing additional physics implemented in simulations, especially multi-scale physics, against observations.With these measurements, we cannot distinguish between which of the above processes is the primary driver of the observed divergence.The central AGN mostly affects gas dynamics within the cool-core and is not influential in the regions considered, however AGN and galaxy-scale outflows originating from other cluster galaxies could affect density fluctuations in relaxed clusters more strongly than in merging ones.Additionally, while cooling times in the regions considered are greater than the Hubble time, radiative physics could still affect (/)  in simulations (e.g., by creating clumpy structures, see Section 4.1.2below).
Finally, the fact that relaxed clusters appear "smoother" in simulations compared to observations could have consequences for cluster cosmology predictions based on hydrostatic mass measurements (e.g., Mantz et al. 2015a;Vikhlinin et al. 2009;Pratt et al. 2019).Namely, the predicted contribution of non-thermal pressure from simulations could be underestimated (e.g., Lau et al. 2012;Nelson et al. 2014b).We leave the quantitative exploration of this issue for future works.

Clumping factor and non-thermal pressure
It is interesting to compare the estimated ICM clumping factor shown in Figure 5 with predictions from numerical simulations.In earlier works based on Omega500 simulations, Nagai & Lau (2011) found that 1.03 <  < 1.04 in relaxed clusters and 1.16 <  < 1.26 in unrelaxed ones within  2500c when radiative processes (e.g., cooling, star formation/feedback, hereafter CSF) are included (see also Zhuravleva et al. 2013).The non-radiative runs predict a higher level of clumpiness, ∼ 1.14 (1.48) for relaxed (unrelaxed) systems, see gray and yellow regions in Figure 5.More recently, Angelinelli et al. (2021) found 1.05 <  < 1.15 in a sample of clusters from the non-radiative Itasca cluster sample (Vazza et al. 2017).Our measurements are broadly consistent with these values.Specifically, we find that nearly all relaxed and two-thirds of unrelaxed clusters in our sample are consistent with the Omega500 CSF predictions.The remaining third of unrelaxed clusters have upper limits lower than predicted in these simulations.Clumping in non-radiative runs in Omega500 is systematically higher and consistent with ∼ 70 per cent of the relaxed and ∼ 45 per cent of unrelaxed clusters from our sample.As for the Itasca simulations, we find that 15 out of 80 clusters in our sample have clumping factors constrained below 1.05, the rest are consistent with the predicted range.
The non-thermal pressure fraction, shown in Figure 6, has also been measured in several cosmological simulations.Battaglia et al. (2012); Nelson et al. (2014b); Shi et al. (2015) find that, broadly,  kin / tot ranges between ∼ 5 per cent and 25 per cent in the bulk ICM within ∼ 2500c .In our sample, 76 clusters (95 per cent of our sample) could fall within this range; among these, 6 clusters with the tightest constraints: A2744, the Bullet cluster, A520, MACS0717, 4C+37.11,and RXC0528.The 4 clusters with kinetic pressure fractions below this range, those mentioned in Section 3.2, have exceptionally low  kin / tot .Our measurements are, therefore, consistent with these simulations.

Turbulent re-acceleration of cosmic rays
Other than the X-ray thermal emission, the ICM also exhibits nonthermal synchrotron emission that is bright in radio (e.g., van Weeren et al. 2019).In particular, volume-filling radio halos are thought to be powered by turbulent re-acceleration of seed relativistic cosmic ray electrons (CRe) (e.g., Petrosian & Bykov 2008;Brunetti & Lazarian 2011;Brunetti 2016;Brunetti et al. 2017).From this model, one may expect a correlation between the observed radio power ( 1.4 GHz ) and ICM turbulence (or kinetic energy) (e.g., Bonafede et al. 2018), likely with significant scatter due to cluster-to-cluster variation of the seed CRe population (Nishiwaki & Asano 2022).
This relationship has been previously investigated (with mixed results) by Eckert et al. (2017) and Zhang et al. (2022).Eckert et al. (2017) found a strong correlation between the  1D,k (at 1/ = 660 kpc) and  1.4 GHz , finding a Pearson correlation coefficient of  : ≈ 0.8 in a sample of 25 radio halo detections (their sample also included 26 clusters with upper limits on  1.4 GHz ).Zhang et al. (2022), however, were unable to find any significant correlation between the same measurements in a sample of 11 halos, with  : ≈ 0.4 (1/ ≈ 400 kpc).Both studies inferred velocities indirectly through a similar method used in this work.The correlation has also been explored within individual clusters.Bonafede et al. (2018) found kinetic energy in MACS0717 ∼ 30 per cent higher in the region with strong radio emission compared to the radio-quiet region.In our sample, 22 clusters have detected radio halos at 1.4 GHz, in addition to 5 with upper limits (see Table B2 for details and references).It is interesting to check the radio power -turbulence relation, especially since we are probing smaller scales that are least affected by the choice of underlying model, merging substructures, and largescale asymmetries (Section 2.5).At multiple length scales, we compare the characteristic velocity to radio power in these 22 clusters, calculating  : at each scale.Representative distributions of  1D,k and  1.4 GHz are plotted in Figure 8, where velocity is measured at 0.3  2500c (blue circles) and 0.6  2500c (pink squares).Galaxy clusters with radio halo upper limits are also plotted as unfilled points.
To calculate uncertainties on  : , we assume the measurement uncertainties on  1D,k and  1.4 GHz are approximately Gaussian and randomly sample this distribution, finding the standard deviation on this randomized correlation coefficient.Additionally, we estimate the sample variance of  : via a bootstrap, combining these uncertainties in quadrature.The latter source of uncertainty is much more significant at all scales.
Across all length scales where velocities are measurable in most clusters, we find a mostly scale-independent correlation coefficient of  : ≈ 0.45 ± 0.3.This correlation is significantly weaker than that found by Eckert et al. (2017), and similar to that of Zhang et al. (2022).Given our larger sample size than the latter work, we can confirm a weak positive correlation at the  < 0.05 level between scales of 1/ ≈ 0.9 to 0.3  2500c (see Figure 9).
We stress that though we see a weak positive correlation in our full sample, this correlation is driven by a single outlier cluster, MACS0717, which has the highest characteristic velocity in our sample at the scales we consider here, and the second highest radio power.Removing MACS0717 from the calculation results in a correlation coefficient consistent with  : = 0 at all wavenumbers.A similar but less drastic reduction in  : is seen when the Bullet cluster is removed from the sample.These two clusters are present in the Eckert et al. (2017) sample, but not in that of Zhang et al. (2022), suggesting that the strong correlation found in the former work may also be significantly driven by these same clusters.Confirmation of this correlation, therefore, requires a larger sample of radio halo and direct velocity measurements.
Given the large uncertainties on  : , a strong correlation cannot be excluded.To further investigate, we attempt to fit the  1D,k - 1.4 GHz distribution to a power-law relationship.Using the Bayesian MCMC routine LinMix (Kelly 2007), which includes clusters with upper limits on  1.4 GHz , we find that the correlation between these variables is not strong enough to recover a best-fit slope.The slope varies significantly, between ∼ 2 and ∼ 10, depending on the numerical method used.This is also the case when MCMC maximum-likelihood fitting and orthogonal distance regression are used.
It is well-established that  1.4 GHz strongly correlates with cluster mass (e.g., Cassano et al. 2023Cassano et al. , 2013)).This correlation could contribute to the radio power -velocity correlation.Therefore, besides radio power, we consider  1.4 GHz / 2500c , which is a proxy for the efficiency of radio emission generation.This characteristic is compared to the turbulent energy dissipation rate  K ∝  3 1D,k  and M 2 1D,k , which traces the ratio of kinetic to thermal energy in the ICM.For each of these cases, we calculate the Pearson correlation coefficient across a range of length scales from 1 to 0.15  2500c , shown in the lower panel of Figure 9.In the upper panel of this figure, we plot the number of clusters in which  1D,k is observable at each scale.We additionally test other combinations of these variables, which do not change our conclusions.
The strongest correlation we found is when comparing radio emission efficiency to the turbulent dissipation rate, with   :/ ≈ 0.55±0.3,significant at the  < 0.05 level at length scales larger than 0.25  2500c .This correlation is also significantly driven by outlier clusters, as described above.The comparison of  1.4 GHz / 2500c to M 2 1D,k is consistent with no correlation at all scales.Further investigation of the relationship between ICM gas motions and non-thermal emission could be improved by larger samples of radio halo detections and velocity measurements, including direct velocity measurements with XRISM that will soon become available.We note that  1.4 GHz is measured across the entire extent of the radio halo (with typical radii of ∼ 1 Mpc), while our characteristic velocities are measured within  2500c (∼ 500 kpc).The effect of this mismatch is expected to be small, as typically ∼ 90 per cent of the total emission from radio halos is found within  2500c (Botteon et al. 2022).Additionally, the radio power measured at 1.4 GHz does not reflect the total luminosity of the halo, as there is evidence for variation in spectral slope between halos (van Weeren et al. 2019).Moreover, in some cases, significant areas around prominent substructures are removed.Thus, the comparison between our X-ray measurements and radio powers found in the literature is not one-to-one and could be improved in future works.
If substructures are included in the estimations of  1D,k , the strength of all correlations increases on length scales smaller than ∼ 0.8 2500c . : and   :/ both rise to ∼ 0.7±0.3.Thus, we find that improper treatment of X-ray substructures may lead to spurious correlations that do not reflect the state of the bulk ICM.

Constraints on effective viscosity
Threaded by a weak, likely tangled ∼ few  magnetic field (e.g., Carilli & Taylor 2002), the ICM is prone to the firehose and mirror plasma instabilities.These produce fluctuations, which can scatter particles, increasing their effective collisional rate and, therefore, decreasing the effective viscosity of the ICM (Schekochihin & Cowley 2006).At the same time, the gyroradius of ions is many orders of magnitude smaller than the typical Coulomb mean free path.Therefore, particles are confined to magnetic field lines, making the transport of heat and momentum anisotropic with respect to the local magnetic field (Kunz et al. 2011).
Observations of prominent structures within the ICM support these predictions.For instance, the width of interfaces between two gas phases in cold fronts is typically narrower than the mean free path of particles, indicating suppressed diffusion, conduction and mixing across the interface (e.g., Markevitch & Vikhlinin 2007;Werner et al. 2016;Zuhone & Roediger 2016).The presence of Kelvin-Helmholtz instabilities in cold fronts (e.g., Werner et al. 2016;Ichinohe et al. 2017;Su et al. 2017;Hu et al. 2021) and the morphology of stripped tails of infalling galaxies and groups (e.g., Kraft et al. 2017) also support suppressed effective viscosity, typically ≲ 5 per cent, with respect to the isotropic Spitzer level ( 0 , Spitzer 1962).Probing the effective viscosity in the bulk ICM has recently become possible through the measurements of density fluctuation power spectra in the Coma cluster on scales comparable to the Coulomb mean free path ( mfp ).The shape of the spectra indicated  eff / 0 ≲ 10 per cent regardless of the level of thermal conduction (Zhuravleva et al. 2019).With our measurements, it is possible to expand the latter analysis to a sample of 16 clusters, for which (/)  are measured across a broad range of scales close to the mean free path.
Using deprojected profiles of ICM temperature and density, we calculate radial profiles of  0 (see Sarazin 1988, for details).From our measurements of the scale-dependent characteristic velocity and under the assumption of Kolmogorov turbulence, we additionally find the turbulent energy dissipation rate (  ) as   = 2 (3  /2) 3/2  3 1D,k , where   ≈ 1.65 is the Kolmogorov constant (Zhuravleva et al. 2014a).Given  0 and   we find the median within our region of interest.The values of   we find are listed in Table B2.
Using these characteristics, we rescale power spectra of all clusters in our subsample4 , see Figure 10 (purple-orange curves), and obtain the weighted mean value (the gray region).We also plot the earliermeasured power spectra in the Coma cluster (gray hatched regions).For comparison with simulations, the equivalent spectra of passive scalars from direct numerical simulations (DNS) of hydrodynamic turbulence (Gauding et al. 2014) are also plotted with different colors.The passive scalar is characterized by the thermal Prandtl number , the ratio of kinematic viscosity and thermal conductivity.We use the DNS simulations for  = 1, 0.25 and 0.11 5 .We note that the Prandtl number of density fluctuations in the ICM is generally unknown, however it is expected to be less than unity (Kunz et al. 2022), thus we vary  in this range.The DNS spectra predict that the passive scalar spectrum should exponentially decline/steepen on scales probed with our observations.This steepening is not observed in any measurement of (/)  from this subsample, confirming earlier results based on the Coma cluster and suggesting that suppressed transport in the bulk of the gas is a general property of the ICM.
Following the Coma analysis, we can use these spectra to place upper limits on  eff .Suppressing  eff lowers   , essentially moving the DNS spectra to higher wavenumbers ( K ) such that the slope of measured (/)  eventually matches the slope of the DNS spectra.We use an MCMC chain with the Kolmogorov microscale and the normalization as free parameters to find the highest value of  eff / 0 to compare these spectra, however this normalization is arbitrary. 5We also compared the observed power spectra with the DNS simulations from Watanabe & Gotoh (2007); Yeung et al. (2002Yeung et al. ( , 2014)), and Ishihara et al. (2016) and confirmed that our conclusions are robust against the specific choice of the DNS simulations as long as the Reynolds number is sufficiently large.that provides a reasonable fit to the observed spectrum.These upper limits are presented for each object as green arrows in Figure 11.Additionally, we find the weighted average and standard deviation of (/)   K ( K  K ) −1/3 in this subsample and find an upper limit for the average spectrum, which is plotted as a dashed black line in the same figure.
We do not find any clusters in this subsample that are consistent with Spitzer-level effective viscosity.With  = 1.0, viscosity is suppressed by at least a factor of 2 in nearly all clusters.With  = 0.11, the suppression is at least a factor of 10 in all clusters, and a factor of 20 in all but 3 clusters.On average, the effective viscosity  eff / 0 is below 13 per cent in the case of  = 1.0, ≲ 5 per cent for  = 0.25, and ≲ 2 per cent for  = 0.11.The strongest effective viscosity suppression in the sample is found in A2390, with  eff / 0 ≲ 3 per cent for all Prandtl numbers we considered here.This demonstrates that, indeed, magnetic fields modify the behavior of particles in the ICM, suppressing transport properties of the bulk of the gas, and that this is a general property of ICM plasmas.
We note that the Δ-variance method is designed to measure smooth, power-law power spectra.The dotted lines in Figure 10 show how exponential cutoffs are modified by the Δ-variance method (for details, see Arevalo et al. 2012;Zhuravleva et al. 2019).Accounting for this correction does not affect the viscosity constraints shown in Figure 11.
Additionally, we note that we use the median measurement of the Kolmogorov microscale.Within the regions we consider in this work,  K typically varies by a factor of ∼ 2, which could increase the upper limit on  eff / 0 .If the minimum value of  K is used instead of the median, we still detect effective viscosity suppression in all systems except the Bullet cluster and A2142 (for the  = 1.0 case).A2390 has the largest variance, where our upper limit on  eff / 0 increases by a factor of ∼ 4.5 if the minimum value of  K is used instead of the median, however  eff / 0 is still ≲ 13 per cent for  = 1.0, even in this conservative limit.

CONCLUSIONS
By measuring the power spectra of X-ray surface brightness fluctuations in a sample of 80 galaxy clusters observed with Chandra, we obtain the amplitude of density fluctuations and, indirectly, characteristic velocities in cluster regions dominated by merger activity, i.e., outside the cool-cores and within the  2500c radius.Each cluster is treated individually to ensure robust measurements at the wavenumbers/length scales we present.Observational effects (e.g., Poisson noise, PSF effects, contribution of unresolved point sources) and systematic uncertainties induced by the choice of unperturbed model and cluster substructures are carefully tested and removed.Galaxy clusters are categorized into relaxed, intermediate, and unrelaxed merger states based on their X-ray morphology and a literature review.Scale-dependent amplitudes of density fluctuations, (/)  , are averaged within each morphological category.Our main findings are summarized below.
• We observe a mild increase in (/)  measured in the bulk ICM from relaxed to unrelaxed clusters.E.g., at 1/ = 0.5 2500c , (/)  ∼ 12 ± 2 per cent in relaxed clusters and ∼ 19 ± 5 per cent in unrelaxed ones.This suggests that fluctuations in the bulk ICM are not strongly dependent on recent merger activity.
• When prominent substructures are included in the measurement, density fluctuations in merging clusters (∼ 40 per cent at 1/ = 0.5 2500c ) are significantly higher than in relaxed clusters (∼ 12 per cent).Between length scales of ∼ 0.2 to 0.8 2500c (∼ 80 to 450 kpc), (/)  measured including substructures could be used as an indicator of cluster dynamic state, complementing other commonlyused methods.
• Power spectra of density fluctuations measured on scales from 1/ ∼ 1 to 0.2 2500c in intermediate and unrelaxed clusters are in good agreement with theoretical predictions from cosmological simulations considered here.
• On length scales smaller than ∼ 0.75 2500c , density fluctuations measured in relaxed clusters are systematically higher by a maximum factor of ∼ 2.5 compared to similar clusters in cosmological simulations.The effects of observational biases and numerical viscosity are likely minimal on these scales, as discussed in Section 4.1.1.The tension could be due to missing radiative and feedback physics in considered simulations.Fluctuations in relaxed clusters are most sensitive to these physics compared to merging clusters, in which the amplitude of fluctuations is likely dominated by ongoing or recent mergers.For this reason, relaxed clusters are the most useful probes for additional, multi-scale physics in simulations compared to clusters in other dynamic states.
• Estimated from power spectra of density fluctuations, clumping factors typically vary between 1.01 ≲  ≲ 1.10 within the  2500c region in our sample, with a few outliers (the Bullet and MACS0717 clusters).The measured clumping factors are broadly consistent with predictions from cosmological simulations, and in good agreement with other existing measurements within a similar region.
• Using a statistical relation between the amplitude of density fluctuations and velocity (with a recently-updated proportionality coefficient for clusters in different dynamical states), we find that the sample-average one-dimensional velocities are ∼ 150 (250) km s −1 at 1/ = 0.5 2500c for relaxed (unrelaxed) clusters.Within the uncertainties, the average power spectrum slope is consistent with a Kolmogorov cascade, although the uncertainties are large.
• Characteristic velocities (scale-independent) in the range of a few 100 km s −1 are found in essentially all clusters.The lowest velocities,  1D < 200 km s −1 , are found in A2597, A1650, A2626, A2199, and A2029.The highest velocities,  1D > 350 km s −1 , are found in the Bullet cluster, A1916, RBS1316 and MACS0717.These velocities await confirmation via direct observation by XRISM.
• We confirm kinetic pressure fractions in the bulk ICM of ≲ 10 per cent in 22 clusters, and ≲ 5 per cent in 6 of these: A2219, A1650, PKS0745, A1835, A2199, and A2029.The lowest is found in A2029.Note these values likely do not represent the total kinetic pressure support unless the cluster is unaffected by merging substructures (relevant to most relaxed clusters in our sample).These values of non-thermal pressure support are consistent with predictions from cosmological simulations.
• Using a subsample of galaxy clusters that host extended radio halos, we measure a weak positive correlation between characteristic velocity and radio power between length scales of 0.9 ≳ 1/ ≳ 0.3 2500c , consistent with Zhang et al. (2022) but weaker than Eckert et al. (2017).We additionally test correlations between the "radio emission efficiency" ( 1.4 GHz / 2500c ), turbulent energy dissipation rate (  ), and the ratio of kinetic to thermal energy (M 2 1D,k ).The strongest correlation we find is between   and radio emission efficiency, with the Pearson   :/ ≈ 0.55 ± 0.3.
• While these correlations are somewhat significant in the full sample, it is important to note that they are driven by a few outlier clusters.Removing the Bullet cluster or MACS0717 from the sample significantly reduces the correlation, such that the Pearson test is consistent with zero.Additionally, the results can be significantly affected by the inclusion/exclusion of substructures and the choice of underlying model.This demonstrates the need for upcoming direct velocity measurements with XRISM to investigate this correlation further.
• Using a subsample of 16 clusters from our sample, we constrain effective viscosity in the bulk ICM by comparing the shape of the power spectra of density fluctuations with DNS simulations.We find that the spectra of all 16 clusters are significantly shallower than an exponential cutoff predicted by DNS simulations on the same scales.This suggests that the effective viscosity in the bulk ICM is suppressed by magnetic fields.The estimated suppression is  eff / 0 ≲ 13 per cent for  = 1.0, ≲ 5 per cent for  = 0.25, and ≲ 2 per cent for  = 0.11.The lowest effective viscosity is found in A2390, with  eff / 0 ≲ 3 per cent for  = 1.0.These results are consistent with previous constraints in the Coma cluster and suggest that suppressed effective viscosity is a general property of the ICM plasmas.This work demonstrates the power of scale-dependent density fluctuations as tracers of a diverse combination of astrophysics and plasma physics.New (XRISM) and future (Athena, LEM) X-ray observatories will be able to directly measure gas velocities and structure functions, calibrating the (/)  - 1D,k relationship we employ in this work to confirm and expand these results.Even in this upcoming era of X-ray astronomy, the density fluctuations method will remain the observationally "cheapest" (most easily obtained) method of characterizing ICM dynamics.Improving our understanding of the relationship between density fluctuations and gas velocities will be a subject of future numerical and observational works.
Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.Additional partial support has been provided by the Alfred P. Sloan Foundation through the Sloan Research Fellowship.IZ was partially supported by a Clare Boothe Luce Professorship from the Henry Luce Foundation.WF acknowledges support from the Smithsonian Institution, the Chandra High Resolution Camera Project through NASA contract NAS8-03060, and NASA Grant 80NSSC19K0116 and Chandra Grant GO1-22132X.RJvW acknowledges support from the ERC Starting Grant ClusterWeb 804208.

APPENDIX A: ASYMMETRIC UNPERTURBED MODELS
We measure perturbations in density relative to an underlying gas distribution to obtain gas velocities in the ICM.In order to do this, we must derive the best-fitting unperturbed model that reflects the global distribution of gas within the cluster.The nominal model we use is an elliptical single or double -model, which has been shown to describe the ICM within our region of interest (Cavaliere & Fusco-Femiano 1978;Mohr et al. 1999;Kay & Pratt 2022).However, as there are many clusters in our sample undergoing mergers, the symmetric -model often fails to describe the gas distribution.One such example is the Bullet cluster with a plane-of-sky merger that has created a strong asymmetry in the SB.The peak of X-ray emission from the main cluster is due east of the eponymous "bullet" cold front.However, fitting the appropriate -model to this cluster creates a strong asymmetry, clearly seen in the top left panel of Figure A1.To create a model that better reflects the large-scale morphology of the cluster, we patch the unperturbed -model.The residual image  =  net /, where  net is the background-subtracted X-ray counts and  is the product of Chandra exposure map and the symmetric model, is smoothed with a large-scale Gaussian kernel , as is the mask .The symmetric model  sym is then modified as The residual image based on this asymmetric model is shown in the top left panel of Figure A1.
In practice, we test many different scales of patching, examining each resultant residual image to find which scale is sufficient to visually remove the large-scale asymmetry without over-correcting the model.Our final spectra only include scales that are largely unaffected by the final choice of patching scale (see the shaded regions in the bottom panel of Figure A1).

Figure 1 .
Figure 1.Our sample of 80 galaxy clusters in which we measure density fluctuation power spectra:  2500 vs. cluster redshift.Marker shapes indicate cluster morphology: "relaxed" (circles, cool-core systems with no signs of recent mergers), "intermediate" (squares, evidence of a minor merger), or "unrelaxed" (triangles, recent or ongoing major merger, no cool-cores), based on a literature review and visual inspection of the X-ray images, see Section 2.2 for details.Color indicates flare-removed (clean) Chandra exposure time in ks.Top: histogram of cluster redshifts.Right: histogram of cluster masses within  2500c .

Figure 2 .
Figure 2. Representative examples of clusters from relaxed, intermediate, and unrelaxed subsamples.Exposure-corrected Chandra X-ray (0.5-3.5 keV band)images are shown, lightly smoothed with a Gaussian for display purposes.Solid white circles represent  2500c , while dashed circles mark the cool-core radius if the cluster has one.A1835 is a relaxed, cool-core cluster that does not show any signs of mergers(Schmidt et al. 2001;Govoni et al. 2009).A401 is a cool-core cluster, however, it is undergoing a merger as evidenced by the sloshing cold front to the south (white arrow), the proximity of the cluster to A399, and the presence of a giant radio halo(Murgia et al. 2010).For these reasons, we place it in the intermediate category.A2034 is a clear example of an unrelaxed, merging cluster as it lacks a cool-core, shows irregular morphology, contains a merger shock north of the cluster center (white arrowOwers et al. 2013), and hosts a radio halo(Botteon et al. 2022).

Figure 3 .
Figure 3.Effect of bright substructures on measurements of ( /)  in the Bullet cluster.Solid lines show the mean value of ( /)  , corrected for Poisson noise and effects of the Chandra PSF, with 1 uncertainties (upper and lower limits).Shaded regions show wavenumbers where the measurement is unaffected by observational effects.Inset shows the exposure-corrected image of the  2500c region.The orange spectrum is measured in the whole region, while blue and green spectra are measured when their respective ellipses (see inset) are excluded.The exclusion of the blue ellipse is sufficient for characterizing fluctuations in the bulk ICM.

Figure 4 .
Figure 4. Middle row: Amplitude of gas density fluctuations vs. wavenumber (bottom x-axis) and length scale (top x-axis), in units of  2500c .Each line represents one cluster, where ( /)  is measured within a region of radius  2500c (excluding cool-cores).Clusters are divided into relaxed (left), intermediate (middle), and unrelaxed (right) subsamples.Only scales largely unaffected by systematics (as discussed in Section 2.7) are shown.Line color corresponds to the measurement's relative uncertainty, derived from 1 Poisson uncertainties (see color scale at right).Subsample-averaged amplitudes at each length scale, plus/minus one standard deviation, are shown in grey.Bottom row: Velocity power spectra for each cluster, with equivalent notations.The Kolmogorov slope is plotted as a guide (dashed black line).Two highlighted clusters, A2029 and A3667, will be observed with XRISM during the PV phase.Top row: Number of clusters with measured ( /)  at each scale within each subsample.Prominent substructures are removed from the analyzed images to measure fluctuations in the bulk ICM.For the version with included substructures, see Figure B2.

Figure 5 .
Figure 5. ICM clumping factors obtained by integrating ( /)  for each cluster in our sample.In pink, we integrate equation (8) from  min to  max , providing a lower limit on each value.We additionally fit the observed spectrum to a power-law via MCMC and integrate from  = 1/ 2500c to  = ∞.The median values, along with upper and lower uncertainties, are presented in blue with unfilled markers.These values could be interpreted as upper limits.Clumping factors measured in the Omega500 cosmological simulations(Nagai & Lau 2011;Zhuravleva et al. 2013) are plotted in shaded regions in their corresponding morphological categories.Grey regions show simulations where radiative processes (e.g., cooling and star formation/feedback) are included, while yellow regions show results from non-radiative simulations.
-scale density fluctuations Scale-dependent amplitudes of density (or any other characteristic) fluctuations can be used as sensitive probes of multi-scale, nonlinear physics and sub-grid prescriptions in numerical simulations.An equivalent method of measuring density fluctuations has been recently applied to a representative sample of 78 galaxy clusters (Zhuravleva et al. 2023) from the Omega500 cosmological simulations (Nagai et al. 2007a,b; Nelson et al. 2014b).Similar to our

Figure 6 .
Figure 6.Characteristic one-dimensional Mach numbers (left y-axis) and kinetic pressure fractions (right y-axis) for each cluster in our sample.Notations are equivalent to Figure 5.

Figure 7 .
Figure 7. Sample-averaged ( /)  measured in the bulk ICM of relaxed, intermediate, unrelaxed, and all clusters.Observational measurements from this work are plotted in blue, with 1 uncertainties in the hatched regions.Measurements from Omega500 cosmological simulations (for power spectra, see Zhuravleva et al. 2023) are plotted in pink.The blue bracket represents length scales where we confirm a divergence between observed and simulated measurements in relaxed systems.Black dashed lines represented a  1/3 power-law relationship for guidance.

Figure 8 .
Figure 8. Radio halo power at 1.4 GHz (see Table B2 for references) versus characteristic velocity  1D,k at 1/ = 0.6 (blue circles) and 0.3  2500c (pink squares).Unfilled squares/circles indicate upper limits on radio halo power.The Pearson correlation coefficients  : at the two different length scales are shown in legends.Uncertainties on  : are calculated via a bootstrap and by randomly sampling measurement uncertainties on  1D,k and  1.4 GHz .For clusters with  1D,k measured at both scales, paired pink and blue points are offset for visual clarity.

Figure 9 .
Figure 9. Pearson correlation coefficient , with 1 uncertainties, as a function of wavenumber/scale, comparing  1D,k to  1.4 GHz (orange),  K to  1.4 GHz / 2500c (blue), and M 2 1D,k to  1.4 GHz / 2500c (green).Uncertainties on  are calculated via a bootstrap and by randomly sampling measurement uncertainties on  1D,k and  1.4 GHz . = 0 is plotted as a solid black line.The number of clusters where  1D,k is measurable at each scale is plotted in the top panel.The p-value of each Pearson test is plotted in the bottom panel, with  = 0.05 plotted as a dashed black line.

Figure 10 .
Figure 10.Observed rescaled ( /)    in the 16 clusters in which we constrain effective viscosity.Notations are equivalent to Figure 4, with the sample-averaged values ±1 uncertainties shown with the grey region.The equivalent passive scalar power spectra from DNS simulations of hydrodynamic turbulence for three thermal Prandtl numbers are shown in solid lines.The dotted lines show how the exponential cutoff of the DNS spectra would be modified if the Δ-variance method of calculating power spectra was used.Hatched regions show upper and lower limits on rescaled ( /)  taken from the central and outer regions of the Coma cluster(Zhuravleva et al. 2019), for comparison.The slope of the observed spectrum is shallower at the equivalent   K than any of the DNS spectra, suggesting a suppressed effective viscosity in the ICM.

Figure 11 .
Figure 11.Upper limits on effective viscosity suppression,  eff / 0 , relative to a Spitzer value in a subsample of 16 galaxy clusters for different thermal Prandtl numbers.Sample-averaged upper limits are plotted as dashed black lines.

Figure A1 .
Figure A1.Effects of an asymmetric surface brightness model on measurements of ( /)  in the Bullet cluster.Top: Residual images within  2500c of the cluster using a spherical (left) and asymmetric (right; patched at a scale of ∼ 400 kpc = 90") -model for the unperturbed model.Prominent substructures have been removed as in Figure 3. Bottom: The orange power spectrum corresponds to the left residual image, while the blue power spectrum corresponds to the right.The green spectrum is patched on a smaller scale to ensure that measured wavenumbers are not strongly affected by the choice of patching scale.Notations are equivalent to Figure 3.

Figure B1 .
Figure B1.One example of a measured power spectrum (orange, solid lines) taken from the Bullet cluster compared to the MCMC-fitted power spectrum (black, dashed lines).1 uncertainties are shown in shaded regions.Notations are equivalent to Figure 3.

Figure B2 .
Figure B2.Equivalent to the top and middle rows of Figure 4, but with substructures included in the calculation of ( /)  .

Table B1 :
Name, redshift, Chandra ObsIDs used in this work, and cleaned exposure time (in ks) for each cluster in our sample.

Table B2 :
Cluster properties: Dynamic state (Section 2.1), chosen center coordinates, median ICM temperature in the region of interest, cluster mass ( 2500 ) and radius ( 2500c ), cool-core radius, total power of radio halo (if present, Section 4.2.1), and estimated Kolmogorov microscale (Section 4.2.2).Temperatures marked with † are measured in this work.Temperatures marked with * are taken from projected temperature profiles.