Issue 
A&A
Volume 559, November 2013



Article Number  A78  
Number of page(s)  18  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201322295  
Published online  18 November 2013 
Constraining turbulence and conduction in the hot ICM through density perturbations^{⋆}
Max Planck Institute for Astrophysics, KarlSchwarzschildStrasse 1, 85741 Garching, Germany
email: mgaspari@mpagarching.mpg.de
Received: 16 July 2013
Accepted: 3 October 2013
Turbulence and conduction can dramatically affect the evolution of baryons in the universe; current constraints are however rare and uncertain. Using 3D highresolution hydrodynamic simulations, tracking both electrons and ions, we study the effects of turbulence and conduction in the hot intracluster medium. We show how the power spectrum of the gas density perturbations (δ = δρ/ρ) can accurately constrain both processes. The characteristic amplitude of density perturbations is linearly related to the strength of turbulence, i.e. the 3D Mach number, as A(k)_{δ,max} = c M, where c ≃ 0.25 for injection scale of 500 kpc. The slope of A_{δ}(k) in turn reflects the level of diffusion, dominated by conduction. In a nonconductive medium, subsonic stirring motions advect density with a similar nearly Kolmogorov cascade, E_{δ}(k) ∝ k^{− 5/3}. Increasing conduction (parametrized via the magnetic suppression f = 10^{3} → 1) progressively steepens the spectrum towards the Burgerslike regime, E_{δ}(k) ∝ k^{2}. The slope is only weakly dependent on M. The turbulent Prandtl number defines the dynamic similarity of the flow; at scales where P_{t} ≡ t_{cond}/t_{turb} < 100, the power spectrum develops a significant decay, i.e. conduction stifles turbulent regeneration. The transition is gentle for highly suppressed conduction, f ≤ 10^{3}, while sharp in the opposite regime. For strong conductivity (f ≥ 0.1), P_{t} ~ 100 occurs on spatial scales larger than the injection scale, globally damping density perturbations by a factor of 2−4, from large to small scales. The velocity spectrum is instead not much affected by conduction. The f ≥ 0.1 regime should also affect the appearance of Xray images, in which KelvinHelmholtz and RayleighTaylor rolls and filaments are washed out. In a stratified system, perturbations are characterized by a mixture of modes: weak/strong turbulence induces higher isobaric/adiabatic fluctuations, while conduction forces both modes towards the intermediate isothermal regime. We provide a general analytic fit which is applied to new deep Chandra observations of Coma cluster. The observed spectrum is best consistent with strongly suppressed effective isotropic conduction, f ≃ 10^{3}, and mild subsonic turbulence, M ≃ 0.45 (assuming injection scale at ~250 kpc). The latter implies E_{turb} ≃ 0.11 E_{th}, in agreement with cosmological simulations and linebroadening observations. The low conductivity corroborates the survival of sharp features in the ICM (cold fronts, filaments, bubbles), and indicates that cooling flows may not be balanced by conduction.
Key words: conduction / turbulence / hydrodynamics / galaxies: clusters: intracluster medium / methods: numerical / plasmas
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2013
1. Introduction
The intracluster medium (ICM) plays central role in the evolution of baryons in the universe. The ICM is the hot plasma filling the gravitational potential of galaxy clusters, the largest virialized structures in the universe. Since most of the cluster baryons reside in the ICM (~85 percent), this gaseous medium represents the crucible out of which essential astrophysical structures condense. It is often assumed that the ICM settles in hydrostatic equilibrium after the initial cosmological collapse in the potential well of the cluster. However, the ICM plasma is a remarkably dynamic entity, continuously perturbed by mergers, feedback processes (active galactic nuclei – AGN, supernovae – SNe), galaxy motions, and cosmological accretion, all shaping a chaotic and turbulent atmosphere.
Current Xray observations have in general hard time detecting surface brightness (SB_{x}) fluctuations, due to the significant level of Poisson noise dominating on small scales (several tens kpc for nearby clusters) and projection smearing effects. This has lead, throughout the past decades, to the common assumption that the ICM is a static entity, both in theoretical and observational work. Only in recent time, few observational investigations have started to focus on the perturbations in the ICM, thanks to deep Chandra or XMMNewton data. Schuecker et al. (2004) found that at least ~10 percent of the total ICM pressure in Coma hot cluster is in turbulent form. The spectrum of pressure fluctuations, in the range 40−90 kpc, appears to be described by a Kolmogorov slope. Very recently, Churazov et al. (2012, and in prep.) have analyzed very deep observations of Coma (650 × 650 kpc; Sect. 2.1). The characteristic amplitude of the relative density fluctuations reaches 5−10 percent, from small (30 kpc) to large scales (500 kpc), again resembling the Kolmogorov trend. Sanders & Fabian (2012) also studied density/pressure perturbations in a coolcore cluster, AWM7, finding an amplitude of ~4 percent, with a largescale spectrum shallower than Kolmogorov.
Cosmological simulations (e.g. Norman & Bryan 1999; Dolag et al. 2005; Nagai et al. 2007, 2013; Lau et al. 2009; Vazza et al. 2009, 2011; Borgani & Kravtsov 2011) indicate that subsonic chaotic motions are ubiquitous, with turbulent pressure support in the range 5−30 percent, from relaxed to merging clusters. On the other hand, largescale simulations have severe difficulty in studying the details of perturbations, due to the limited resolution, the adaptive mesh (AMR) derefinement, or the smoothed particle hydrodynamics viscosity. From a theoretical point of view, little attention has thus been paid to studying the role of density perturbations, and in particular the associated power spectrum, down to kpc scale. In idealized periodic boxes, Kim & Ryu (2005) showed that isothermal turbulence produce a Kolmogorov spectrum, progressively flattening with increasing Mach number. Even in the presence of weak magnetic fields, the power spectrum seems to retain the Kolmogorov slope, at least in ideal magnetohydrodynamics (MHD) simulations (Kowal et al. 2007). In multiphase flows (e.g. the interstellar medium – ISM), thermal instability generate a more complex nonlinear dynamics, inducing high smallscale density perturbations, still increasing with M (Kissmann et al. 2008; Gazol & Kim 2010).
In the current work, we intend to carefully study the power spectrum (or better, the characteristic amplitude) of 3D density fluctuations, driven by turbulent motions (Sect. 2.3) in a real galaxy cluster, Coma (Sect. 2.1). This run will set the reference model. We pay particular attention in modeling a realistic hot ICM plasma, as tracking both electrons and ions (2T; Sect. 2.5), and avoiding restrictive assumptions on the equation of state (e.g. isothermality). Ideal hydrodynamics is however not enough to study a consistent evolution of an astrophysical plasma as the ICM. The very high temperatures (~10^{8} K), combined with the low electron densities (~10^{3} cm^{3}), warn that thermal conduction may have a profound impact in shaping density inhomogeneities (Sect. 2.4).
The electron thermal conductivity of the ICM is a highly debated topic in astrophysics, and currently poorly (or not) constrained. In the standard picture of a uniformly magnetized plasma, classic Spitzer conduction (Sect. 2.4) is suppressed perpendicular to the Bfield lines (by a factor f ≲ 10^{12}), due to electron scattering limited by the Larmor radius. However, turbulent plasmas develop tangled magnetic fields with a chaotic topology. According to Rechester & Rosenbluth (1978) and Chandran & Cowley (1998), after ~30 times a random walk of the Bfield coherence length l_{B}, an electron in the ICM could be fully isotropized, leading to an effective isotropic suppression f ~ 10^{2}−10^{3}, still a substantially stifled conductive flux. Narayan & Medvedev (2001) and Chandran & Maron (2004) further argued that, in a turbulent plasma, field lines can be chaotically tangled even on scales <l_{B}, possibly restoring the effective conductivity up to f ~ 0.1 − 0.4 (the Spitzer value).
Past investigations have mainly focused on the role of conduction in balancing radiative losses. In order to prevent the cooling catastrophe, the level of thermal conduction requires to be substantial, f ≳ 0.1, or even impossible for several observed clusters, f > 1 (Kim & Fabbiano 2003; Zakamska & Narayan 2003; Voigt & Fabian 2004). Simulations also confirm the inefficiency of conduction (Dolag et al. 2004; Parrish et al. 2009), requiring other heating mechanisms, as AGN feedback (Churazov et al. 2000, 2001; Ruszkowski & Begelman 2002; Brighenti & Mathews 2003; Gaspari et al. 2011a,b, 2012a, 2013a) or turbulent mixing (Ruszkowski & Oh 2010, 2011). Overall, observations lean towards highly suppressed conduction (f ≲ 10^{3}), given the ubiquitous presence of cool cores, along with sharp temperature gradients linked to cold fronts (Ettori & Fabian 2000; Roediger et al. 2013; ZuHone et al. 2013), Xray cavities, or cold filaments (Forman et al. 2007; Sanders et al. 2013). However, these constraints do not provide the effective isotropic conductivity in the bulk of the ICM; in fact, the magnetic field lines tend to naturally align perpendicular to the temperature gradient in a turbulent medium (Komarov et al. 2013), hence strongly preventing the heat exchange between sharp fronts.
Scope of this work is to provide, for the first time, a global constrain on the conductive and turbulent state of the ICM, instead of relying on local features. This will be possible exploiting the power spectrum of density perturbations. After setting the physical and numerical framework (Sect. 2), we proceed step by step with controlled experiments (Sect. 3), assessing first the role of turbulence (weak, moderate, strong), and then, gradually increasing the effective thermal conductivity. The three features of the power spectrum unveil each a crucial aspect of the ICM state. The normalization results to be linearly related to the turbulent Mach number. The slope of the spectrum steepens from Kolmogorov to Burgers trend, with rising conductivity. The decay/cutoff of the spectrum is provided by a key recurrent threshold, t_{cond}/t_{turb} ≲ 100 (Sect. 3). In Sect. 4, we discuss important properties of the models, and provide a simple model to assess the effective conductivity and turbulence. We apply the prescription to new very deep observations of Coma cluster, constraining the actual state of the hot ICM. The results are summarized in Sect. 5. In the Appendices, we compare different methods to calculate the spectrum, and analytically study the βprofile in Fourier space. The increasing quality of future Xray data will provide a big opportunity to exploit and perfect this new modeling, and hopefully to lead to highprecision measurements of the ICM (“ICMology”).
2. Physics and numerics
2.1. Initial conditions: Coma galaxy cluster
Hot galaxy clusters are optimal systems to study the effects of thermal conduction and turbulence, due to the fairly low ICM densities and the substantial level of dynamical activity. The archetypal noncoolcore system is Coma cluster (Abell 1656). Given its proximity, brightness and flat Xray core, it is ideal to study density perturbations. Churazov et al. (2012) retrieved the characteristic amplitude of density fluctuations in Coma from deep XMM and Chandra observations, finding significant values up to 10 percent, while resolving scales of tens kpc.
In this study, we adopt Coma cluster as fiducial astrophysical laboratory, setting the density and temperature profile according to the most recent XMM observations (Fig. 1). An excellent fit to the radial electron density distribution is given by a single βmodel profile: (1)with central density n_{e,0} = 3.9 × 10^{3} cm^{3}, core radius r_{c} = 272 kpc, and β = 0.75. In Appendix B, we discuss the properties of the βprofile in Fourier space. The gas temperature is roughly isothermal in the core, declining at large radii as observed for the majority of clusters (e.g. Vikhlinin et al. 2006): (2)where T_{0} = 8.5 keV and r_{t} = 1.3 Mpc. The electron and ion temperature (Sect. 2.5) are initially in equilibrium. The combined property of high temperature and low density sets a perfect environment to study the role of conduction, since the thermal diffusivity is ∝ (Sect. 2.4). In addition, radiative cooling is ineffective, t_{cool} ~ 2.5 t_{H}. The simulated system is initialized in hydrostatic equilibrium. The gas temperature and density allow to retrieve the gravitational acceleration (dominated by dark matter). The resulting potential is appropriate for a massive cluster in the ΛCDM universe, with virial mass M_{vir} ~ 10^{15}M_{⊙} and r_{vir} ~ 2.9 Mpc (r_{500} ~ 1.4 Mpc).
Fig. 1
Initial conditions for our reference hot galaxy cluster, Coma. From top to bottom panel: electron temperature, number density and entropy (). The radius is normalized to r_{500} ≃ 1.4 Mpc. We extracted the red data points from recent deep XMM observations, reaching ~0.5 r_{500} (Lyskova et al., in prep.). 
Since deep observations of density perturbations reach at most 0.5 r_{500} (Fig. 1), we adopt a 3D box with a diagonal of ~2.4 Mpc. As turbulence is volume filling and since we are interested in the power spectrum, the best numerical approach is to use a fixed grid, without adaptive refinement. In fact, there is no trivial AMR criterium to apply due to the uniformly chaotic dynamics. Moreover, when the cube is derefined by more than 50 percent, we found that that the density spectrum has a significant decrease in power towards the small scales (producing a mock diffusivity), by over a factor of 2. Based on these tests, we warn that using largescale (cosmological) simulations will numerically steepen the slope of power spectra (density, velocity, etc.), as also shown by Vazza et al. (2009). Albeit computationally challenging, we thus run all the models with fixed grid and high resolution of 512^{3} (considering the implemented physics). We tested also 256^{3} runs, finding a very similar evolution and spectrum, though with double dissipation scale. The simulations are thus in the convergence limit.
The resolution is Δx ~ 2.6 kpc, i.e. roughly on the scale of the (unmagnetized) plasma mean free path. Going below this scale would formally require a kinetic approach. This also means that numerical viscosity is on the scale of the physical Spitzer viscosity, further corroborating the use of such resolution (Sect. 2.6). The total evolution time is typically ≳2 eddy turnover times, roughly the statistical steady state after the turbulent cascade. Boundary zones have Dirichlet condition, with value given by the largescale radial profile. Inflow is prohibited, in order to avoid any spurious wave altering the dynamics inside the domain.
2.2. Hydrodynamics
We use a modified version of the grid code FLASH4 (Fryxell et al. 2000) in order to integrate the 3D equations of hydrodynamics for a 2 temperature (electronion; Sect. 2.5) plasma, with the addition of turbulence (Sect. 2.3) and electron thermal conduction (Sect. 2.4):
where ρ is the gas density, v the velocity, e_{i} and e_{e} the specific internal energy of ions and electrons, P_{tot} the total pressure (ions and electrons), γ = 5/3 the adiabatic index. The mean atomic weight of electrons and ions is μ_{i} ≃ 1.32, μ_{e} ≃ 1.16, providing a total gas μ ≃ 0.62, appropriate for a totally ionized plasma with ~25% He in mass. The atomic weight determines also the specific (isochoric) heat capacity, c_{V} = k_{B}/ [(γ − 1) μm_{p}]. For numerical consistency, the total energy (e_{tot} = e_{i} + e_{e} + v^{2}/2) equation, in conservative form, is also integrated.
In order to integrate the hyperbolic part of the hydrodynamics equations, we use a robust third order reconstruction scheme (piecewise parabolic method) in the framework of the unsplit flux formulation with hybrid Riemann solver (Lee & Deane 2009). Albeit computationally expensive, this setup keeps at minimum the numerical diffusivity. We tested different Riemann solvers (e.g. HLLC, Roe), characteristic slope limiters (MinMod, Van Leer, Toro), and other parameters (e.g. Courant number, interpolation order). They give comparable results, although we note that lower order reconstruction schemes (e.g. MUSCL) are more diffusive and thus truncate the turbulent cascade on roughly two times larger scale.
2.3. Turbulence driving
Continuous injection of turbulence is modeled with a spectral forcing scheme that generates statistically stationary velocity fields (Eswaran & Pope 1988; Fisher et al. 2008; Gaspari et al. 2013b). This scheme is based on an OrnsteinUhlenbeck (OU) random process, analogous to Brownian motion in a viscous medium. The driven acceleration field is timecorrelated, with zero mean and constant root mean square, an important feature for modeling realistic driving forces. In the OU process, the value of the gas acceleration at previous timestep a^{n} decays by an exponential damping factor f = exp( − Δt/τ_{d}), where τ_{d} is the correlation time. Simultaneously, a new Gaussiandistributed acceleration with variance is added as (8)where G^{n} is the Gaussian random variable, ϵ^{∗} is the specific energy input rate, and is the updated acceleration. The six amplitudes of the acceleration (3 real and 3 imaginary) are evolved in Fourier space and then directly converted to physical space. In this approach, turbulence can be driven by stirring the gas on large scales and letting it cascade to smaller scales. This is an efficient approach as the alternative would involve executing fast Fourier transforms (FFTs) for the entire range of scales, where the vast majority of modes would have small amplitudes. Since ICM turbulence is always subsonic, we impose a divergencefree condition on acceleration, through a Helmholtz decomposition in Fourier space.
The physical quantity of interest is the resultant 3D turbulent velocity dispersion, σ_{v}, which drives the ICM dynamics. The driving of turbulence is intentionally kept simple as our goal is not to consider any specific stirring source, but to keep the calculation fairly general. For example, the (combined) sources of turbulence may be major or minor mergers, galaxy motions, AGN feedback or supernovae. Observations (Schuecker et al. 2004; Churazov et al. 2008; de Plaa et al. 2012; Sanders & Fabian 2013) and simulations (Norman & Bryan 1999; Lau et al. 2009; Vazza et al. 2009, 2011; Gaspari et al. 2012b; Borgani & Kravtsov 2011 for a review) show that ICM turbulent energies are in the range few −30 percent of the thermal energy, from very relaxed to merging clusters.
We test therefore three regimes of ICM turbulence, weak (M ~ 0.25), mild (M ~ 0.5), and strong (M ~ 0.75), corresponding to a ratio of turbulent to thermal energy of 3.5, 14 and 31 percent (E_{turb} ≃ 0.56 M^{2}E_{th}). This is achieved by adjusting the energy per mode^{1}ϵ^{∗} and correlation time τ_{d} (ϵ^{∗} ~ 5 × 10^{5} − 10^{3} cm^{2} s^{3} and 200 Myr, respectively). As long as different choices of these parameters result in the same velocity dispersion, the dynamics of the flow remains unaffected. We stir the gas only on large scales, with typical injection peak L ~ 600 kpc (in the last set of runs L′ ~ 300 kpc), letting turbulence to naturally cascade. This allows us to exploit the entire dynamic range of our box (512^{3}) and to better appreciate the effect of conduction, without being strongly affected by numerical diffusion. Notice that in few t_{eddy}, largescale turbulence is not able to eject a substantial amount of mass outside the box. Since turbulence is kept subsonic, dissipational heating, which is proportional to , is also secondary. Turbulent diffusion can instead effectively flatten the global entropy gradient, especially in the nonconductive runs. We are nevertheless interested in the relative variations of δρ/ρ, removing the underlying profile.
The characteristic time of turbulence is defined by the eddy turnover time at a given physical scale l (e.g. Fig. 8, magenta line). Extrapolating from the injection scale via the Kolmogorov scaling, σ_{v} = σ_{v,inj} (l/L)^{1/3}, yields (9)Finally, we note turbulence can be expressed as a diffusion process acting on entropy on sufficiently large scales, ≳l, albeit the equations are intrinsically hyperbolic. We define the effective turbulent diffusivity as (10)The transport of heat due to turbulent diffusion can be written as (11)where s = c_{V}ln(P/ρ^{γ}) is the entropy. We remark that turbulence diffuses entropy, while seeding perturbations in density and temperature. In our analysis and discussion, we assume a diffusion constant c_{t} = 1, but in real plasmas this uncertain value could be much lower (Dennis & Chandran 2005, and references therein).
2.4. Thermal conduction
In ionized plasmas such as the ICM, electrons conduct internal energy with a heating rate per unit volume given by (12)The thermal conductivity can be written as (Spitzer 1962; Cowie & McKee 1977) (13)where lnΛ_{ei} = 37.8 + ln [(T_{e}/10^{8} K) (n_{e}/10^{3} cm^{3})^{− 1/2}] is the Coulomb logarithm (the ratio of the largest to smallest impact parameter; e.g. Voigt & Fabian 2004), and f is the magnetic suppression factor (Sect. 2.4.1). The previous conductivity derives from the more significant expression (14)which points out that the characteristic length scale and speed of conduction is the electron mean free path and the electron thermal speed v_{e} = (3k_{B}T_{e}/m_{e})^{1/2}, respectively. Another important quantity is the isochoric diffusivity (cm^{2} s^{1}), defined as D_{Sp} = κ/c_{V,e} ρ = κ/1.5 n_{e} k_{B} ≃ 0.5 v_{e}λ_{e}.
2.4.1. Magnetic field suppression
The intracluster plasma is likely magnetized. Although the ICM magnetic field (few μG) appears dynamically unimportant compared to the thermal pressure, electrons and ions are anchored to the Bfield lines. The gyroradius or Larmor radius is many orders of magnitude lower than the Coulomb mean free path, hence the charged particles can diffuse only along the Bfield lines. On the other hand, a real atmosphere is always characterized by some degree of chaotic motions (e.g. AGN/supernovae feedback, sloshing, mergers, galaxy motions). Ruszkowski & Oh (2010) showed through 3D MHD simulations that very weak subsonic turbulence always induces a tangled magnetic field with small coherence length (l_{B} ≲ 10 kpc; Fig. 2). Faraday rotation measures also indicate small l_{B} (e.g. Kim et al. 1990). Therefore, on scales larger than the Bfield coherence length, the average suppression due to anisotropic conduction can be parametrized with the socalled f factor, as in our study. If l_{B} is instead very large (as in the idealized case of a parallel topology), this prescription is less accurate: parallel and perpendicular conductivity should be studied separately. The spectral analysis should be nevertheless not strongly affected by localized features (as fronts and filaments), weighing more volumefilling properties; Fourier spectra provide thus the conductivity in the bulk of the ICM.
Highresolution MHD^{2} simulations are required to assess the role of anisotropic conduction and will be tackled in a future work. Assuming a sufficiently tangled field, f is geometrically expected to be ⟨ cos^{2}θ ⟩ ≈ 1/3 (θ is the angle between the B field and the T gradient), as confirmed by turbulent MHD runs (Ruszkowski & Oh 2010). Microscopic effects and plasma instabilities can severely suppress the conductive flux down to f ~ 10^{3} (Rechester & Rosenbluth 1978; Chandran & Cowley 1998), although Narayan & Medvedev (2001) argue that in a highly tangled and chaotic field conduction may be restored to f ~ 0.2 (Sect. 1). Considering these uncertainties, we test a wide range of f values, ranging from the strongly to weakly suppressed regime, f ~ 10^{3} − 1.
From an observational point of view, deep data (e.g. Markevitch & Vikhlinin 2007) show sharp contact discontinuities in the ICM, leading to the conclusion that conduction is severely suppressed, f ~ 10^{3} (Ettori & Fabian 2000). However, these estimates focus only on local features. Simulations of sloshing motions (e.g. ZuHone et al. 2013) show that the magnetic field tends to remain perpendicular to the temperature gradient (a natural outcome of the frozenin property; Komarov et al. 2013), with a strong suppression across the front. It is still unknown what is the effective global conductivity of the turbulent ICM (and hence global heat diffusion), which we aim to constrain with our spectral method, comparing simulated and observed density perturbations (Sect. 4.3).
2.4.2. Saturation
Since the conductive flux depends also on the temperature gradient, electrons may conduct heat much faster than their thermal velocity in an unphysical way. This happens whenever the temperature scale height is smaller than the electron mean free path, l_{T} ≡ T/  ∇T  ≲ λ_{e}. In this regime, the conductive flux saturates at a value given by (Cowie & McKee 1977; Balbus & McKee 1982) (15)where α is an uncertainty factor representing microscopic processes in the magnetized plasma (as instabilities). Following the indication of Balbus (1986) based on plasma experiments, we set α ~ 0.1, although the exact value has no great impact on the dynamics. Saturation changes the nature of the equations, from parabolic to hyperbolic, yet we can define an effective diffusivity of the form D_{sat} =  ∇T/F_{sat} . Numerically, saturation is implemented via a smooth flux limiter on the diffusion coefficient: . We tested other types of limiters, as harmonic or min/max, without finding relevant differences. Saturation is only relevant in the unsuppressed regime, with typically less than 10 percent of zones saturated, while negligible for f ≲ 0.1 models.
The final diffusivity is important to determine the characteristic timescale of conduction (see Fig. 8, black line): (16)In order to integrate the diffusion equation, we initially used an explicit fluxbased scheme. However, the computational time becomes prohibitive since it is strongly limited by Eq. (16), allowing to integrate only few 100 Myr. It is thus essential to adopt the (unsplit) implicit solver, allowing for a fast yet accurate execution. The solver efficiently uses the HYPRE linear algebra package to solve the diffusion equation linked to electron thermal conduction (cf. FLASH4 guide for the validation tests). The associated boundaries are set in outflow or zerogradient mode.
2.5. 2T plasma: electron–ion equilibration
In astrophysical simulations, it is widely assumed that the plasma has one single temperature, T_{e} ≈ T_{i}. However, this approximation is only good for a relatively cold medium. For hot clusters, especially noncoolcore systems, the ion–electron equilibration time due to Coulomb collisions can be t_{ei} ≳ 50 Myr. Since conduction operates on the Myr scale and only on electrons, it is important to follow the evolution of both the electron and ion temperature (or internal energy; Eqs. (5), (6)). The heat exchange rate (erg s^{1}) between ions and electrons is given by (17)where we choose the widely used Spitzer electron–ion equilibration time for a fully ionized plasma (Huba 2009): (18)The equilibration time is dominated by ∝(m_{i}T_{e})^{3/2}. Therefore, in systems with characteristic temperature ≲6 keV (and dense cores), the equilibration time is comparable or less than the unsuppressed conduction timescale (∝T^{− 5/2}), even on few kpc, t_{ei} ≲ 1 Myr. Neglecting equilibration in the strongly conductive runs, the ions would be forced to be quickly isothermal, inducing spurious features. In the more realistic evolution, turbulent motions displace the ions before having time to fully equilibrate with electrons, leading to a gentler equilibration. To our knowledge, this is the first attempt to study ICM turbulence and conduction with a twotemperature (2T) simulation, in analogy with highenergy physics studies.
2.6. Viscosity
Viscosity may in principle not be neglected when conduction operates, since both transport processes are intimately connected on the microscopic scale, and both altered by Bfield instabilities and topology. In the present work, we do not implement a direct physical viscosity, yet we point out two important arguments. Since viscosity is the transport of momentum due to ions, while the conductive flux is associated with electrons, viscous stresses are slower by at least 1−2 orders of magnitude compared with conduction (thermal v_{i} ≃ v_{e}/43), and should have a secondary role in damping density perturbations (see also Sect. 3). Further, we choose numerical resolution to be on the scale of the ion mean free path (λ_{i} ≃ λ_{e}); the flow velocity is also comparable to the characteristic velocity of viscosity. This implies that numerical diffusivity approximately reflects Spitzer viscous diffusion (∝λ_{i} v_{i}) – the unsuppressed dynamic viscosity in Coma is g cm^{1} s^{1} (cf. Reynolds et al. 2005). For example, the nonconductive spectrum in Fig. 2 indicates that the typical Reynolds number at injection is Re_{L} ≃ (L/λ_{T})^{2} ~ 500, where the Taylor scale λ_{T} marks the incipit of dissipation. Comparing the simulated Taylor scale with the Spitzer value, λ_{T} ≃ (ν_{Sp} L/σ_{v})^{1/2} ~ 80 kpc (ν_{Sp} = η_{Sp}/ρ), we see that the grid acts as an effective viscosity with f_{η} ~ 1/4. In future, we will test the effects of a varying and anisotropic physical viscosity.
2.7. Power spectrum of density perturbations
To study density perturbations, we adopt the characteristic amplitude, instead of the power spectrum P(k) (or energy spectrum E(k)), defined as (19)where (we typically use l = 1/k in kpc units). The characteristic amplitude is insightful, since the units are the same of the variable in real space. Since we are interested in the relative perturbations of density, δ ≡ δρ/ρ, A(k)_{δ} represents the typical level of fluctuations at a given scale k. The peak of A(k)_{δ} provides a good and simple estimate for the total amount of perturbations; the exact total variance can be computed integrating P(k) 4πk^{2}dk over the whole range of scales (the difference may reach a factor of 1.4, depending on the power slope).
We retrieve the relative density perturbations δ, dividing ρ by the background profile, δ = ρ/ρ_{b} − 1. For each snapshot, we execute a bestfitting routine to compute the new underlying βprofile (e.g. strong turbulent diffusion can lower the central density by a factor of 2 over few Gyr), minimizing the deviation between the data and the model. We note that the βprofile removal affects only the spectrum on very large scales. In principle, A(k)_{δ} could be studied even without removing the background, since perturbations start to dominate on large k (see Appendix B). The power spectrum of δ is finally computed with the “Mexican Hat” filtering (Arévalo et al. 2012) instead of performing Fourier transforms, which can lead to spurious features due to the nonperiodicity of the box (see Appendix A).
3. Results
Key properties of A_{δ}(k) for the simulated models: normalization (maximum), slope, and scale of significant decay.
We now describe the results of the simulated models, reminding that the main goal of the present investigation is to understand the role of turbulence and conduction in shaping the power spectrum of density perturbations, A(k)_{δ} (Sect. 2.7; we defer to future work the study of other statistics, as PDFs and structure functions). We are interested in the characteristic level of δρ/ρ perturbations, the slope of the spectrum, as well as any evident decline (or cutoff). Table 1 summarizes the key retrieved properties and serves as a guide for the analysis of A(k)_{δ}.
A key quantity for describing the evolution of perturbations is the ratio of the conduction and turbulence timescale (Eqs. (9) and (16)), which normalized to reference values of the unsuppressed conductive run results to be (20)We do not consider saturated conduction for this estimate, since the temperature gradient is not steep for the majority of the zones (the interpolated D_{cond} shall be used in this case for higher accuracy; Sect. 2.4). This key timescale ratio can be also seen as the Prandtl number applied to turbulence, instead of to the kinematic viscosity, which we define as (21)where the thermal and turbulent diffusivities are provided in Sects. 2.3, 2.4. The reference value is usually taken at the injection scale, l = L. Remarkably, the qualitative evolution of a very complex nonlinear dynamics can be approximately predicted via the dominant timescale ratio (or dimensionless number; see discussion in Sect. 4.1). For instance, different values of f and M can lead to the same P_{t}, hence to a similar qualitative dynamics and power spectrum of density perturbations (cf. Sects. 3.2 and 3.4).
3.1. Weak turbulence: M ~ 0.25
Fig. 2
Characteristic amplitude of δρ/ρ, , for the models with weak turbulence M ~ 0.25 and varying conduction, after reaching statistical steady state (~2 t_{eddy}) with the same level of continuous stirring. The driving is initiated only above 550 kpc. From top to bottom curve (black to bright red), the suppression of conduction is f = 0, 10^{3}, 10^{2}, 10^{1}, 1. First column of Table 1 summarizes the key properties. Strong conduction globally damps δ perturbations by a factor of 2−4, substantially steepening the spectrum and departing from the Kolmogorov slope of the noconduction run. Weak conduction is able to induce the steep decay only near the scale linked to P_{t} ~ 100. 
The first set of models implements a low level of stirring, with average massweighted 3D Mach number M ~ 0.25 (~370 km s^{1}). Observations and simulations suggest in fact that turbulence in the ICM typically remains subsonic (Nagai et al. 2007, 2013; Vazza et al. 2011; Gaspari et al. 2012b; Sanders & Fabian 2013). The turbulent energy is ~3.5 percent of the total thermal energy, E_{turb} ≃ 0.5 γ(γ − 1) M^{2}E_{th} ≃ 0.56 M^{2}E_{th}. Dissipational heating is thus negligible. In the current models, t_{eddy} = L/σ_{v} ~ 1.6 Gyr. We always analyze the system as soon as it establishes statistical steady state, i.e. after ~2 t_{eddy}.
In the purely hydrodynamic run (f = 0, i.e. no conduction or P_{t} → ∞) the driven stirring motions generate a turbulent cascade in the δρ/ρ power spectrum analogous to that of the turbulent velocities. The characteristic amplitude shows the typical injection peak at low k (l ~ 600 kpc), followed by the inertial range and the final steepening due to dissipation (Fig. 2, black line). The characteristic level of perturbations, given by the maximum of A(k)_{δ}, is 6.5 percent (6.7 percent using FFTs). An important result is that the inertial range of the density perturbations is remarkably similar to the Kolmogorov slope, A(k)_{δ} ∝ k^{− 1/3} (or E(k) ∝ k^{− 5/3}), slightly flattening towards the injection scale. Stratification has overall a secondary impact on δ (Sect. 4.1). Dissipation via numerical viscosity (mimicking Spitzer viscosity; Sect. 2.6) becomes substantial below 6 resolution elements. Notice that current observations are limited to scales ≳30 kpc due to Poisson noise and projection (Fig. 9), which are well resolved by the current runs. The dynamics is driven by key hydrodynamical instabilities, as KelvinHelmholtz (KH) and RayleighTaylor (RT), inducing the characteristic rolls, curls, and edges in both δ (Fig. 3) and SB_{x} maps (Fig. 4, top).
Overdense and underdense regions are associated with a mixture of isobaric and adiabatic perturbations, which we analyze in Sect. 4.2. With weak stirring the former dominates, while strong turbulence enhances the adiabatic mode (analogous to pressure waves). Unlike in the strongly conductive runs, the entropy gradient becomes progressively shallower, inducing a lower central density (30 percent) and higher temperature (10 percent). The transport of heat due to turbulent diffusion is ∝D_{turb}∇s (Eq. (11)). Therefore, turbulence seeds perturbations in density and temperature, while diffusing entropy. Conduction diffuses instead temperature and density fluctuations.
In the next experiment, we enable electron thermal conduction. We start analyzing the unmagnetized case (f = 1), though even in the unsuppressed regime, conduction can be limited by the saturated flux (Sect. 2.4). The overall dynamics and A(k)_{δ} is however not affected by saturation; the fraction of saturated zones is less than 10 percent (becoming ≲1 in the f = 0.1 run). The discrepancy between the purely turbulent and the conductive run is evident in the density perturbations (Fig. 2, red line), allowing to put critical constraints on the physical properties of the ICM. Three are the key modifications imparted by conduction. First, density perturbations are significantly damped over the whole range, by a factor of 4 on small scales to a factor of 2 on large scales, where the peak of perturbations reaches 3.5 percent. Second, the slope after the injection hump is considerably steeper than Kolmogorov, following a Burgerslike^{3} spectrum k^{− 1/2} (E(k) ∝ k^{2}). Third, there is no evident cutoff, meaning that conduction efficiently operates on all scales, while turbulence is not able to consistently regenerate perturbations. Considering the dimensionless turbulent Prandtl number (Eq. (21)), at the injection scale P_{t} ~ 1 (Fig. 8). Albeit thermal diffusion is a smallscale process (P_{t} ∝ l), it is ubiquitous and quick enough to efficiently stifle the full turbulent cascade in the whole cluster. The spectrum does not sharply decline to zero value, since turbulent regeneration is continuous. This is in analogy with the observed radio spectra, where relativistic electrons are regenerated by turbulence, preventing a dramatic cutoff.
The δρ/ρ maps of the strongly conductive runs (Fig. 3, bottom) clearly show the absence of significant density perturbations, especially on small scales. Since turbulent diffusion is severely inhibited, the cluster strongly retains the initial spherical symmetry and radial profiles, as indicated by the Xray surface brightness map in Fig. 4 (bottom row). Only the electron temperature is able to become quickly isothermal both locally and globally; due to the nonnegligible ion–electron equilibration delay, ~50 Myr (Sect. 2.5), the ion temperature has instead difficulty in becoming fully isothermal as T_{e}. The discrepancy is in the range (T_{e} − T_{i})/T_{i} ~ 1 − 15 percent, from the inner to external radial shells (in particular beyond r_{c}, where the decreasing densities increase the lagging; cf. Sect. 4.1). In the opposite regime, the purely turbulent run shows numerous filaments and depressions in density. In the SB_{x} image (top row), the perturbations are partially veiled by the line of sight integration; the ideal location to observe perturbations is at r ≳ r_{c} (see Appendix B).
Fig. 3
Midplane cuts of δρ/ρ for the models with M ~ 0.25. From top to bottom: f = 0, 10^{3}, 10^{2}, 10^{1} (the latter very similar to f = 1 run). The color coding is blue → white → red: −40% → 0% → 40%. 
We apply now a suppression factor of 10 on thermal conduction. The regime f ~ 0.1 − 0.3 is widely adopted in astrophysical studies (e.g. Narayan & Medvedev 2001; Voigt & Fabian 2004; Dennis & Chandran 2005; Ruszkowski & Oh 2010, 2011; see Sect. 2.4.1). As shown in Fig. 2, the characteristic amplitude is similar to the unsuppressed case. This is a key result, telling us that density perturbations are suppressed regardless of commonly adopted suppression factors (f ≳ 0.1). Observations could hence put strong limits on the suppression f, based on the steepness and/or normalization of A(k)_{δ} spectrum (Sect. 4). Compared with f = 1 model, the density fluctuations increase by a factor of 0.3 on large scales, while small perturbations have similar power. The absence of a dramatic decline is mainly due to the fact that on small scales the sound crossing time becomes greater than the conduction time (e.g. Fig. 8), hence the tiny bubbles do not have time to find a new pressure equilibrium. Besides global diffusion, strong conduction can thus promote minor stirring motions on small scales, preventing an abrupt decay of A(k)_{δ}. In this run, the spectrum slope in the inertial regime is steep, A(k)_{δ} ∝ k^{− 2/3}, significantly different from the noconduction run. Radial profiles and SB_{x} maps (Fig. 4, bottom) are very similar to the f = 1 model, retaining their initial spherical morphology. The P_{t} number is roughly 10 at the injection scale. Albeit turbulent regeneration starts to be more effective on large scales, the key P_{t} threshold appears to be an order of magnitude higher, as shown by the next model.
We suppress further the conductive flux by f = 10^{2}, a value advocated by several plasma physics theories (e.g. Rechester & Rosenbluth 1978; Chandran & Cowley 1998). The Prandtl number is 100 at the injection scale: turbulence can now restore part of the perturbations, though only near L (the normalization rises again to ~6 percent). This marked discrepancy between large and small scales induces a remarkably steep slope, A(k)_{δ} ∝ k^{− 4/5} (E(k) ∝ k^{2.5}), which should emerge in observed data in a clear way, if f ~ 10^{2} is the conductive regime of the ICM. The δρ/ρ map (Fig. 3, third panel) visualizes well the regeneration of turbulent eddies on large scales, while the smallscale flow remains considerably smooth, as corroborated by the SB_{x} map (Fig. 4, third row). Since this model shows a clear cutoff, it represents the cleanest case to retrieve the key threshold for the suppression of density perturbations, which we find to be P_{t} ~ 100. This is not a strict demarcation line, but rather a transition layer.
Only when conduction is substantially suppressed, f = 10^{3} (the typically lowest suppression factor adopted in theories), the turbulent cascade is significantly restored, generating the same peak and density spectrum down to ~L/2. Since thermal diffusion is too week, KelvinHelmholtz rolls and RayleighTaylor instabilities can develop again over a large range, defining the entire flow dynamics (Fig. 3) and perturbing the Xray surface brightness (Fig. 4, second row). Turbulent diffusion is able to efficiently mix the entropy profile, again lowering/increasing the central density/temperature (the discrepancy between T_{e} and T_{i} is now ≲1 percent; Sect. 4.1). Conduction can affect only the scales smaller than 100 kpc, creating a gentle exponential decrease in the logarithmic A(k)_{δ}. The suppression of δ reaches a factor of 2 near 30 kpc. Notice how conduction still dominates the diffusivity, overcoming (numerical) viscosity. When turbulent regeneration is efficient, it is not trivial to define an exact cutoff. Nevertheless, the threshold P_{t} ~ 100 (l ~ 100 kpc) appears a robust criterium: at that scale we see the beginning of a substantial decay of the density spectrum (changing slope to k^{− 1/2}).
Fig. 4
Xray surface brightness for the models with weak (M ~ 0.25; left), mild (M ~ 0.5; middle), and strong turbulence (M ~ 0.75; right). From top to bottom: f = 0, 10^{3}, 10^{2}, 10^{1} (f = 1 maps are similar to the latter images). The color coding is black − blue − pale green, in the range 10^{7} − 4 × 10^{5} erg s^{1} cm^{2}. Conduction prevents the development of the full turbulent cascade, especially if f ≥ 0.1. KH and RT rolls and filaments are thus suppressed, and the cluster retains the spherical, smooth shape. Strong turbulence is instead able to deform the cluster, flattening the entropy profile and inducing a fainter (more rarefied) core. Perturbations are best observed at r > r_{c}, or over the whole cluster if M > 0.5. 
3.2. Mild turbulence: M ~ 0.5
We now increase the level of turbulent motions by a factor of two, M ~ 0.5 (σ_{v} ~ 750 km s^{1}). Turbulent energy is thus ~14 percent of the thermal energy, still within the range retrieved by ICM observations and cosmological simulations. The characteristic eddy turnover time is t_{eddy} ~ 0.8 Gyr.
Figure 5 shows that the overall behavior of A(k)_{δ} is similar to the previous set of models, with differences laying in the details. The purely turbulent case (f = 0) forms the usual injection peak, with maximum at ~12 percent (Table 1), i.e. two times that of the previous run with half the turbulent velocity. Therefore, we infer that δρ/ρ ∝ M, or more precisely δρ/ρ ≃ 1/4 M (with l = L; the 0.25 factor is likely related to γ – see Sect. 4.1). This is a key result that will be further corroborated by the runs with stronger turbulence and smaller injection scale. In the nonconductive runs, the direct relation between the density and velocity power spectrum holds also on intermediate scales.
Fig. 5
Characteristic amplitude of δρ/ρ, , for the models with mild turbulence M ~ 0.5 and varying conduction. From top to bottom curve (black to bright red): f = 0, 10^{3}, 10^{2}, 10^{1}, 1. The level of density perturbations is linearly related to the Mach number, δ ~ A(k)_{δ,max} ≃ 1/4 M. Strong conduction globally damps perturbations by a factor of 2−3, steepening the spectrum towards ∝k^{− 1/2} (from slightly shallower than Kolmogorov in the f = 0 run). In the weakly conductive regime, the threshold for A(k)_{δ} decay is again P_{t} ~ 100. 
The inertial range of A(k)_{δ} in the purely turbulent run is slightly shallower than the Kolmogorov spectrum (∝k^{− 1/5}; cf. Kim & Ryu 2005). The decrease due to viscous dissipation occurs again within 15−20 kpc, while it is dominated by thermal diffusion in all the conductive runs. Given the significant level of turbulent diffusion, the entropy profile flattens twice more rapidly, doubling its central value to 700 keV cm^{2} in less than 2 Gyr (the bestfit profile has almost twice lower normalization, with 20 percent larger core radius). Turbulent dissipational heating still plays a secondary role, since the temperature at large radii is increased by just 15 percent (its timescale is ~4 t_{eddy}; Gaspari et al. 2013b). Overall, it is remarkable that the density perturbations follow in a fairly strict way the Kolmogorov velocities, even in the presence of a stratified atmosphere.
We proceed testing thermal conduction, starting from the highest suppression factor. As before, when f = 10^{3}, a gentle exponential cutoff develops. The Prandtl number is ~2000 at the injection scale (twice the previous models, P_{t} ∝ M), signaling that turbulent regeneration is extremely efficient, as highlighted by the filamentary SB_{x} image (Fig. 4, second column). The dynamics is dominated by KH and RT instabilities. We see a significant decline in density perturbations only at l ≲ 60 kpc, i.e. P_{t} ~ 100. This threshold appears again a robust indicator for the suppression of density inhomogeneities. When P_{t} < 100, the slope steepens from nearly Kolmogorov to Burgers spectrum (A(k)_{δ} ∝ k^{− 1/2}). The threshold and spectrum slope are independent of L (see Sect. 3.4), since shaped by the specific physical parameters of the plasma conductivity.
The effects of conduction can be best appreciated when f = 10^{2}. Using the key threshold P_{t} ~ 100, we predict a sharp decline around ~330 kpc, and indeed we see the k^{− 1/2} spectrum appearing on this large scale. The suppression of δ on small scales is a factor of 3, while near the injection scale turbulent regeneration is unhindered (A(k)_{δ} ~ 11 percent). The SB_{x} maps also reveal that only the large filaments, edges, and rolls are retained. Assuming a deep exposure, Xray imaging will thus help to assess the role of conduction and turbulence in the ICM.
Restoring the conductive flux to strong levels (f ~ 0.1 − 1, i.e. P_{t} ~ 20 − 2 at L) promotes a global suppression of δ, by a factor of ~2−3, from large to small scales (Fig. 5, red lines; the maxima of A(k)_{δ} show less separation compared with the models with weak turbulence). The slope is slightly shallower than Burgers spectrum, A(k) ∝ k^{− 4/9} (although turbulence has now twice more strength). As in Sect. 3.1, strong conduction can drive smallscale stirring, due to the fast conduction timescale relative to the gas/ion sound speed (Fig. 8). These models confirm that when P_{t} ≲ 20 (f ≳ 0.1) the density perturbations are strongly suppressed, producing a similar A(k)_{δ} and smooth/spherical SB_{x} maps (Fig. 4), regardless of the exact suppression factor.
A relevant result, which will be dissected in future work, is that conduction does not dramatically alter the power spectrum of the velocity field (nearly Kolmogorov), both in slope and normalization. Therefore, the quoted M or velocity dispersion (roughly v(k) at the injection scale) does not change between simulations with different suppression f.
3.3. Strong turbulence: M ~ 0.75
In the next experiments, we test the case with strong (still subsonic) turbulence, with three times higher turbulent velocities compared with the reference model (Sect. 3.1), M ~ 0.75 (σ_{v} ~ 1100 km s^{1}). Statistical steady state is reached fast, since t_{eddy} ~ 0.5 Gyr. The turbulent energy is ~0.31 the thermal energy. Approaching the transonic regime becomes progressively unrealistic: it has been shown both in observations and cosmological simulations that the turbulent pressure support in the ICM should remain in the range few −30 percent (e.g. Churazov et al. 2008; Lau et al. 2009; Vazza et al. 2009, 2011), the upper envelope defined by unrelaxed systems. Dissipational heating also starts to be relevant (∝M^{2}) increasing the overall temperature by ≳30 percent, though we are now interested in the relative variations.
The current set of simulations corroborates the previous key findings. According to our modeling, δρ/ρ ≃ 1/4 M, we predict a characteristic level of density perturbations of 18.75 percent (for the weakly conductive models, f < 0.1). Measuring the normalization of the A(k)_{δ} peak (f < 0.1), we retrieve a value of 18.8 percent (Table 1), an excellent match. Albeit the evolution is strongly chaotic and nonlinear, it is remarkable that we can convert the Mach number and density perturbations in a linear and simple way. In fact, the normalization of A(k)_{δ} is 2.9 (1.6) times that of the models with weak (mild) turbulence. Therefore, we expect real clusters to show density perturbations in the same range of the observed (low) Mach numbers. This could also explain why clusters do not show exaggerated clumpiness.
Fig. 6
Characteristic amplitude of δρ/ρ, , for the models with strong turbulence M ~ 0.75 and varying conduction: f = 0, 10^{3}, 10^{2}, 10^{1}, 1. These models corroborate the previous findings: the level of density perturbations is given by A(k)_{δ,max} ≃ 1/4 M (with l = L); conduction damps density fluctuations by at least a factor of 2, steepening the spectrum towards k^{− 1/2} (again from the Kolmogorov spectrum); the threshold for the A(k)_{δ} decay is P_{t} ~ 100. 
The other features of the density power spectrum, i.e. the slope and decay (Fig. 6), follow the analysis described in the previous two sections. We summarize the key points. The A(k)_{δ} slope of the f = 0 run is again slightly shallower than k^{− 1/3}, although not shallower than the M ~ 0.5 model, indicating that the Kolmogorov spectrum is indeed the “saturated” regime, at least for subsonic turbulence^{4}. Enabling conduction leads to the steeping of the characteristic amplitude towards the k^{− 1/2} regime, with the clearest manifestation in the f = 10^{2} case, where turbulence is able to regenerate perturbations near L (notice the short formation of the gentle exponential decline), while small scales are substantially damped. In the strong conductive regime (f ≳ 0.1) conduction inhibits the perturbations up to a factor of ~2. There is no drastic decline at small scales (A(k)_{δ} ∝ k^{− 4/9}) since conduction also promotes minor stirring, albeit the increased gas sound speed/entropy (due to turbulent heating and diffusion) alleviates this process. The commonly adopted suppression factors f ≥ 0.1 do not have a diverse impact on the density perturbations, which are damped over the whole range.
The turbulent Prandtl number is increased by a factor of 3 compared with the reference run, P_{t} ~ 3−3000 (f ~ 1−10^{3}) at the injection scale. This translates in a less marked – though not globally different – suppression between the purely turbulent and conductive case. More important, P_{t} ~ 100 corresponds to ~45 and 240 kpc for the f = 10^{3} and 10^{2} model, respectively (in the f ≳ 0.1 runs, the transition is well over L). Indeed, the substantial decline occurs around these scales. This threshold is thus a robust criterium to estimate the impact of conduction in a realistic cluster atmosphere.
The surface brightness maps (Fig. 4, third column) highlight the substantial level of perturbations, which are now clearly manifest even in the cluster core (especially when f ≳ 10^{2}). Such a configuration could be easily spotted in Xray images even by eye inspection. Notice how in the weakly suppressed runs, the KH and RT instabilities translate into extended filaments and curved edges/fronts in Xray brightness. The strong conductive flux instead smears out most of these features, preventing the ICM to depart from its spherical and smooth hydrostatic shape. This might be seen as an effective equation of state with a lower adiabatic index, implying a less reactive pressure during thermodynamic transformations (P ∝ ρ versus P ∝ ρ^{5/3}). Remarkably, the gas velocities remain ~1000 km s^{1} on large scales, with the v(k) spectrum not much affected by conduction (an indepth comparison will be presented in a companion work).
3.4. Small injection scale: M ~ 0.25, L/2
The final set of simulations includes testing turbulence and conduction diminishing the injection scale. To carry on a proper comparison, we set the new injection scale as L′ = L/2 ~ 300 kpc, but maintaining the turbulent Mach number of the reference run, M ~ 0.25 (via the scaling σ_{v} ~ (N L′ ϵ^{∗})^{1/3}; Sect. 2.3). The characteristic eddy turnover time is t_{eddy} ~ 0.8 Gyr, i.e. the same characteristic timescale of the previous models with M ~ 0.5 but larger injection scale (Sect. 3.2). The typical turbulent energy is 3.5 percent of E_{th} (as in Sect. 3.1).
The key result is that the overall A(k)_{δ} spectra are analogous to the models with M ~ 0.5 (Fig. 7). The characteristic Prandtl numbers, at the respective injection scales, are in fact the same. The major differences are twofold. First, the current spectra are truncated exactly at the smaller injection scale L′. Second, the normalization is defined only by the characteristic Mach number (and not P_{t}). In fact, using our previous A(k)_{δ} modeling and rescaling it with the new injection scale, we estimate A(k)_{δ,max} ≃ 1/4 M (L′/L)^{1/5} ≃ 0.054. The measured peak in the simulated run is 5.3 percent (Table 1), confirming well the prediction. Notice that the A(k)_{δ,max} fit is very weakly dependent on the injection scale, implying that the linear conversion between δ and M is fairly general (see Sect. 4.1).
Similar to previous M ~ 0.5 models, the strongly conductive runs (f ≳ 0.1) show density perturbations damped by a factor of 2−3, from large to small scales. The A(k)_{δ} slope follows also the previous trends, slightly shallower than Kolmogorov in the f = 0 run (∝k^{− 1/5}; higher D_{turb} means faster δ regeneration), while steepening towards the Burgerslike spectrum for the conductive runs. Using the threshold P_{t} ~ 100, we retrieve the region where A(k)_{δ} substantially decays. Interestingly, the estimated cutoff for the f = 10^{2} run is slightly larger than the injection scale (~330 kpc) and indeed we see the decrease in normalization (by 10 percent), instead of a progressive decline from the f = 0 baseline. In the f = 10^{3} run, the threshold is still below 100 kpc, allowing the exponential cutoff to partially develop, although entering quickly the regime of substantial suppression on few tens kpc. The hydrodynamics is similar to the M ~ 0.5 models, with KT and RT instabilities producing filaments and edges in SB_{x} (if f < 0.1), but now with maximum size ~100 kpc. On the other hand, strong conduction has the ability to preserve the smooth and spherical shape of the galaxy cluster (cf. last row in Fig. 4).
Fig. 7
Characteristic amplitude of δρ/ρ, , for the models with weak turbulence M ~ 0.25 and half the reference injection scale, L′ = L/2. From black to bright red line, the suppression of conduction is f = 0, 10^{3}, 10^{2}, 10^{1}, 1. The previous key features (e.g. slopes and P_{t} ~ 100 threshold) are valid also for this set of models. The system is dynamically similar to the runs with M ~ 0.5, since the P_{t} number (at injection scale) is the same (t_{turb} ~ 0.8 Gyr). Notice, however, that the A(k)_{δ} cascade is truncated at the new L′ and that the normalization is defined only by the Mach number (not P_{t}). 
4. Discussion
We now explain and discuss the key features of the models, and then constrain ICM conduction and turbulence in real Coma cluster through density perturbations.
4.1. Dominant timescales and simple predictions
The qualitative evolution of a dynamical system can be described via the timescale ratio of the physical processes involved, i.e. via characteristic dimensionless numbers. Let us consider as reference Fig. 8, showing the typical timescales for M ~ 0.25 and strong conduction (f = 1) run, within the core radius. The first key quantity is the (turbulent) Mach number, M = t_{cs}/t_{turb} (<1), which we found to be linearly related to the relative density perturbations δ. The scaling may be justified by simple physical arguments (see also Zhuravleva et al. 2013, and in prep.). Perturbations tend to settle back towards shells of similar entropy. Assuming an isobaric blob, its entropy varies as K ∝ ρ^{− γ} ∝ (1 + δ)^{− γ} – for an isothermal blob K ∝ (1 + δ)^{1 − γ}. Since the cluster entropy typically scales with the radius, the displacement is directly linked to the density contrast, Δr/r_{init} ≈ ΔK/K_{init} = 1 − (1 + δ)^{− γ} ≈ γ δ. In a stratified medium, the balancing force is F ∝ g δ, again directly related to the displacement. Since the specific kinetic energy is v^{2} ∝ F Δr ∝ (Δr)^{2}, we infer that M ∝ Δr ∝ δ, with the normalization shaped by the effective γ and the K(r) slope. In the limiting case of negligible buoyancy, density may be instead roughly treated as a passive scalar. In the classic theory of Obukhov & Corrsin (Warhaft 2000 for a review), the spectrum of passive scalars is expected to be directly related to v(k). The scaling δ ∝ M is thus fairly general.
Fig. 8
Typical timescales in the core of Coma (n_{e} and T_{e} are roughly constant) as function of l, for the run with M ~ 0.25 and f = 1: conduction (black; Eq. (16)), turbulence (magenta; Eq. (9)), sound (red; t_{cs} = l/c_{s}), electronion equilibration (cyan; Eq. (18)). The conduction timescale for the saturated zones corresponds to the upper black envelope, limited by the electron sound speed (≃43 times that of protons). The yellow line is the BruntVäisälä (buoyancy) timescale, as function of r, using an average g ~ 5 × 10^{9} cm s^{2} and shallow core entropy slope (~0.1); t_{BV} is similar to the freefall time, within a factor of a few. 
The second key quantity is the conduction timescale (Fig. 8), which is typically the shortest one (≲50 Myr). The black line shows that conduction dominates the dynamics over the whole range of scales, even in the saturated regime (upper black line). Increasing the level of suppression (f ≪ 1) induces the crossover of the conduction and turbulence (magenta) timescale. The predicted threshold for the turbulent regeneration of perturbations is P_{t} ≡ t_{cond}/t_{turb} ~ 100. If P_{t} ~ 100 occurs within the injection scale, the normalization of the spectrum does not decrease, and a decay appears on intermediate scales (Table 1). In this work, we used a turbulent diffusion coefficient c_{t} = 1, yet the actual value could be as low as 0.1−0.01 (Kim & Fabbiano 2003; Dennis & Chandran 2005), rescaling the threshold to P_{t} ~ 1.
Other timescales, linked to the relevant physics, are involved in the complex dynamics, although inducing secondary effects. The BruntVäisälä timescale tracks the impact of the restoring buoyant forces along the radial direction: t_{BV} = [γ (r/g) (d lnK/d lnr)^{1}] ^{1/2}. Buoyancy is thus defined by gravity g and the slope of the entropy profile^{5}. Within the core radius, ∇K is very shallow (and gravity roughly constant), hence t_{BV} > t_{turb}, or Froude number Fr > 1 (Fig. 8, yellow line), implying that the ICM is approaching neutral buoyancy and turbulent motions can easily retain isotropy. On large scales, the entropy profile follows the selfsimilar slope ∝r, as most clusters, but it is balanced by the decreasing gravity: buoyancy can again not efficiently inhibit radial motions. Substantially weaker stirring may progressively force the chaotic flow to follow tangential streamlines (Fr ≪ 1). Notice that the buoyancy timescale is essentially the freefall time, within a factor of a few: t_{ff} = (2r/g)^{1/2}. The ratio of the freefall time and the gas soundcrossing time (red line) is always >1, i.e. the cluster is in global hydrostatic equilibrium. As seen in Fig. 8, in the strongly conductive runs (f ≳ 0.1) the soundcrossing time may slightly fall below the conduction time, approaching the small scales. Albeit saturation prevents electrons to conduct faster than their thermal speed, t_{cs}/t_{cond} ≲ 1 warns that small fluctuations do not have time to readjust via pressure equilibrium. Progressively stronger conduction is thus not able to completely stifle A_{δ} at large k.
The electronion equilibration introduces another relevant timescale (Fig. 8, cyan). On small scales, t_{e − i} is usually larger than t_{cond} and comparable to the eddy turnover time. This means that the proton temperature can not become as quickly isothermal as T_{e} (lagging by 50−100 Myr), a feature aggravated by the continuously chaotic ambient. The discrepancy resides in the range 1−15 percent, from the inner to outer radial bins (for model M ~ 0.25, f = 1), decreasing to a maximum of 10, 3, and 1 percent, for f = 10^{1},10^{2}, and 10^{3}, respectively. Higher Mach number linearly increases the maximum discrepancy (occurring at r > r_{c}, where the medium is more rarefied). The estimate of  T_{e} − T_{i}  /T_{i} can be used to properly constrain the turbulent velocities from the observed spectral lines (e.g. Inogamov & Sunyaev 2003). In order to remove thermal line broadening (∝), it is assumed T_{i} = T_{e}, hence an effectively lower ion temperature would imply higher turbulent σ_{v}.
4.2. Modes of perturbations: isobaric, isothermal, isentropic
It would be tempting to characterize chaotic fluctuations with a single mode, as adiabatic (isentropic), isobaric, or isothermal. As experiment, we set up the cube with pure isobaric fluctuations (fixing Kolmogorov spectrum). The hydro run develops the turbulent cascade, preserving the isobaric fluctuations, which gradually decay in few t_{eddy} due to no forcing. However, as soon as conduction is enabled, the mode of perturbations suddenly change from isobaric to isothermal/adiabatic. The result is the generation of acoustic oscillations and, strikingly, a density power spectrum similar to that of the cosmic microwave background (CMB; Planck Collaboration 2013), although the two processes are inherently different – CMB evolution is driven by (compressive) gravity. The spectra observed in clusters do not present such an imprint, meaning that the actual type of fluctuations is a mixture of desynchronized modes between the adiabatic and isobaric extremes^{6}.
We analyze our models to understand the correlation between relative density (δ_{ρ}) and temperature fluctuations (δ_{T} = T_{e}/T_{b} − 1, where T_{b} is the underlying radial profile). The correlation is positive for adiabatic fluctuations (δ_{ρ}/δ_{T} = γ − 1; constant entropy), while negative for the isobaric mode (δ_{ρ}/δ_{T} = −1); the isothermal mode is the intermediate regime (δ_{ρ}/δ_{T} = ∞). Analyzing the linear regression coefficient of the δ_{T}δ_{ρ} diagram (512^{3} points) is not much meaningful, since there is too high dispersion. We use thus Pearson coefficient^{7} to assess the degree of positive and negative linear correlation. In the noconduction run with M ~ 0.25, the correlation is r_{P} ≃ − 0.86, implying that weak turbulence prefers isobaric fluctuations (adiabatic sound waves are negligible). Since stirring is subsonic, the fluctuations can not deviate much from hydrostatic equilibrium. When conduction is enabled, though partially suppressed, the scatter in the phase diagram progressively increases: fluctuations start filling also the adiabatic region, showing r_{p} = −0.82, − 0.78, − 0.65, for the runs with f = 10^{3},10^{2},10^{1}, respectively. With full Spitzer conduction (f = 1), the correlation drops to r_{p} = −0.30: the perturbations are now almost equally isobaric and adiabatic.
Increasing the strength of stirring (M ~ 0.75), rises the level of turbulent mixing, leading to shallower entropy gradients. In fact, even in the hydro run, the adiabatic and isobaric modes are balanced (r_{P} ≃ − 0.20). Conduction tries to push the correlation towards the isothermal regime, albeit there is no clear preference in the mixture of perturbations (r_{P} ~ − 0.1; f = 1 run). Smaller injection scales also induce a more efficient turbulent cascade, promoting mixing and leading to higher adiabatic fluctuations (the latter are also enhanced in supersonic^{8} turbulence). Overall, we suggest to exploit the correlation between δ_{T} and δ_{ρ} to assess the dominant physics in the ICM.
4.3. Comparison with real Coma cluster
In order to constrain the transport properties in the bulk of the ICM, we take as exemplary case Coma cluster, due to the availability of new deep observations (~500 ks; Churazov et al., in prep.). We refer to Churazov et al. (2012) for the extraction method of A_{δ}, accounting for Poisson noise, point sources, and deprojection. The hot ICM, as in Coma, is best suited for this analysis, due to the potential strong conductivity and negligible cooling or AGN heating. In cold, dense systems (T_{vir} < 5 keV) conduction is at least 100−1000 times weaker, thus its impact is expected to be marginal (Dolag et al. 2004; Smith et al. 2013), even if unsuppressed, and very difficult to isolate (the analysis of other systems will be provided in a subsequent work).
The previous simulations provided a clear picture of the interplay between conduction and turbulence in the hot ICM. We can thus formulate a simple, general model linking conduction and turbulence to the level of density fluctuations at a given physical scale l (=k^{1}): (22)The term in brackets represents the normalization of δ which must be rescaled to the current injection scale L_{inj}, using the slope given by the purely hydrodynamic cascade, α_{h} (≈0.2−0.3; Table 1). As shown before, the characteristic spectrum slope α_{c} is dictated by the level of conduction, and can be much steeper than α_{h} (α_{c} ≳ 0.45, for f ≳ 10^{2}). A steep slope combined with globally smooth residual/SB_{x} maps indicates that the normalization should be reduced by 1.3−2, due to very strong conductivity (f ≳ 0.1; Table 1). This prescription is sufficient to assess the physical state of the hot ICM. For a more precise modeling, we suggest to insert two (exponential) cutoffs linked to the injection, exp[− ^{(}l/a_{1}L^{)}^{η1}], and dissipation scale, exp[− (a_{2}Δx/l)^{η2}]; typical parameters fitting our simulated spectra are η_{1} = 4, η_{2} = 3/2, a_{1} = 2, a_{2} = 3.
In Fig. 9, we present the characteristic amplitude of density perturbations observed in real Coma cluster (within statistical uncertainties; cyan envelope), compared with the prediction based on our modeling (black line). It is clear that the observed spectrum is shallow. A powerlaw with slope α_{c} ≃ 0.36, slightly steeper than Kolmogorov (0.33), fits well the inertial range of the spectrum. Coma observations of pressure perturbations further corroborate this trend (Schuecker et al. 2004). Any plasma with significant conduction would produce instead a steeper spectrum. We can therefore exclude a strong or mild conductive state of the ICM, f ≥ 10^{2}. The inertial range develops unimpeded down to tens kpc, indicating that turbulent regeneration is substantial and that the Prandtl number at the injection scale is high (P_{t} ≳ 2000). The fact that we do not see a sharp cutoff, but only a gentle decay, indicates that the suppression factor is very low, f ≃ 10^{3}, although not f = 0 (the nonconductive cascade would be significantly higher below 60 kpc; P_{t} ~ 100). In addition, the observed δ map (Churazov et al. 2012) is not smooth, revealing a variegated morphology of local features, reminiscent of Fig. 3 (top two panels). The relation between the velocity and density spectrum can finally break any ambiguity of the proper regime (Gaspari et al., in prep.).
As reviewed in Sects. 1 and 2.4.1, the level of suppression in turbulent plasmas is a matter of debate. It strongly depends on the topology of magnetic field lines and the relation between the electrons mean free path and the field correlation length, l_{B} (e.g. Rechester & Rosenbluth 1978; Chandran & Cowley 1998; Malyshkin & Kulsrud 2001; Narayan & Medvedev 2001; Chandran & Maron 2004). For example, the latter authors suggested a suppression f ~ 0.1 if λ_{e} ≪ l_{B}, while in the opposite regime magnetic mirrors can further suppress heat transport. On top of this, kinetic plasma instabilities can substantially lower f. If the field is more ordered, the conductivity across the field can be heavily suppressed, f ≪ 1, as found by observational and theoretical analyses of sharp features in the ICM, e.g. cold fronts (Ettori & Fabian 2000; Roediger et al. 2013; ZuHone et al. 2013) and linear filaments (Forman et al. 2007; Sanders et al. 2013). Interestingly, the low value of f could explain the survival of positive T gradients in the center of many clusters since z ~ 1.2 (McDonald et al. 2013). It is also clear that such a low conductivity can not provide sufficient heating for quenching radiative cooling, leaving AGN feedback as the main solution for the cooling flow problem (Sect. 1).
Fig. 9
Cyan envelope: characteristic amplitude of δρ/ρ in real Coma, extracted from deep Chandra observations (Churazov et al., in prep.; ~500 ks), within the statistical uncertainties. Black line: A_{δ}(k) prediction based on our modeling. The state of the ICM is best matched by our simulations with highly suppressed conduction, f ≃ 10^{3}, since the spectrum shows a shallow slope, α_{c} ≃ 0.36, slightly steeper than Kolmogorov. The normalization (10 percent) indicates that the level of turbulence is M ≃ 0.45. The spectrum also suggests an injection scale of ~250 kpc, although the precise value is uncertain due to the selection of the unperturbed model. 
The A_{δ}(k) peak, i.e. the normalization, is ≃10 percent, which can be achieved by a level of turbulence with M ≃ 0.45. The ratio of turbulent to thermal energy is thus 11 percent, concordant with observations and cosmological simulations of clusters (Sect. 2.3). The spectrum suggests an injection scale ~250 kpc, similar to the estimated impact parameter of merging clusters (Sarazin 2002).
We warn that the observed low frequencies have significant uncertainty, due to the stochastic nature of perturbations and modeling errors, as the deprojection and the underlying profile removal. For instance, Churazov et al. (2012) tested asymmetric models, finding that nonsphericity can vary the deprojected amplitude above 300 kpc up to a factor of 2, while smaller scales remain unaffected. The current simulations also do not model mergers or external accretion. Density inhomogeneities due to the substructures in the potential and accreted filaments may come on top of the largescale perturbations directly related to the turbulent cascade (Churazov et al. 2012; Sanders et al. 2013). Multiple injection scales may also flatten the spectrum, although they require finetuning to be masked as a single spectral profile. While the turbulent cascade can be formed by a variety of drivers, its appearance, modeled in the present work, is nevertheless expected to be largely universal. It is no coincidence that the cascade in cosmological simulations is nearly Kolmogorov (e.g. Vazza et al. 2011), similar to our hydro runs. We remark that strong conduction would be clear in the spectrum, regardless of the abovementioned mechanisms, uniquely stifling the cascade.
On very small scales, the observed spectrum is limited by Poisson noise, overcoming the fluctuations below ~30 kpc. The range 30−300 kpc is yet sufficient to determine the state of the ICM, since the effects of diffusion would clearly emerge within this range (Sect. 3). A steep slope can indeed not be produced by changing M or L_{inj}, but only increasing f. The SB_{x} map can further clarify the proper regime. Notice that if the precise L_{inj} is larger, the estimate of M should not vary significantly, since the more extended cascade provides a higher A_{δ}(k) peak. With even deeper data, the A_{δ}(k) slope could be traced down to smaller scales, further constraining the role of conduction (notice that low resolution smoothens the δ field). Future deeper Xray observations might further improve these values, perhaps combining several clusters, and providing a highprecision spectrum of ICM perturbations.
In passing, we note that the δ prescription (Eq. (22)) can be very useful for those (subgrid) simulations, semianalytic or analytic works, which intend to model, ab initio and in a simple way, the interplay of conduction, turbulence, and inhomogeneities in the ICM. For instance, the above equation can be used either to set the level of density (and T; Sect. 4.2) perturbations from the Mach number (particularly useful for clumpiness studies), or to impose the effective Mach number from the amplitude of density fluctuations. The current method can be readily applied to other areas of astrophysics, as the interstellar or intergalactic medium. Interestingly, also the ISM shows a nearly Kolmogorov slope over five orders of magnitude (Armstrong et al. 1995), though steepening in dwarf galaxies (Dutta et al. 2009); conduction seems thus substantially suppressed in different astrophysical atmospheres.
5. Conclusions
Using 3D highresolution hydrodynamic simulations, we thoroughly tested the physics of thermal conduction (f = 0 → 1) and turbulence (M = 0.25 → 0.75) in the hot ICM, tracking both electrons and ions (T_{e} − T_{i}  /T_{i} ≲ 0.15). We showed how to exploit the power spectrum of the relative gas density perturbations δ = δρ/ρ (normalization, slope, decay), in order to accurately constrain the effective conductive and turbulent state of the ICM. We chose Coma cluster as reference laboratory. The main results are as follows.

The normalization of the characteristic amplitude spectrum, , determines the strength of turbulence, and vice versa. The peak amplitude is linearly related to the average 3D Mach number as A(k)_{δ,max} = c M, where c ≃ 0.25 for injection scale L_{inj} ≃ 500 kpc. In the hydro run, the spectrum of δ tracks that of velocities even on intermediate scales. After rescaling, the steady spectra are remarkably similar despite different M or L_{inj}, allowing to build a general theory of ICM perturbations.

The slope of A_{δ}(k) determines the level of diffusion dominated by conduction. In a nonconductive ICM (f = 0), subsonic stirring motions generate a nearly Kolmogorov cascade, E_{δ}(k) ∝ k^{− 5/3}. Similar to velocities, the inertial range of density perturbations peaks at the injection scale, and decays below 10 kpc due to “viscous” dissipation. Increasing the level of conduction, with magnetic suppression f = 10^{3} → 1, progressively steepens the spectrum slope towards the Burgerslike regime (E_{δ}(k) ∝ k^{2}; Table 1), a feature that would be manifest in observations. The velocity spectrum is not much affected by conduction. Even with f = 1, the δ amplitude does not drop to zero, since strong conduction also promotes fluctuations. The A_{δ}(k) slope is only weakly dependent on M, becoming slightly shallower above M ≳ 0.5.

The dominant dimensionless number or timescale ratio, shaping the flow dynamic similarity (Sect. 4.1), is the turbulent Prandtl number: P_{t} = D_{turb}/D_{cond} = t_{cond}/t_{turb}. The threshold P_{t} ≲ 100 approximately indicates where A_{δ}(k) has a significant decay due to diffusion. The transition is very gentle for strong suppression of conduction, f ≲ 10^{3}, becoming a sharp decay – though not a cutoff – in the opposite regime. If P_{t} ~ 100 occurs above the injection scale, density perturbations are inhibited over the whole range of scales, inducing a decrease in normalization up to a factor of ~2 (on small scales the suppression can reach a factor of 4). This state occurs only with strong conductivity, f ≥ 0.1 and would be pinpointed by the SB_{x} or residual δ images, in which KH/RT rolls and filaments are washed out, preserving the smooth and spherical shape of the cluster.

Realistic perturbations, in a stratified system, are characterized by a mixture of modes, shaped by the dominant physics. Weak turbulence drives higher isobaric perturbations (δ_{ρ}/δ_{T} = −1); strong turbulence enhances the adiabatic modes (δ_{ρ}/δ_{T} = γ − 1), while increasing conduction forces both modes towards the intermediate isothermal regime. Although density statistics is much better constrained by observations, unveiling temperature perturbations could substantially advance our knowledge of the ICM physics.

Based on our experiments, we provided a general analytic model to constrain density perturbations, conduction, and turbulence in the bulk of the ICM: δ(l) ≃ [0.25 M (L_{inj}/L_{500})^{αh}] (l/L_{inj})^{αc} (see Sect. 4.3). We apply it to new very deep Chandra observations of Coma. The observed A_{δ}(k) spectrum is best consistent with strongly suppressed conduction, f ≃ 10^{3}, and mild subsonic turbulence, M ≃ 0.45 (L_{inj} ~ 250 kpc). The latter (E_{turb} ≃ 0.11 E_{th}) is in agreement with cosmological simulations of clusters formation and linebroadening observations. The low conductivity corroborates the survival of local sharp features (cold fronts, filaments, bubbles), and indicates that cooling flows may not be balanced by conduction, leaving AGN feedback as the main driver of heating.
The increasing quality and sample size of future Xray data will provide a key opportunity to exploit this new spectral modeling, and hopefully to open the path to highprecision physics of the ICM, as has been done for the CMB.
Via simple dimensional analysis (Ruszkowski & Oh 2010), σ_{v} ∝ (N L ϵ^{∗})^{1/3}; the number of modes is typically N < 1000.
In reality, only kinetic 3D simulations, solving the Vlasov equations, can exactly determine the effective f, tightly linked to magnetic instabilities, as firhose, mirror, etc. (Schekochihin & Cowley 2007), which is out of reach for the current computing power.
In the supersonic regime (Federrath 2013), the development of filamentary structures via strong shocks may further flatten the δ power spectrum (e.g. Kim & Ryu 2005).
In the case of anisotropic conduction, buoyancy is constrained by ∇T (Sharma et al. 2009), which is usually shallower than ∇K.
Cosmological runs also show mixed perturbations (Zhuravleva et al. 2013), further corroborating the realism of our turbulence driving.
Acknowledgments
The FLASH code was in part developed by the DOE NNSAASC OASCR Flash center at the University of Chicago. We acknowledge the MPA, RZG, and CLS center for the availability of highperformance computing resources. M.G. is grateful for the financial support provided by the Max Planck Fellowship. E.C. acknowledges support by the Research Program OFN17 of the Division of Physics, Russian Academy of Science. We thank R. Sunyaev, A. Schekochihin, F. Vazza, and the anonymous referee for helpful comments.
References
 Arévalo, P., Churazov, E., Zhuravleva, I., HernándezMonteagudo, C., & Revnivtsev, M. 2012, MNRAS, 426, 1793 [NASA ADS] [CrossRef] [Google Scholar]
 Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A. 1986, ApJ, 304, 787 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & McKee, C. F. 1982, ApJ, 252, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Borgani, S., & Kravtsov, A. 2011, Adv. Sci. Lett., 4, 204 [CrossRef] [Google Scholar]
 Brighenti, F., & Mathews, W. G. 2003, ApJ, 587, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Chandran, B. D. G., & Cowley, S. C. 1998, Phys. Rev. Lett., 80, 3077 [NASA ADS] [CrossRef] [Google Scholar]
 Chandran, B. D. G., & Maron, J. L. 2004, ApJ, 602, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Churazov, E., Forman, W., Jones, C., & Böhringer, H. 2000, A&A, 356, 788 [NASA ADS] [Google Scholar]
 Churazov, E., Brüggen, M., Kaiser, C. R., Böhringer, H., & Forman, W. 2001, ApJ, 554, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Churazov, E., Forman, W., Vikhlinin, A., et al. 2008, MNRAS, 388, 1062 [NASA ADS] [CrossRef] [Google Scholar]
 Churazov, E., Vikhlinin, A., Zhuravleva, I., et al. 2012, MNRAS, 421, 1123 [NASA ADS] [CrossRef] [Google Scholar]
 Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135 [NASA ADS] [CrossRef] [Google Scholar]
 de Plaa, J., Zhuravleva, I., Werner, N., et al. 2012, A&A, 539, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dennis, T. J., & Chandran, B. D. G. 2005, ApJ, 622, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Dolag, K., Jubelgas, M., Springel, V., Borgani, S., & Rasia, E. 2004, ApJ, 606, L97 [NASA ADS] [CrossRef] [Google Scholar]
 Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Dutta, P., Begum, A., Bharadwaj, S., & Chengalur, J. N. 2009, MNRAS, 398, 887 [NASA ADS] [CrossRef] [Google Scholar]
 Eswaran, V., & Pope, S. B. 1988, Comput. Fluids, 16, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Ettori, S., & Fabian, A. C. 2000, MNRAS, 317, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C. 2013, MNRAS, in press [arXiv:1306.3989] [Google Scholar]
 Fisher, R. T., Kadanoff, L. P., Lamb, D. Q., et al. 2008, IBM J. Res. Dev., 52, 127 [Google Scholar]
 Forman, W., Jones, C., Churazov, E., et al. 2007, ApJ, 665, 1057 [NASA ADS] [CrossRef] [Google Scholar]
 Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., Melioli, C., Brighenti, F., & D’Ercole, A. 2011a, MNRAS, 411, 349 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., Brighenti, F., D’Ercole, A., & Melioli, C. 2011b, MNRAS, 415, 1549 [Google Scholar]
 Gaspari, M., Ruszkowski, M., & Sharma, P. 2012a, ApJ, 746, 94 [Google Scholar]
 Gaspari, M., Brighenti, F., & Temi, P. 2012b, MNRAS, 424, 190 [NASA ADS] [CrossRef] [Google Scholar]
 Gaspari, M., Brighenti, F., & Ruszkowski, M. 2013a, Astron. Nachr., 334, 394 [NASA ADS] [Google Scholar]
 Gaspari, M., Ruszkowski, M., & Oh, S. P. 2013b, MNRAS, 432, 3401 [NASA ADS] [CrossRef] [Google Scholar]
 Gazol, A., & Kim, J. 2010, ApJ, 723, 482 [NASA ADS] [CrossRef] [Google Scholar]
 Huba, J. D. 2009, NRL Plasma Formulary (NRL/PU/6790−11551) [Google Scholar]
 Inogamov, N. A., & Sunyaev, R. A. 2003, Astron. Lett., 29, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, D.W., & Fabbiano, G. 2003, ApJ, 586, 826 [Google Scholar]
 Kim, J., & Ryu, D. 2005, ApJ, 630, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, K.T., Kronberg, P. P., Dewdney, P. E., & Landecker, T. L. 1990, ApJ, 355, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Kissmann, R., Kleimann, J., Fichtner, H., & Grauer, R. 2008, MNRAS, 391, 1577 [Google Scholar]
 Komarov, S., Churazov, E., & Schekochihin, A. 2013 [arXiv:1304.1857] [Google Scholar]
 Kowal, G., Lazarian, A., & Beresnyak, A. 2007, ApJ, 658, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Lau, E. T., Kravtsov, A. V., & Nagai, D. 2009, ApJ, 705, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, D., & Deane, A. E. 2009, J. Comput. Phys., 228, 952 [NASA ADS] [CrossRef] [Google Scholar]
 Malyshkin, L., & Kulsrud, R. 2001, ApJ, 549, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McDonald, M., Benson, B. A., Vikhlinin, A., et al. 2013, ApJ, 774, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, D., Vikhlinin, A., & Kravtsov, A. V. 2007, ApJ, 655, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, D., Lau, E. T., Avestruz, C., Nelson, K., & Rudd, D. H. 2013, ApJ, 777, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., & Medvedev, M. V. 2001, ApJ, 562, L129 [NASA ADS] [CrossRef] [Google Scholar]
 Norman, M. L., & Bryan, G. L. 1999, in The Radio Galaxy Messier 87, eds. H.J. Röser & K. Meisenheimer (Berlin: Springer Verlag), Lect. Notes Phys., 530, 106 [Google Scholar]
 Parrish, I. J., Quataert, E., & Sharma, P. 2009, ApJ, 703, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration 2013 [arXiv:1303.5075] [Google Scholar]
 Rechester, A. B., & Rosenbluth, M. N. 1978, Phys. Rev. Lett., 40, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, C. S., McKernan, B., Fabian, A. C., Stone, J. M., & Vernaleo, J. C. 2005, MNRAS, 357, 242 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Roediger, E., Kraft, R. P., Forman, W. R., Nulsen, P. E. J., & Churazov, E. 2013, ApJ, 764, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., & Begelman, M. C. 2002, ApJ, 581, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., & Oh, S. P. 2010, ApJ, 713, 1332 [NASA ADS] [CrossRef] [Google Scholar]
 Ruszkowski, M., & Oh, S. P. 2011, MNRAS, 414, 1493 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. S., & Fabian, A. C. 2012, MNRAS, 421, 726 [NASA ADS] [Google Scholar]
 Sanders, J. S., & Fabian, A. C. 2013, MNRAS, 429, 2727 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. S., Fabian, A. C., Churazov, E., et al. 2013, Science, 341, 6152 [CrossRef] [Google Scholar]
 Sarazin, C. L. 2002, in Merging Processes in Galaxy Clusters, eds. L. Feretti, I. M. Gioia, & G. Giovannini, Astrophys. Space Sci. Lib., 272, 1 [Google Scholar]
 Schekochihin, A. A., & Cowley, S. C. 2007, in Turbulence and Magnetic Fields in Astrophysical Plasmas, ed. S. Molokov, R. Moreau, & H. K. Moffatt (Springer), 85 [Google Scholar]
 Schuecker, P., Finoguenov, A., Miniati, F., Böhringer, H., & Briel, U. G. 2004, A&A, 426, 387 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sharma, P., Chandran, B. D. G., Quataert, E., & Parrish, I. J. 2009, in AIP Conf. Ser. 1201, eds. S. Heinz, & E. Wilcots, 363 [Google Scholar]
 Smith, B. D., O’Shea, B. W., Voit, G. M., Ventimiglia, D., & Skillman, S. W. 2013, ApJ, submitted [arXiv:1306.5748] [Google Scholar]
 Spitzer, L. 1962, Physics of Fully Ionized Gases (New York: Interscience, Wiley) [Google Scholar]
 Vazza, F., Brunetti, G., Kritsuk, A., et al. 2009, A&A, 504, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vazza, F., Brunetti, G., Gheller, C., Brunino, R., & Brüggen, M. 2011, A&A, 529, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691 [NASA ADS] [CrossRef] [Google Scholar]
 Voigt, L. M., & Fabian, A. C. 2004, MNRAS, 347, 1130 [NASA ADS] [CrossRef] [Google Scholar]
 Warhaft, Z. 2000, Ann. Rev. Fluid Mech., 32, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Zakamska, N. L., & Narayan, R. 2003, ApJ, 582, 162 [NASA ADS] [CrossRef] [Google Scholar]
 Zhuravleva, I., Churazov, E., Kravtsov, A., et al. 2013, MNRAS, 428, 3274 [NASA ADS] [CrossRef] [Google Scholar]
 ZuHone, J. A., Markevitch, M., Ruszkowski, M., & Lee, D. 2013, ApJ, 762, 69 [NASA ADS] [CrossRef] [Google Scholar]
Online material
Appendix A: Mexican Hat versus FFT spectrum
Since the data cube is non periodic, computing the power spectrum via FFTs is in principle inconsistent. As a consequence, the unresolved large scale power can leak into the available frequency range, distorting the spectrum. We use thus a modified Δvariance method, known as “Mexican Hat” filtering (MH; cf. Arévalo et al. 2012). For each spatial scale σ, the method consists of three steps:
 1.
the realspace cube C is convolved with two Gaussian filters having slightly different smoothing lengths: and , where ϵ ≪ 1;
 2.
the difference of the two cubes is computed, resulting in a cube dominated by the fluctuations at scales ≈σ (the difference of two Gaussian filters is simply the Mexican Hat filter, F(x) ∝ ϵ [1 − x^{2}/σ^{2}] exp [ − x^{2}/2σ^{2}], characterized by a positive core and negative wings);
 3.
the variance V_{σ} of the previous cube is calculated and recast into the estimate of the power, knowing that
Fig. A.1
Comparison of the characteristic amplitude spectra (for the run with M ~ 0.5 and f = 10^{2}), computed with two different methods: Mexican Hat filtering (black) and fast Fourier transforms (blue). The retrieved spectrum is consistent in both cases, without major differences. 
In Fig. A.1, we show the comparison between the MH and FFT spectrum, for the run with M ~ 0.5 and f ~ 10^{2}. In our study, there is no dramatic difference between the two methods. The slope in the inertial range is almost identical. At very small scales, the FFT spectrum produces a characteristic hook, in part due to the numerical noise near the maximum resolution, but also due to the contamination of jumps at the nonperiodic boundaries. The MH spectrum shows instead a gentle decline. In the opposite regime, the MH filter tends to smooth the scales greater than the injection scale, while the FFT spectrum shows a steeper decrease. The FFT peak is slightly higher, typically by 2−3 percent, likely affected by the nonperiodic box. Progressively trimming the box increases the relative normalization of the FFT spectrum, even by 20 percent, while distorting the lowfrequency slope; the MH spectrum is instead unaltered.
Appendix B: βprofile in Fourier space
Fig. B.1
Analytic 1D power spectra: βprofile (red), Kolmogorov noise (blue), βprofile perturbed by the noise (black; ). The spectrum is normalized to the value at k_{0} = 1/L = 0.01 (dimensionless units; 2π is dropped for clarity). The core radius is r_{c} = 20, i.e. L/5. The relative amplitude of the noise is ~10 percent. The noise clearly emerges beyond the core radius (k > 0.05), regardless of largescale structures. 
We present here the analytic conversion of the βprofile to Fourier space, and its interplay with a powerlaw noise. Using the notation , the Fourier transform of the βprofile (Eq. (1)) results to be (B.1)where ξ ≡ 3β/2, K_{1/2 − ξ} is a modified Bessel function of the second kind, and Γ is the Gamma function. The (1D) power spectrum is as usual retrieved as . Assuming β = 2/3 ≃ 0.66 (a typical value for galaxy clusters), the power spectrum of the βprofile reduces to (B.2)The previous equation strikes for its simplicity, and can be readily used in semianalytic studies. Changing β in the range 0.5−1 does not significantly alter P(k), hence Eq. (B.2) is an excellent approximation for the majority of clusters (Fig. B.1, red line). A remarkable feature is that the transition from real to Fourier space does not dramatically deform the profile, in tight analogy with Gaussian functions (∝exp [ − k^{2}]). The spectrum is dominated by the power on large scales, with the core radius playing a crucial role; a progressively rising r_{c} leads to an increase in both the normalization and steepness of the spectrum.
For our study, it is useful to analyze the superposition of the βprofile and a powerlaw Kolmogorov noise (with 1D power ∝k^{− 5/3}), n_{p} = n_{β} (1 + δ). Using the convolution theorem, the power spectrum of the perturbed density profile is given by . The cross terms cancel out since the
δ field is random and the phases are uncorrelated. In Fig. B.1, we show three power spectra: βprofile (red), noise with ~10 percent relative amplitude (blue), and the superposition of both (black). Beyond the core radius (k ≳ 0.05), the noise clearly starts to dominate. It is thus not essential to remove the underlying profile or largescale coherent structures, in order to unveil density perturbations, especially with substantial turbulence.
All Tables
Key properties of A_{δ}(k) for the simulated models: normalization (maximum), slope, and scale of significant decay.
All Figures
Fig. 1
Initial conditions for our reference hot galaxy cluster, Coma. From top to bottom panel: electron temperature, number density and entropy (). The radius is normalized to r_{500} ≃ 1.4 Mpc. We extracted the red data points from recent deep XMM observations, reaching ~0.5 r_{500} (Lyskova et al., in prep.). 

In the text 
Fig. 2
Characteristic amplitude of δρ/ρ, , for the models with weak turbulence M ~ 0.25 and varying conduction, after reaching statistical steady state (~2 t_{eddy}) with the same level of continuous stirring. The driving is initiated only above 550 kpc. From top to bottom curve (black to bright red), the suppression of conduction is f = 0, 10^{3}, 10^{2}, 10^{1}, 1. First column of Table 1 summarizes the key properties. Strong conduction globally damps δ perturbations by a factor of 2−4, substantially steepening the spectrum and departing from the Kolmogorov slope of the noconduction run. Weak conduction is able to induce the steep decay only near the scale linked to P_{t} ~ 100. 

In the text 
Fig. 3
Midplane cuts of δρ/ρ for the models with M ~ 0.25. From top to bottom: f = 0, 10^{3}, 10^{2}, 10^{1} (the latter very similar to f = 1 run). The color coding is blue → white → red: −40% → 0% → 40%. 

In the text 
Fig. 4
Xray surface brightness for the models with weak (M ~ 0.25; left), mild (M ~ 0.5; middle), and strong turbulence (M ~ 0.75; right). From top to bottom: f = 0, 10^{3}, 10^{2}, 10^{1} (f = 1 maps are similar to the latter images). The color coding is black − blue − pale green, in the range 10^{7} − 4 × 10^{5} erg s^{1} cm^{2}. Conduction prevents the development of the full turbulent cascade, especially if f ≥ 0.1. KH and RT rolls and filaments are thus suppressed, and the cluster retains the spherical, smooth shape. Strong turbulence is instead able to deform the cluster, flattening the entropy profile and inducing a fainter (more rarefied) core. Perturbations are best observed at r > r_{c}, or over the whole cluster if M > 0.5. 

In the text 
Fig. 5
Characteristic amplitude of δρ/ρ, , for the models with mild turbulence M ~ 0.5 and varying conduction. From top to bottom curve (black to bright red): f = 0, 10^{3}, 10^{2}, 10^{1}, 1. The level of density perturbations is linearly related to the Mach number, δ ~ A(k)_{δ,max} ≃ 1/4 M. Strong conduction globally damps perturbations by a factor of 2−3, steepening the spectrum towards ∝k^{− 1/2} (from slightly shallower than Kolmogorov in the f = 0 run). In the weakly conductive regime, the threshold for A(k)_{δ} decay is again P_{t} ~ 100. 

In the text 
Fig. 6
Characteristic amplitude of δρ/ρ, , for the models with strong turbulence M ~ 0.75 and varying conduction: f = 0, 10^{3}, 10^{2}, 10^{1}, 1. These models corroborate the previous findings: the level of density perturbations is given by A(k)_{δ,max} ≃ 1/4 M (with l = L); conduction damps density fluctuations by at least a factor of 2, steepening the spectrum towards k^{− 1/2} (again from the Kolmogorov spectrum); the threshold for the A(k)_{δ} decay is P_{t} ~ 100. 

In the text 
Fig. 7
Characteristic amplitude of δρ/ρ, , for the models with weak turbulence M ~ 0.25 and half the reference injection scale, L′ = L/2. From black to bright red line, the suppression of conduction is f = 0, 10^{3}, 10^{2}, 10^{1}, 1. The previous key features (e.g. slopes and P_{t} ~ 100 threshold) are valid also for this set of models. The system is dynamically similar to the runs with M ~ 0.5, since the P_{t} number (at injection scale) is the same (t_{turb} ~ 0.8 Gyr). Notice, however, that the A(k)_{δ} cascade is truncated at the new L′ and that the normalization is defined only by the Mach number (not P_{t}). 

In the text 
Fig. 8
Typical timescales in the core of Coma (n_{e} and T_{e} are roughly constant) as function of l, for the run with M ~ 0.25 and f = 1: conduction (black; Eq. (16)), turbulence (magenta; Eq. (9)), sound (red; t_{cs} = l/c_{s}), electronion equilibration (cyan; Eq. (18)). The conduction timescale for the saturated zones corresponds to the upper black envelope, limited by the electron sound speed (≃43 times that of protons). The yellow line is the BruntVäisälä (buoyancy) timescale, as function of r, using an average g ~ 5 × 10^{9} cm s^{2} and shallow core entropy slope (~0.1); t_{BV} is similar to the freefall time, within a factor of a few. 

In the text 
Fig. 9
Cyan envelope: characteristic amplitude of δρ/ρ in real Coma, extracted from deep Chandra observations (Churazov et al., in prep.; ~500 ks), within the statistical uncertainties. Black line: A_{δ}(k) prediction based on our modeling. The state of the ICM is best matched by our simulations with highly suppressed conduction, f ≃ 10^{3}, since the spectrum shows a shallow slope, α_{c} ≃ 0.36, slightly steeper than Kolmogorov. The normalization (10 percent) indicates that the level of turbulence is M ≃ 0.45. The spectrum also suggests an injection scale of ~250 kpc, although the precise value is uncertain due to the selection of the unperturbed model. 

In the text 
Fig. A.1
Comparison of the characteristic amplitude spectra (for the run with M ~ 0.5 and f = 10^{2}), computed with two different methods: Mexican Hat filtering (black) and fast Fourier transforms (blue). The retrieved spectrum is consistent in both cases, without major differences. 

In the text 
Fig. B.1
Analytic 1D power spectra: βprofile (red), Kolmogorov noise (blue), βprofile perturbed by the noise (black; ). The spectrum is normalized to the value at k_{0} = 1/L = 0.01 (dimensionless units; 2π is dropped for clarity). The core radius is r_{c} = 20, i.e. L/5. The relative amplitude of the noise is ~10 percent. The noise clearly emerges beyond the core radius (k > 0.05), regardless of largescale structures. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.