© ESO, 2014

1. Introduction

The power spectrum of perturbations and velocity in a given fluid has historically represented one of the crucial tools to understand and constrain the dominant astrophysical processes. In the cosmology field, the temperature fluctuations in the cosmic microwave background have allowed us to put precise constraints on the geometry and composition of the universe (e.g. through the acoustic spectral peaks; Planck Collaboration XVI 2014). Closer to our case, the observed electron density perturbations in the interstellar medium (ISM) have revealed a highly turbulent medium, showing a Kolmogorov power-law spectrum spanning more than 10 decades (Armstrong et al. 1981, 1995, and references therein). The observed power spectrum of the solar wind density, which is also consistent with the famous −5 / 3 slope, has further proven that turbulent processes are a key component shaping the dynamics of astrophysical plasmas (e.g. Woo & Armstrong 1979; Marsch & Tu 1990). In a similar way, the wealth of information contained in the power spectrum extracted from the hot plasma filling galaxy clusters can help us to significantly advance our knowledge of the intracluster medium (ICM) astrophysics.

In the context of galaxy clusters, Gaspari & Churazov (2013; hereafter GC13) have shown for the first time that the power spectrum of the ICM density fluctuations linearly rises with the level of turbulent motions. Diffusive processes, as thermal conduction, instead fight to damp the cascade of perturbations. Many questions still remain to be tackled. In this work, we focus on the statistics and features of the velocity field in the stratified ICM, such as the power spectrum, the real-space and projected maps, and, in particular, its intimate relation with the thermodynamic perturbations. The ICM power spectrum can be viewed in various forms (e.g. Schuecker et al. 2004; Churazov et al. 2012) through the lenses of gas velocities (δv/c_s), density (δρ/ρ), entropy (δK/K), or pressure (δP/P) fluctuations, thus offering multiple joint constraints. Each physical process leaves marked imprints behind. The normalization of the perturbation spectrum is tied to the combined action of gravity and sound waves excited by large-scale (100s kpc) turbulence. At intermediate scales, the thermodynamic perturbations act as “tracers” of the eddy inertial cascade (in line with the classic advection theory by Obukhov 1949 and Corrsin 1951), while rising diffusivity conspire to decouple the tight relation. Faster gas motions alter the thermodynamic mode (from isobaric to adiabatic), changing the interplay between different fluctuations.

The linear relation between the density variance and the turbulent Mach number M has been also observed in simulations of supersonic isothermal turbulence in homogeneous and periodic boxes (e.g. Padoan et al. 1997; Konstandin et al. 2012, and references therein), in connection with ISM studies. However, this relation purely arises from the high compression imparted by supersonic turbulence, creating shocks and sharp peaks (e.g. Kim & Ryu 2005). In the subsonic regime, the compressibility drastically diminishes, and density perturbations fade as M² (Kowal et al. 2007). In the case of stratified galaxy clusters, our novel linear relation is instead tied to the radial gradients of entropy and pressure (Sect. 5), which is already developing in the subsonic regime – the realistic state of ICM turbulence (M ~ 0.2–0.7; e.g. Norman & Bryan 1999; Lau et al. 2009; Vazza et al. 2009).

The velocity statistics of the diffuse medium is notoriously difficult to assess through observations. On the contrary, X-ray surface brightness images can robustly constrain the ICM density. Being able to convert between the spectra (or even just the normalization) of perturbations and velocity, is a powerful tool, which can be exploited by theoretical studies and by the large amount of available Chandra and XMM data. For instance, the quick estimate of the ICM turbulent velocities allows us to study the level of hydrostatic equilibrium in the hot halo (e.g. Vikhlinin et al. 2006), the transport and dilution of metals (e.g. Rebusco et al. 2005), the deposition of energy imparted by the active galactic nucleus (AGN) outflows (e.g. Churazov et al. 2004; Gaspari et al. 2011b, 2012b), the evolution of filaments and bubbles (e.g. Scannapieco & Brüggen 2008), or the reacceleration of cosmic rays (e.g. Brunetti & Lazarian 2007). On the other hand, being able to quickly assess the conductive state of the plasma allows us to constrain the survival of the cold/warm gas, which is crucial for star formation (e.g. McDonald & Veilleux 2009) and black hole accretion (e.g. Gaspari et al. 2013), to study the quenching of cooling flows (e.g. Kim & Fabbiano 2003) and the evolution of cosmic structures (e.g. Dolag et al. 2004).

We could soon take advantage of the inverse process, albeit more expensive. The upcoming Astro-H mission (Takahashi et al. 2010) and the future Athena (Nandra et al. 2013) will provide unprecedented detections of the turbulent velocity dispersion in the ICM (via line broadening) and bulk motions (via line shift; e.g. Inogamov & Sunyaev 2003; Zhuravleva et al. 2012; Nagai et al. 2013; Tamura et al. 2014) down to Mach numbers ~0.1 for massive clusters (see Sect. 4). Reliable constraints on the gas motions allow us to accurately calibrate the above relation and to assess the level of density fluctuations, if the imaging is poor. For instance, apparently “relaxed” systems may host ≳10 percent density fluctuations, which may significantly alter the formation or regeneration of cool cores in clusters and bias the estimate of radial profiles to name a few interesting applications. The same spectral analysis can be extended to the gaseous halos of massive galaxies and groups.

The physics of the ICM is a strongly debated topic. Turbulence has been mainly studied by means of cosmological and isolated simulations (e.g. Norman & Bryan 1999; Dolag et al. 2005; Kim & Ryu 2005; Nagai et al. 2007; Lau et al. 2009; Vazza et al. 2009; Valdarnini 2011; Borgani & Kravtsov 2011; Miniati 2014; Schmidt et al. 2014; Shi & Komatsu 2014). Similarly, diffusion processes as conduction have been mainly investigated via theoretical studies (Chandran & Cowley 1998; Narayan & Medvedev 2001; Zakamska & Narayan 2003; Ruszkowski & Oh 2010, 2011; Voigt & Fabian 2004; Roediger et al. 2013; Smith et al. 2013; ZuHone et al. 2013). Instead, observations have a hard time in resolving and constraining such processes through local features, granting in the last decade only a few estimates (Ettori & Fabian 2000; Markevitch & Vikhlinin 2007; Forman et al. 2007; Eckert et al., in prep.). However, we are now able to retrieve the statistics of density/pressure fluctuations in the ICM (Schuecker et al. 2004; Churazov et al. 2012; Sanders & Fabian 2012), which allows us to probe the gas physics without the need to resolve local structures.

In Churazov et al. (2012), we outlined possible effects that can contribute to density fluctuations: turbulence (via the Bernoulli term ∝ M² or via sound waves), entropy variations (due to mergers or turbulence), perturbations of the gravitational potential, metallicity variations, and AGN bubbles. Using 3D high-resolution plasma simulations, we focus in this work on the role of turbulence and thermal diffusivity and the driven thermodynamic perturbations, including entropy and pressure variations. Controlled experiments allow us to discriminate the exact contribution of each included physics. In a companion paper (Zhuravleva et al. 2014, hereafter Z14), we analyze the density perturbations in cosmological AMR simulations and we focus on the role of gravity waves. At the price of lower resolution, we are thus able to include the turbulence driving led by mergers (Sect. 5.3).

This work is structured as follows. In Sect. 2, we review the main physical and numerical ingredients of the simulated models. In Sect. 3, we present the power spectrum of velocities and density fluctuations by focusing on their tight connection and key features (normalization, cascade, and damping). In Sect. 4, we analyze the real-space properties of the velocity and perturbation field, showing what X-ray observations are able to detect. In Sect. 5, we throughly discuss the physical interpretation of the power spectrum, the interplay of g-waves and p-waves, along with the development of the spectral cascade of all thermodynamic perturbations, in relation to the advection theory of tracers. In Sect. 6, we summarize the results and remark how the ICM power spectrum can be exploited by future observations and theoretical studies to probe the physics of the diffuse medium with high precision.

2. Physics and numerics

The implemented physics and numerics are described in depth in Gaspari & Churazov (2013; Sect. 2), to which the reader is referred for the complete details. Here we summarize the essential features.

The initial conditions for the hot gas are modeled following the latest XMM observed temperature and density radial profiles of Coma cluster (β-model with core radius r_c = 272 kpc and index 0.75). Given the high ICM temperature, T ~ 8.5 keV, and low electron number density, n_e ~ 4 × 10^-3, Coma serves as an excellent laboratory to study the effects of conduction and turbulence without being strongly influenced by radiative cooling or AGN feedback. The hot gas is initialized in hydrostatic equilibrium by providing a gravitational potential appropriate for a massive cluster in the ΛCDM universe with a virial mass M_vir ~ 10¹⁵M_⊙ (r₅₀₀ ~ 1.4 Mpc, which is covered by the width of the 3D box).

The conservative hydrodynamics equations are integrated with the Eulerian code FLASH4 (Fryxell et al. 2000) by using a third order scheme (piecewise parabolic method) in the framework of the unsplit flux formulation. The ICM plasma has an adiabatic index γ = 5/3 and mean atomic weight μ ≃ 0.62. We chose the numerically expensive uniform grid (512³), instead of adaptive cells, to remove any substantial, spurious diffusivity due to the refinement/derefinement. The resolution is Δx≃2.6 kpc, which is roughly on the scale of the (unmagnetized) plasma mean free path, and mimics a slightly suppressed Spitzer viscosity. Boundary zones have a Dirichlet condition, which are fixed by the large-scale radial profile; inflow is prohibited. In addition to hydrodynamics, we add turbulence driving, thermal conduction, and electron-ion equilibration as source terms.

Injection of subsonic, solenoidal turbulence is modeled with a spectral forcing scheme that generates a statistically stationary velocity field (GC13, Sect. 2.2), which is based on an Ornstein-Uhlenbeck random process. The amplitudes of the driven acceleration are evolved in Fourier space and then directly converted to physical space. Since 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; Schmidt et al. 2014; Shi & Komatsu 2014) show that the ICM turbulent energies are ~3–30 percent of the thermal energy, we test subsonic Mach numbers in the range of M ≡ σ_v/c_s ~ 0.25–0.75, where σ_v is the 3D¹ velocity dispersion (average sound speed of Coma is c_s ≃ 1500 km s^-1). The source of turbulence can be various, which includes cosmological flows/mergers, galaxy motions, and feedback processes. The former usually dominate by affecting large volumes beyond the core (e.g. Shi & Komatsu 2014; Z14) and they are typically related to a solenoidal flow (e.g. Miniati 2014). We thus stir the gas on large scales with a typical injection peak L ~ 600 kpc (in a few runs ~300 kpc), and we let turbulence to naturally cascade. The turbulence timescale as a function of physical scale l is the eddy turnover time t_turb ≃ (L^{1 / 3}/σ_v,L)l^{2 /3}, using the Kolmogorov scaling σ_v ∝ l^1/3. Since turbulence is kept subsonic, dissipational heating is subdominant t_diss,heat ~ M^-2t_turb (e.g. Ruszkowski & Oh 2011) on timescales on order of the eddy turnover time. We recall that turbulence acts as an effective diffusivity on entropy with a coefficient D_turb ~ σ_vl (e.g. Dennis & Chandran 2005).

The conduction of thermal energy, due to the plasma electrons, is driven by a flux F_cond = − fκ_S∇T_e with a conductivity ${\mathit{\kappa }}_{\mathrm{S}}\mathrm{\simeq }\mathrm{5}\mathrm{\times }{\mathrm{10}}^{-7}\hspace{0.17em}{\mathit{T}}_{\mathrm{e}}^{\mathrm{5}\mathit{/}\mathrm{2}}$erg s^-1 K^-1 cm^-1 (Spitzer 1962). The related diffusivity and timescale is D_cond = fκ_S/ 1.5 n_ek_B and t_cond = l²/D_cond, respectively. The conductive flux saturates as ${\mathit{F}}_{\mathrm{sat}}\mathrm{\propto }{\mathit{n}}_{\mathrm{e}}{\mathit{T}}_{\mathrm{e}}^{\mathrm{3}\mathit{/}\mathrm{2}}$, whenever the temperature scale height is smaller than the electron mean free path. We use an advanced implicit solver, which allows for long, Gyr integration times. The magnetohydrodynamic (MHD) simulations (e.g. Ruszkowski & Oh 2010) show that the outcome of subsonic turbulence is a tangled magnetic field with small kpc coherence length (see also constraints in Kim et al. 1990). On scales larger than the coherence length, the average suppression due to anisotropic conduction and magnetic microinstabilities can be parametrized with the so-called f factor, commonly known as f ~ 10^-3–10^-1 (GC13, Sect. 2.1.1). Using the effective isotropic conductivity has the advantage of modeling any level of suppression affecting the bulk of the ICM. The MHD runs only provide the geometric suppression above the plasma mean free path λ; in a chaotic atmosphere, it is typically f ~ 1/3 (Ruszkowski & Oh 2010; see also Narayan & Medvedev 2001). However, line wandering, magnetic mirrors, and other plasma microinstabilities, which are well below λ can strongly suppress the transport of heat down to f ~ 10^-3 (Rechester & Rosenbluth 1978; Chandran & Cowley 1998; Komarov et al. 2014). The survival of cold fronts, bubbles, and cold gas (Sect. 1), in conjunction with our GC13 spectral analysis, point toward strongly suppressed values, f ~ 10^-3.

We integrate both the electron and ion temperature equation, since equilibration times can become considerable in a hot plasma (t_ei ≳ 50 Myr). The heat exchange rate is ∝(T_e − T_i) /t_ei, using a Spitzer equilibration time ${\mathit{t}}_{\mathrm{ei}}\mathrm{\propto }{\mathit{T}}_{\mathrm{e}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{/}{\mathit{n}}_{\mathrm{e}}$ (cf. GC13). The 2T modeling allows us to prevent the formation of spurious perturbations due to the unphysical instantaneous transfer of heat.

3. The ICM power spectrum: velocity and δρ/ρ

We now describe the results of the simulated models, which focus on the spectral and real-space properties of the turbulent velocity, in relation with the statistics of gas density perturbations, δρ/ρ. Gas density is indeed the primary astrophysical observable, which is directly extracted from the X-ray surface brightness. As thoroughly discussed in Sect. 5, the perturbations would be actually more evident through entropy (for low M) or pressure (for high M), and then by retrieving δρ/ρ via the main thermodynamic mode (isobaric, isothermal, or adiabatic). Unfortunately, K and P are difficult X-ray observables to constrain. Nevertheless, the spectrum of the leading “tracer” is tied to velocities in a very similar manner, although the underlying cause differs, granting a fairly universal M − δρ/ρ relation (Sects. 5.1.1–5.1.2).

We first retrieve the characteristic amplitude of total velocity normalized to c_s ≃ 1500 km s^-1. It is convenient to use the characteristic amplitude, instead of the power spectrum P(k) or energy spectrum E(k), since its units are the same of the variable in real space. The amplitude spectrum is defined as $\mathit{A}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\equiv }\sqrt{\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathrm{4}\mathit{\pi }{\mathit{k}}^{\mathrm{3}}}\mathrm{\equiv }\sqrt{\mathit{E}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathit{k}}\mathit{,}$(1)where $\mathit{k}\mathrm{=}\sqrt{{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{\equiv }{\mathit{l}}^{-1}$ (kpc^-1). No major bulk motion is present in our box (average velocity ~0); the velocity dispersion is strictly associated with the turbulence driving. The relative perturbations are divided by the underlying background radial profile, for example δρ/ρ = ρ/ρ_b −1 (GC13, Sect. 2.7). Except for mild deviations (Sect. 5.1.3), the turbulence field can be considered isotropic as a first order approximation, allowing us to use the conversion ${\mathit{v}}_{\mathrm{1}\mathrm{D}}\mathrm{~}\mathit{v}\mathit{/}\sqrt{\mathrm{3}}$. All the power spectra are computed with “Mexican Hat” filtering (Arévalo et al. 2012) instead of performing Fourier transforms (GC13, Appendix A for a comparison), which can lead to spurious features due to the box non-periodicity.

In Fig. 1, we show the retrieved characteristic amplitude of v/c_s (blue; we notice that ${\mathit{A}}_{\mathit{v}}\mathrm{=}\sqrt{{\mathit{A}}_{{\mathit{v}}_{\mathrm{x}}}^{\mathrm{2}}\mathrm{+}{\mathit{A}}_{{\mathit{v}}_{\mathrm{y}}}^{\mathrm{2}}\mathrm{+}{\mathit{A}}_{{\mathit{v}}_{\mathrm{z}}}^{\mathrm{2}}}$), which is superposed to the amplitude of density perturbations (red), after reaching a statistical steady state (≳2 t_turb). From the top left panel, the models have increasing turbulence: weak (M ~ 0.25), mild (M ~ 0.5), and strong (M ~ 0.75). The ratio of turbulent to thermal energy is 3.5, 14, and 31 percent (E_turb ≃ 0.56 M²E_th), respectively. The last model (bottom right panel) tests weak turbulence with half the reference injection scale (~300 kpc). The overall behavior of the spectra related to density, velocities, and their ratio (Fig. 2) is fairly self-similar over different Mach numbers and injection scales. We covered the details of the δρ/ρ spectrum in GC13. We focus here on the turbulent velocity and their relationship.

The normalization of the velocity spectrum sets the level of perturbations, which is tied to the A(k) peak. In Fig. 2, we better illustrate the A_v(k) /A_ρ(k) ratios. The hydro f = 0 runs (top) show a converging maximum value A_v/A_ρ ≈ 2.3 (using the 1D velocity A_v1D/A_ρ ≈ 1.3), implying that a stronger turbulence linearly induces larger density fluctuations. This is a key result that allows us to quickly estimate ICM perturbations via the leading turbulent motions and vice versa. For low M, gravity waves mainly produce entropy perturbations, as δK/K ∝ v_1D (see Sect. 5.1.1). For slow motions, the isobaric mode is respected, hence v_1D ∝ γδρ/ρ, as also simulated. For M ≳ 0.5 (see Sect. 5.1.2), compressive sound waves start to significantly contribute: entropy perturbations remain constant, while δP/P increases and sustains the same linear relation, albeit the thermodynamic mode smoothly shifts toward the adiabatic regime. In Sect. 5, we thoroughly discuss all the thermodynamic perturbations (Figs. 8, 9) and the underlying physical interpretation.

For applications in other studies, it is important to note that, proper attention should be paid to the conversion to adopt, given the initial quantity or observable. Shifting from the spectral to physical integrated² quantities, the velocity/density ratio remains the same. However, if we relate the real-space Mach number (i.e. the total variance) to the spectral peak, as in GC13, the conversion to use is M ≈ 4 A(k)_ρ,max (assuming L ~ 600 kpc). Furthermore, turbulence cannot create relative density perturbations with an amplitude higher than roughly the 1D Mach number (Sect. 5). If significantly violated, this would indicate that the perturbation or velocity field has been contaminated by the unfiltered background profile or laminar flows (which are particularly complex in unrelaxed systems). Similarly, strong inhomogeneities must be properly removed: the linear relation applies to relatively small perturbations δρ/ρ< 1 and not to features as cold fronts or buoyant bubbles.

Below the injection scale, A_v and A_ρ continue to be tightly related in the non-diffusive models (Fig. 2). The ratio is independent of M, as indicated by the tight scatter. It is remarkable that the density acts as an effective “tracer” of the velocity field, developing a similar inertial cascade. This is also true for the leading entropy/pressure perturbations. The phenomenon can be explained via the classic theory of advection of passive tracers in turbulent media, A(k)_ρ ∝ A(k)_v (Obukhov 1949; Corrsin 1951; see Sect. 5.2). On the other hand, in the hydro runs we observe that the ratio steadily declines as l^0.13, which reaches about unity near the dissipation scale. The decrease is associated with a shallower cascade of perturbations due to the initial radial gradients, compressive features, and differences in the diffusivity of the “tracer” (Sect. 5.2.1). In Fig. 3 (top), we show a test with 2 × lower resolution (i.e. 2× higher effective viscosity). Besides the good large-scale convergence, the run clarifies that the density/tracer is susceptible to diffusivity in a slightly different way compared with v; hence, we expect departures from the classic tracers theory. We notice also how the cascade of density perturbations is not a perfect power law but tends to exponentially decline, even in the hydro run.

Using half the reference injection scale, the A_v peak is analogous to the that of the reference run (conserving M ~ 0.25), while the density perturbations slightly decrease by ~2^1/3, that is the previous cascade truncated at L/ 2. Stirring smaller scales reduces the influence of gravity waves, since the zone where the turbulence frequency (∝L^{−1 / 3}) is shorter than the buoyancy frequency shrinks (Fig. 7). By varying the injection scale between L′ = L ≡ 600 kpc, L′ = L/ 2, and L′ = L/ 3 we find that correcting the ratio by a factor ~(L′/L)^{1 /3} restores the normalization to a universal value (Fig. 2, top). The development of the shorter cascade is also hindered by the progressive proximity to the dissipation scale.

A key result is the substantial decoupling of velocities and density perturbations, as we increase the level of conduction (dark to bright line color: f = 0 ,10^-3, 10^-2, 10^-1, 1; Figs. 1, 2). In other words, the quick transfer of heat damps density fluctuations (the forming overdensities quickly re-expand due to the temperature increase), while it leaves the turbulent velocity cascade or momentum transfer unaltered. The rising A_v/A_ρ ratio (Fig. 2, top to bottom) is a crucial result that provides a new constraint on the conductive state of the ICM, in addition to the slope of the spectrum. It also breaks any minor degeneracy that strong conduction (f ≳ 0.1) may induce in the spectrum, due to the damping of power at large scales (slightly flattening the cascade). The upcoming Astro-H mission should provide important constraints on the ICM turbulent velocities; combined with high-quality determinations of density perturbations via Chandra/XMM (and Athena), our knowledge of the ICM physics could significantly improve, exploiting the A_v/A_ρ diagnostic.

The velocity cascade follows the Kolmogorov index (A(k) ∝ k^−1/3 or E(k) ∝ k^−5/3) in all runs, except with weak turbulence, where it becomes slightly steeper, though with increased scatter. Considering the substantial stratification, it is remarkable that classic Kolmogorov theory consistently applies to a cluster atmosphere (see also Vazza et al. 2011; Valdarnini 2011). The density spectrum instead displays a steep decay toward A_ρ ∝ k^−1/2, as the turbulent Prandtl number P_t ≡ t_cond/t_turb ≲ 100 (see GC13³). Such a steepening induces A_v/A_ρ to become gradually shallower⁴, inverting the trend for f ≳ 10^-2 (Fig. 2, third panel). This can be explained in terms of the advection theory of tracers (Sect. 5.2), where the diffusivity of the scalar has no effect on the velocity cascade. The scatter of A_v/A_ρ rises with decreasing Prandtl number, that is with stronger conduction and weaker turbulence (P_t ∝ M/f). Conductivity with f ≳ 0.1 can indeed stifle the regeneration of perturbations by a factor of 2–4 over a large range of scales. Figure 2 confirms that any substantial conductivity in the ICM clearly emerges in the A_v/A_ρ diagnostic, showing values up to ≈5 and 3 for weak and strong turbulence, respectively, even at scales of 100s kpc.

4. Real-space properties and X-ray constraints

Before delving into the theoretical interpretation (Sect. 5), it is worth understanding the real-space properties related to turbulence and perturbations of density (or the “tracer”), and what X-ray observations can detect through the spectral line broadening and the projected images. As reference, we consider the models with M ~ 0.5 (E_turb/E_th ~ 0.14, a common cluster regime; e.g. Schuecker et al. 2004; Lau et al. 2009; GC13 – Sect. 4.3).

In Fig. 4, we compare the mid-plane cross-sections of δρ/ρ (left) and magnitude of total velocity (middle). The key result is the progressive smoothing of density fluctuations raising the level of conduction, while the turbulent velocity field remains unaltered (top to bottom panels: f = 0, 10^-3, 10^-2, 10^-1). In the hydro run, the perturbation field shows a complex morphology of filamentary structures, which are produced by the turbulent velocity field and later deformed by Kelvin-Helmholtz and Rayleigh-Taylor instabilities. The rolls and filaments are almost washed out in the presence of strong conduction (f ≳ 0.1), transforming the perturbations from isobaric to isothermal (Sect. 5.1.1). On the other hand, strong conductivity cannot completely wipe out fluctuations (bottom maps). While entropy fluctuations decrease, compressive pressure perturbations still maintain the same level (Fig. 9, bottom). In other words, strong conduction can also promote minor fluctuations, due to the fast transfer of heat and change of compressibility in the medium.

All the velocity maps (middle) are remarkably similar, both statistically and locally, with minute differences only if f ≳ 0.1. As density fluctuations, the turbulent velocities do not show any major difference within or outside the cluster core, signaling a significant level of homogeneity. To be more quantitative, we extracted the spectra only from the cluster center (r<r₅₀₀/ 4), where stratification is less prominent. As shown in Fig. 3 (bottom), both the δρ/ρ and velocity spectra are similar to those computed in the full box. At the largest scale, δρ/ρ experiences a minor decline (~10%) partly because the entropy/pressure profile is shallower in the core and partly due to the limited statistics related to the smallest modes. Concerning isotropy, the real-space maps also do not highlight major deviations. However, transforming the velocity field in spherical coordinates, we retrieve mild anisotropies (Fig. 10). As M< 0.5, the large-scale velocities become slightly more tangential, due to the stronger influence of stratification (Sect. 5.1.3). Instead, stronger turbulence (and conduction) increases the relevance of sound waves, restoring isotropy.

Figure 4 points out that the phases of the perturbation field are not coincident with that of the velocity field, although the two fields are strongly correlated in amplitude. The Pearson coefficient indicates a negligible anticorrelation between δρ/ρ and total velocity (R< −0.2), in all models. In other words, the density filaments are not strictly tied to a local high velocity. This is expected, since the cause of fluctuations is the turbulence driving, while δρ/ρ plays the role of the tracer in a continuously chaotic environment.

We analyzed the volumetric PDF of the logarithmic density fluctuations, ln(1 + δρ/ρ). The thermodynamic fluctuations can be generally described by a log-normal distribution (see example in Fig. 5) with small corrections due to high-order moments (skewness and kurtosis), thus strengthening the role of the power spectrum. The log-normal distribution and weak non-Gaussian contributions are consistent with classic turbulence studies that test solenoidal stirring, albeit in non-stratified and controlled boxes (Federrath et al. 2010, and references therein), and with ICM studies (Kawahara et al. 2007; Zhuravleva et al. 2013). Significant deviations start to arise with highly compressive turbulence, mainly affecting the wings (e.g. Kowal et al. 2007). We defer the study of high-order moments to future work.

What can be inferred from X-ray observations? Besides X-ray imaging (see surface brightness maps in GC13, Fig. 4, or Churazov et al. 2012), X-ray energy spectra can provide crucial constraints on turbulence, in particular by considering the forthcoming Astro-H mission (e.g. Inogamov & Sunyaev 2003; Zhuravleva et al. 2012; Tamura et al. 2014). We have two important tools to exploit: one is the broadening of the spectral line, and the other is the line shift. In Fig. 4 (right), we show the observed line broadening due to turbulent motions along the y-axis view. The Doppler broadening relative to the line rest energy, E₀, is defined⁵ as $\mathrm{\Delta }\mathit{E}\mathit{/}{\mathit{E}}_{\mathrm{0}}\mathrm{\equiv }\sqrt{\mathrm{2}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{1}\mathrm{D}\mathit{,}\mathrm{ew}}\mathit{/}\mathit{c}$. The observed 1D velocity dispersion is computed as ${\mathit{\sigma }}_{\mathrm{1}\mathrm{D}\mathit{,}\mathrm{ew}}^{\mathrm{2}}\mathrm{=}\mathit{E}\mathrm{\left[}{\mathit{v}}_{\mathrm{1}\mathrm{D}}^{\mathrm{2}}\mathrm{\right]}\mathrm{-}{\mathit{E}}^{\mathrm{2}}\mathrm{\left[}{\mathit{v}}_{\mathrm{1}\mathrm{D}}\mathrm{\right]}$, where E [ x ] is the X-ray emission-weighted average along the line of sight, using emissivity n_en_iΛ(T) with X-ray threshold T_x ≳ 0.3 keV (Gaspari et al. 201a). The related full width of the line at half maximum is FWHM$\hspace{0.17em}\mathrm{\equiv }\mathrm{2}\sqrt{\ln \mathrm{2}}\hspace{0.17em}\mathrm{\Delta }\mathit{E}\mathrm{\simeq }\mathrm{1.66}\hspace{0.17em}\mathrm{\Delta }\mathit{E}$. The projected maps appear partly different from the cross-sections, but the statistics, like the velocity dispersion, is the same after deprojection. Astro-H will be able to resolve ~4 eV (and Athena half this value); using the bright Fe XXV line at 6.7 keV, the lower detection limit becomes ΔE/E₀ ~ 0.06 percent (the dark regions in the right-hand panels of Fig. 4). Assuming good statistics, the maps show that Astro-H could detect turbulence in most of the cluster, where σ_1D,ew ≳ 200 km s^-1 or 1D Mach number ≳0.13 (the non-black regions), covering our entire simulated range (see Nagai et al. 2013 for synthetic maps using Astro-H response). The Fe XXV is an excellent line, since the associated thermal broadening is just σ_th = (k_bT/ 56 m_p)^0.5 ~ 120 km s^-1; thereby, the turbulent dispersion typically dominates the contribution to the total broadening of this line.

While the projected velocity dispersion σ_1D,ew highlights the small-scale motions via the line broadening, the projected X-ray emission-weighted velocity field probes the large-scale motions via the line shift (e.g. Zhuravleva et al. 2012). We notice that the driven velocity field still has an average 3D laminar motion ~0. Figure 6 shows the large eddies that are several 100 kpc, which carry most of the specific kinetic energy (${\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{~}{\mathit{A}}_{\mathit{v,}\max }^{\mathrm{2}}$) and would be absent in atmospheres stirred only at small scales. The relative strength of the line broadening and shift thus carries important information about the nature of the driven turbulence, which is a solid proxy corroborated by its insensitivity to conduction.

In passing, we note that X-ray observations are not the single channel to probe ICM fluctuations. We propose using the thermal Sunyaev-Zel’dovich (SZ) effect to independently extract pressure fluctuations. Current X-ray maps still have ≳4 × higher spatial resolution. However, future observations (e.g. ALMA, CCAT, SKA) will allow us to constrain the ICM power spectrum even in the submillimeter/radio band, avoiding the use of expensive X-ray spectroscopy.

5. Discussion and physical interpretation

In this section, we discuss the physical interpretation of the ICM power spectrum; in particular, we are concerned with the tight relation between the velocity and the other primary thermodynamic quantities (entropy, pressure, and density). It shall be kept in mind that the reason why we perform 3D simulations is that it is impossible to analytically solve a chaotic, nonlinear system. The following arguments arise from first order perturbations or dimensional theories and shall be regarded as simple estimates. The simulations generally confirm the ansatz presented below, albeit with relevant differences, which are critically discussed. We first focus on the normalization of the spectra (l ~ L; Sect. 5.1), and then we analyze the spectral cascade (l<L; Sect. 5.2), along with the alterations imparted by conduction from the ideal evolution.

5.1. Spectra normalization

The normalization of the spectra is likely related to the relative importance of gravity waves and sound waves. Linearizing the perturbed hydrodynamic equations, it is possible to describe the propagation of a general wave in a spherical and gravitationally stratified atmosphere in terms of the following dispersion relation (see Balbus & Soker 1990, for the WKBJ analysis): ${\mathit{\omega }}^{\mathrm{4}}\mathrm{-}{\mathit{\omega }}^{\mathrm{2}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{\omega }}_{\mathrm{BV}}^{\mathrm{2}}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}\mathit{,}$(2)where ${\mathit{k}}^{\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}$ (k_r and k_⊥ are the radial and azimuthal components of the wavenumber vector k, respectively), c_s = (γk_bT/μm_p)^{1 /2} is the adiabatic sound speed, and ω_BV is the Brunt-Väisälä (buoyancy) frequency, defined as ${\mathit{\omega }}_{\mathrm{BV}}\mathrm{\equiv }{\left[\frac{\mathit{g}}{\mathit{\gamma }}\hspace{0.17em}\frac{\mathrm{d}\ln \mathit{K}}{\mathrm{d}\mathit{r}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(3)where g is the gravitational acceleration. Equation (2) includes the action of two key waves. The middle term is associated with pressure waves, or simply sound waves (p-waves), while the last term represents gravity waves (g-waves) driven by the restoring buoyant force. In the next sections, we show that small perturbations driven by the two waves are tied to the Mach number, ~δv/c_s. For both g- and p-waves, this holds within order unity (as in the simulated nonlinear regime), but the leading perturbations and dynamical modes differ: g-waves mainly drive entropy perturbations (δK/K), increase the gas vorticity, and induce a tangential bias in the turbulent velocity field; p-waves are instead associated with compressive pressure fluctuations (δP/P), a preferentially irrotational field, and isotropic turbulence (or with slightly radial bias).

In Fig. 7, we show the frequency of the simulated turbulent motions, ${\mathit{t}}_{\mathrm{turb}}^{-1}$, compared with the Brunt-Väisälä frequency (black; Eq. (3)) in the full radial range (at variance with GC13, where we focused on the properties of the central region). At large scales, the two frequencies tend to be roughly comparable (Froude number ~1); hence, both g-waves and p-waves can be excited. This is a typical condition for most clusters, since large-scale profiles (entropy and pressure) are fairly self-similar, and turbulence follows our simulated subsonic range. Let us first analyze the two limiting regimes to understand better both processes.

5.1.1. Low frequencies (low M): g-waves

The stratification of the ICM atmosphere allows the excitation of gravity waves (e.g. Lufkin et al. 1995; Ruszkowski & Oh 2010, 2011). In the low frequency regime, the dispersion relation in Eq. (2) can be written as ${\mathit{\omega }}^{\mathrm{2}}\mathrm{\simeq }{\mathit{\omega }}_{\mathrm{BV}}^{\mathrm{2}}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{\perp }}^{\mathrm{2}}}{{\mathit{k}}^{\mathrm{2}}}\mathit{,}$(4)which tells us that gravity waves are evanescent for ω>ω_BV, since k_r must be imaginary. Therefore, wherever ω<ω_BV, g-waves are excited. For a cluster atmosphere, ω_BV declines at large r and waves are trapped within the radius, such that ω ≃ ω_BV (Balbus & Soker 1990) as k ≃ k_⊥. More importantly, buoyant oscillations damp the radial component of turbulence, inducing a tangential bias in the gas velocity field (Froude < 1). In the limiting case, the chaotic motions should collapse in azimuthal shells (e.g. Ruszkowski & Oh 2010). The profile of the velocity anisotropy parameter, $\mathit{\beta }\mathrm{\equiv }\mathrm{1}\mathrm{-}{\mathit{\sigma }}_{{\mathit{v}}_{\mathrm{\perp }}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{r}}}^{\mathrm{2}}$ would show β ≪ 0. Tangentially-biased vorticity is thus a good marker of the g-waves influence (see Sect. 5.1.3).

Gravity waves are mainly tied to entropy perturbations (see also Z14). Using ω_BV as dominant frequency, the buoyant acceleration over a displacement δr can be described with a simple harmonic oscillator: $\begin{array}{c}\mathrm{¨}\\ {\mathit{r}}_{\mathrm{b}}\end{array}\mathrm{=}\mathrm{-}\hspace{0.17em}{\mathit{\omega }}_{\mathrm{BV}}^{\mathrm{2}}\mathit{\delta r}\mathrm{=}\mathrm{-}\hspace{0.17em}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{\gamma }}^{\mathrm{2}}\hspace{0.17em}\mathrm{\left|}{\mathit{H}}_{\mathit{P}}\mathrm{\right|}\hspace{0.17em}{\mathit{H}}_{\mathit{K}}}\mathit{\delta r,}$(5)where ω_BV is written⁶ as a function of both the scale height of entropy, H_K ≡ dr/ dlnK, and pressure, H_P ≡ dr/ dlnP. Physically, g-waves occur because an entropy element is displaced from its equilibrium position, r₀, thus inducing an opposite force that acts to restore the blob back to where the radial entropy is the same (clusters are convectively stable, ∇K> 0). The small displacement is thus linked to δK/K ≃ (dlnK/ dr)₀δr, i.e. δr ≃ (δK/K) H_K0. The specific potential energy of the harmonic oscillator is ${\mathit{E}}_{\mathrm{b}}\mathrm{=}{\mathit{\omega }}_{\mathrm{BV}}^{\mathrm{2}}\mathrm{(}\mathit{\delta r}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{2}$. Substituting ω_BV and δr, we can write ${\mathit{E}}_{\mathrm{b}}\mathrm{\simeq }\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\gamma }}^{\mathrm{2}}}\frac{{\mathit{H}}_{{\mathit{K}}_{\mathrm{0}}}}{\mathrm{\left|}{\mathit{H}}_{{\mathit{P}}_{\mathrm{0}}}\mathrm{\right|}}{\left(\frac{\mathit{\delta K}}{\mathit{K}}\right)}^{\mathrm{2}}\mathrm{·}$(6)Using E_b as an estimate for the average specific kinetic energy ($\mathit{v̅}\begin{array}{c}\mathrm{2}\\ \mathrm{1}\mathrm{D}\end{array}\mathit{/}\mathrm{2}\mathrm{~}{\mathit{E}}_{\mathrm{b}}$) finally yields $\mathit{M̅}\mathrm{1}\mathrm{D}\mathrm{~}\frac{\mathrm{1}}{\mathit{\gamma }}\sqrt{\frac{{\mathit{H}}_{{\mathit{K}}_{\mathrm{0}}}}{\mathrm{\left|}{\mathit{H}}_{{\mathit{P}}_{\mathrm{0}}}\mathrm{\right|}}}\left|\frac{\mathit{\delta K}}{\mathit{K}}\right|\mathrm{\simeq }\mathrm{𝒪}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\hspace{0.17em}\left|\frac{\mathit{\delta K}}{\mathit{K}}\right|\mathrm{·}$(7)The last step arises because the ratio of the scale heights for an isothermal atmosphere is constant, H_K/ | H_P | = 1 / (γ −1) = 1.5. In general, H_K/ | H_P | = | α_P | /α_K ≃ 1–2 (within r₅₀₀), where α_P and α_K are the slopes of the logarithmic radial profiles of pressure and entropy, respectively; α_P steepens with increasing radius (Arnaud et al. 2010), while α_K ≃ 1.

Figure 8 shows the power spectra of all the thermodynamic quantities, including turbulent velocities, for the hydro and a conductive run in the low Mach regime. Considering the hydro model (top), the ratio of ${\mathit{M}}_{\mathrm{1}\mathrm{D}}\mathrm{\simeq }{\mathit{M}}_{\mathrm{3}\mathrm{D}}\mathit{/}\sqrt{\mathrm{3}}$ and δK/K near the injection scale is 0.85, in line with the estimate in Eq. (7). Evidently, only one component of velocity is acting as an efficient mixer, since cluster gradients are functions of r. Figure 8 clarifies that the low frequency regime corresponds to $\left|\frac{\mathit{\delta P}}{\mathit{P}}\right|\mathrm{\ll }\left|\frac{\mathit{\delta K}}{\mathit{K}}\right|\mathit{,}$(8)since slow motions tend to be in pressure equilibrium with the surroundings, which is helped by the convective stability of the ICM. The isobaric behavior also manifests in the relation between density and temperature, or entropy and density (both anticorrelated): $\left|\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right|\mathrm{\approx }\left|\frac{\mathit{\delta T}}{\mathit{T}}\right|\mathrm{and}\left|\frac{\mathit{\delta K}}{\mathit{K}}\right|\mathrm{\approx }\mathit{\gamma }\left|\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right|\mathrm{\left[}\mathrm{isobaric}\mathrm{\right]}\mathit{.}$(9)Both relations are followed within ≲10 percent, as shown in Fig. 8 (consistently with the Pearson analysis in GC13). Recall that if we halve the injection scale (Sect. 3), this slightly reduces the strength of g-waves as L^{1 /3}, since ω_turb ∝ L^−1/3 approaches ω_BV (Fig. 7). The ratio A_v/A_ρ thus increases by the same factor (sound waves are still too weak to contribute).

In the presence of mild conduction, the normalization is unaltered for f< 10^-2, although the intermediate cascade is damped by the increased diffusivity of the “tracer” (Sect. 5.2). For f ≳ 10^-2, the large-scale perturbations are also progressively driven toward the isothermal regime. Gravity waves can still induce a significant entropy contrast, but buoyancy is weakened by the increased diffusivity, tracing the shallower temperature gradient instead of ∇K (Ruszkowski & Oh 2010). The normalization of perturbations can thus decrease by a factor ~2. The relation between the different perturbations changes as $\left|\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right|\mathrm{\gg }\left|\frac{\mathit{\delta T}}{\mathit{T}}\right|\mathrm{and}\left|\frac{\mathit{\delta K}}{\mathit{K}}\right|\mathrm{\approx }\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\left|\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right|\mathrm{\left[}\mathrm{isothermal}\mathrm{\right]}\mathit{.}$(10)Since turbulent regeneration is continuous and the ion-electron equilibration time is non-negligible, the pure isothermal regime is impossible to achieve. Nevertheless, as shown in Fig. 8 (bottom panel), the gap between entropy and density perturbations starts to shrink, as T_e fluctuations gradually lower (see also Fig. 9, bottom).

5.1.2. High frequency (high M): p-waves

In the opposite regime, at high frequency (M> 0.5), the dispersion relation (Eq. (2)) is shaped by the contribution of sound waves (p-waves), which can be now written as ${\mathit{\omega }}^{\mathrm{2}}\mathrm{\simeq }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\hspace{0.17em}{\mathit{k}}^{\mathrm{2}}\mathit{.}$(11)The azimuthal component of the wavenumber scales as k_⊥ ∝ r^-1. Therefore, p-waves excited in the cluster central regions become mainly radial further out, k ~ k_r. Moreover, disturbances to the vorticity results to be proportional to (Lufkin et al. 1995) $\frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{\delta }\mathrm{(}{\nabla }\mathrm{\times }{{v}}^{\mathrm{)}}\mathrm{\propto }{k}\mathrm{\times }{\nabla }\ln \mathit{K,}$(12)implying that p-waves are preferentially⁷ characterized by an irrotational velocity field, in contrast with the more tangential g-waves. Overall, stronger p-waves tend to restore isotropic turbulence, or to induce a slightly radial bias (depending on how many sound waves are excited in the central regions; Sect. 5.1.3).

The characteristic injection frequency of p-waves is ${\mathit{\omega }}_{\mathit{s}}^{\mathrm{2}}\mathrm{\simeq }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}{\mathit{L}}^{\mathrm{2}}$. In this high frequency regime $\left|\frac{\mathit{\delta P}}{\mathit{P}}\right|\mathit{>}\left|\frac{\mathit{\delta K}}{\mathit{K}}\right|\mathit{,}$(13)meaning that pressure perturbations drive the dynamics and fluctuations that follow the adiabatic regime (constant entropy). Following the same arguments provided in the previous section, the displacement magnitude is now tied to pressure variations as | (δP/P) H_P0 |. Using the potential energy as an estimate for the average kinetic energy (cf. Eqs. (6)–(7)) now yields $\mathit{M̅}\mathrm{1}\mathrm{D}\mathrm{~}\frac{\mathrm{\left|}{\mathit{H}}_{{\mathit{P}}_{\mathrm{0}}}\mathrm{\right|}}{\mathit{L}}\left|\frac{\mathit{\delta P}}{\mathit{P}}\right|\mathrm{\simeq }\mathrm{𝒪}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\hspace{0.17em}\left|\frac{\mathit{\delta P}}{\mathit{P}}\right|\mathit{,}$(14)where | H_P0 | ≃ L/ | α_P |, since we are analyzing the large-scale power. We note that stratification allows a more efficient generation of sound waves compared with uniform media (e.g. Stein 1967), as shown by the simulations, due to the partial conversion of solenoidal turbulence in more compressive motions.

Comparing the last estimate with Eq. (7), it is clear that the (1D) Mach number drives the spectrum normalization of perturbations in both cases, as found in the simulations. The transition must be smooth, as hinted by the general dispersion relation (Eq. (2)). However, the driving perturbations change character. Figure 9 shows that δP/P rises linearly with M for M ≳ 0.5, while δK/K remains constant (A_K,max ~ 0.15 in both hydro runs). In the M ~ 0.75 flow (middle panel), δP/P has reached δK/K; hence, p-waves do not yet fully overcome g-waves, which would happen for M ≳ 1 (as ω_turb>ω_BV at each radius, or Froude > 1; Fig. 7). Moreover, while gravity waves tend to accumulate in the system, p-waves may leave it in a few sound-crossing times. Again, turbulence with high Mach number is required to see a system fully dominated by p-waves.

According to Eq. (13), the perturbations of the other thermodynamic quantities start to shift from the isobaric to adiabatic regime, implying the following conversion: $\left|\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right|\mathrm{\approx }\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\left|\frac{\mathit{\delta T}}{\mathit{T}}\right|\mathrm{and}\left|\frac{\mathit{\delta P}}{\mathit{P}}\right|\mathrm{\approx }\mathit{\gamma }\left|\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right|\mathrm{\left[}\mathrm{adiabatic}\mathrm{\right]}\mathit{,}$(15)as signaled by the increasing gap between density and temperature fluctuations (up to ≃1.5). Adding conduction (Fig. 9, bottom) shifts the adiabatic mode again toward a partial isothermal regime (Eq. (10)). Interestingly, pressure fluctuations are now more clearly the driver of density fluctuations (δP/P ~ δρ/ρ), since δK/K has degraded by over 50% compared with the hydro run.

In the low M regime, we note that the presence of weak δP/P can be mainly attributed to the conservation of the Bernoulli parameter (although exactly valid only in the steady state). For a constant potential, it can be written as $\frac{{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathrm{=}\mathrm{const}\mathit{.}$(16)Differentiating and taking v₀ = 0 as the reference velocity yields δP/P ∝ M², which is indeed followed by the runs with M< 0.5 (compare Figs. 8 and 9, top panels), while models with stronger turbulence follow δP/P ∝ M (Fig. 9, top and middle panels).

5.1.3. Real systems: g- and p-waves interplay

The estimates in Eqs. (7) and (14) are crude approximations to reality. In the realistic cluster evolution, the 3D hydrodynamic equations and the related perturbations are nonlinear, with the addition of chaotic stirring (Sect. 5.2.1 for other deviations). Multiple frequencies act at the same time, describing waves at different scales and radii. Although the simulations confirm that the spectra normalization is provided by M_1D, we expect a combination of g- and p-waves that shape the dynamics of a turbulent cluster. This is visually highlighted by the maps in Fig. 4. More quantitatively, we can discriminate the action of both channels through an important marker, which is the anisotropy of velocities along the radial direction (Sects. 5.1.1–5.1.2). We notice that computing the velocity spectra over the scale l obfuscates the radial anisotropy.

In Fig. 10, we show the $\mathit{\beta }\mathrm{\equiv }\mathrm{1}\mathrm{-}{\mathit{\sigma }}_{{\mathit{v}}_{\mathrm{\perp }}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\sigma }}_{{\mathit{v}}_{\mathit{r}}}^{\mathrm{2}}$ parameter for all the hydro models. For the low M ~ 0.25 flow, above ~100 kpc, the motions are mildly⁸ tangential, reaching a minimum value β ~ −1. Analyzing the frequencies in Fig. 7, this corresponds to the region where ω<ω_BV, that is where g-waves tend to be more dominant and damp the radial component of turbulent motions. Raising M increases the turbulence frequency up to ~ω_BV, even near the injection scale at several 100s kpc. Therefore, the large-scale motions progressively lose the tangential bias, drifting toward isotropy (β ~ 0). In the M ~ 0.75 flow (red), p-waves start to dominate: motions show no sign of the tangential bias with instead a slightly radial bias (Sect. 5.1.2). Within r< 100 kpc, the frequencies of the turbulent cascade are always greater than ω_BV; hence, p-waves start to have a major influence. The transition occurs at smaller radii for lower M, as suggested by Fig. 10. On the other hand, large-scale g-waves tend to be trapped within the cluster core; the combination of the two effects results in a quasi isotropic β in the core (or slightly radial β ~ 0.2). In the presence of conduction, the same scenario applies, but the β factor is overall reduced toward the isotropic value (by about half the previous values). As noted in Sect. 5.1.1, thermal conduction indeed inhibits buoyancy.

Overall, we suggest using the β parameter to discriminate the effects of g- and p-waves, in conjunction with the thermodynamic mode. We remark that typical cluster conditions are expected to show at best mild anisotropic motions, which emerge only in spherical coordinates; most clusters are not in a regime in which ω ≪ ω_BV (or the opposite), where radial motions are dramatically suppressed. The recurrent stirring also promotes β values drifting toward isotropy. Nevertheless, the ICM spectrum normalization is always comparable to M_1D within order unity, regardless of which one is the driving wave, for both the linear approximation (Eq. (2)) and the nonlinear simulations.

5.2. Spectral cascade: advection of tracers

We have analyzed the physical interpretation of the spectra normalization (l ~ L) so far. The next question is why perturbations show an inertial cascade similar to that of velocities (l<L). The main thermodynamic variables can be crudely considered as “tracers” of the velocity field, whose spectra are explained with the advection theory of passive scalars (e.g. Sreenivasan 1991; Monin & Yaglom 1975; Warhaft 2000, for a review). This is particularly relevant in the subsonic regime, since the compressive term ∇·v → 0.

The equation governing the advection of a passive incompressible scalar C is given by (Warhaft 2000, Sect. 1) $\frac{\mathit{DC}}{\mathit{Dt}}\mathrm{\equiv }\frac{\mathit{\partial C}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}{\nabla }\mathit{C}\mathrm{=}{\nabla }\mathrm{·}\mathrm{(}\mathit{\kappa }\hspace{0.17em}{\nabla }\mathit{C}\mathrm{)}\mathit{,}$(17)where κ is the diffusivity of the tracer and D/Dt is the Lagrangian derivative. According to classic Kolmogorov-Obukhov-Corrsin theory (KOC; Obukhov 1949; Corrsin 1951; Warhaft 2000), the energy spectrum of the scalar linearly traces that of velocities, such as ${\mathit{E}}_{\mathit{C}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\propto }{\mathit{E}}_{\mathit{v}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\propto }{\mathit{k}}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}}\mathit{.}$(18)Considering a 2D vortical motion, the vorticity is conserved (see also Kelvin’s theorem), D(∇ × v) /Dt = 0, in analogy to Eq. (17). The turbulent eddies and the scalar are thus expected to share similar properties, like the spectral cascade, although the diffusive term introduces some discrepancies (Sect. 5.2.1).

Entropy is an excellent example of “passive tracer”. The quantity S ≡ k_B/ [ (γ − 1)μm_p ] lnK has indeed the advantage of being insensitive to adiabatic compressions or expansions. Aside from diffusion, the Lagrangian derivative of entropy is thus conserved, if no irreversible heating (ℋ) or cooling (ℒ) occurs, such as $\mathit{\rho T}\frac{\mathrm{D}\mathit{S}}{\mathrm{D}\mathit{t}}\mathrm{\equiv }\mathit{\rho T}\left(\frac{\mathit{\partial S}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}{\nabla }\mathit{S}\right)\mathrm{=}\mathrm{\mathscr{H}}\mathrm{-}\mathrm{ℒ}\mathrm{\simeq }{\nabla }\mathrm{·}\mathrm{(}{\mathit{D}}_{\mathrm{turb}}\hspace{0.17em}\mathit{\rho T}\hspace{0.17em}{\nabla }\mathit{S}\mathrm{)}\mathit{.}$(19)The hot ICM has negligible radiative cooling (ℒ ≃ 0); the only source of heating can be turbulent diffusion or dissipation. The latter is subdominant for subsonic flows, t_diss,heat ≃ M^-2t_turb. The entropy S can thus replace the scalar C in Eq. (17) with a diffusivity tied to the turbulent field, D_turb ~ σ_vl (Sect. 2). In the low M regime (Sect. 5.1.1), entropy fluctuations tend to lead the dynamics of perturbations, which are linked to the large-scale g-waves. The δK/K cascade then develops over l<L, tracing the velocity inertial regime (Fig. 8) in line with KOC theory (Eq. (18)). Since the fluctuations of density and temperature follow from the dominant mode – isobaric (Eq. (9)), isothermal (Eq. (10)), or adiabatic (Eq. (15)) – also their cascade traces that of δK/K. As discussed before, increasing the Mach number boosts the impact of p-waves. The leading tracer gradually shifts toward δP/P fluctuations. In the M ~ 0.75 hydro run (Fig. 9), density perturbations start to track more closely the δP/P cascade. Pressure is affected by adiabatic processes, yet the simulations show that its cascade is slightly steeper than that of velocities, signaling that it can be used as a crude tracer of the (subsonic) flow. Interestingly, pressure is not directly affected by turbulent mixing (which acts on entropy) and, thus, displays a tighter cascade with velocities⁹.

Adding conduction, increases the effective diffusivity of the tracer (κ → fκ_S), leading to a cascade steeper than that of Kolmogorov (Fig. 8, bottom; notice how T_e fluctuations are strongly damped). The decline of the tracer spectrum starts to occur as D_cond>D_turb/ 100, i.e. P_t< 100. Physically, the quick increase in temperature induces a rapid re-expansion of the forming overdensity, which overcomes the action of turbulent regeneration. In the conductive regime, the A_v/A_ρ ratio is thus expected to gradually increase as a function of f (Fig. 2). The ratio widens up to a factor of 3–5 over all scales in the presence of strong conduction (f ≳ 0.1).

It is interesting to note that in the simple KOC picture: $\mathit{E}\mathrm{(}\mathit{k}{\mathrm{)}}_{\mathit{C}}\mathrm{=}\mathit{b}\mathit{ϵ̅}\begin{array}{c}\hspace{0.17em}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}\\ \mathit{C}\end{array}\hspace{0.17em}\mathit{ϵ̅}\mathit{v}\hspace{0.17em}\hspace{0.17em}{\mathit{k}}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}}\mathit{,}$(20)where b is a universal constant. The average dissipation rate of the passive tracer and velocity is $\mathit{ϵ̅}{}_{\mathit{C}}\mathrm{~}{\mathit{\sigma }}_{\mathit{C}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathit{L}\mathit{/}{\mathit{\sigma }}_{\mathit{v}}\mathrm{)}$ and $\mathit{ϵ̅}{}_{\mathit{v}}\mathrm{~}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{3}}\mathit{/}\mathit{L}$, respectively, implying that the normalization of Eq. (20) is independent of σ_v. In reality, the compressive term sustains a relation between the density variance and the Mach number, even in homogeneous media: in the subsonic range, δρ/ρ fade as M² (cf. Sect. 5.1.2), while the relation becomes linear for supersonic turbulence (see Kowal et al. 2007 and Sect. 1). The latter is, however, conceptually different from our retrieved linear relation developing in the subsonic state of galaxy clusters. In the stratified ICM plasma, increasing the Mach number implies larger coherent displacement, leading to a larger contrast of entropy/pressure defined by the cluster gradients or scale heights (Eqs. (7) and (14)).

5.2.1. KOC departures and the radial gradients

Although KOC theory can explain the overall picture of the spectral cascade, its arguments are purely based on dimensional analysis. Our retrieved slope of entropy/density is typically shallower than the Kolmogorov index (in the non-diffusive models). Physical experiments (e.g. Fig. 5 in Sreenivasan 1991) show that the tracer slope approaches the Kolmogorov cascade only for very high Reynolds numbers and in a slow asymptotic way. For low Reynolds numbers, as in our ICM simulations (R_L ≲ 500; GC13, Sect. 2.6), the spectral index is expected to be shallow (cf. Fig. 4 in Warhaft 2000). Even pressure fluctuations, albeit in line with the Kolmogorov cascade, are shallower than the classic expectation E_P ∝ k^−7/3 (A_P ∝ k^{−2 /3}; Schuecker et al. 2004), in the M ≳ 0.5 runs. A slope as shallow as A_k ∝ k^{−1 /5} signals that the timescale for transferring the tracer variance from large to small scales is ∝k^{−4 / 5}, instead of the Kolmogorov ∝k^{−2 / 3}, due to diffusion effects affecting the transfer process in different ways (see the viscosity test in Fig. 3). This is remarked by the uncorrelated phases between velocity and the tracer (Fig. 4).

Another departure from the KOC cascade may be associated to compressive features. In the extreme case of highly supersonic turbulence, shocks induce very thin peaks in gas density (Kim & Ryu 2005, Fig. 2). Sharp peaks can be seen as delta functions, which generate a flat spectrum, P_δ ~ k⁰, in Fourier space. The ICM turbulence is, however, subsonic; thus, the contribution of thin compressive features to our observed flattening is limited.

The incomplete similarity with the Kolmogorov cascade more likely depends on the initial entropy/pressure gradients. Eqs. (7) and (14) are only valid for small displacements. For a nonlinear evolution, the injected eddy experiences H_K and H_P which varies with radius, given that the initial cluster profiles are self-similar power laws (Sect. 5.1.1): $\frac{\mathit{\delta K}}{\mathit{K}}\mathrm{=}{\mathit{\alpha }}_{\mathit{K}}\hspace{0.17em}\frac{\mathit{\delta r}}{\mathit{r}}\mathrm{and}\frac{\mathit{\delta P}}{\mathit{P}}\mathrm{=}{\mathit{\alpha }}_{\mathit{P}}\hspace{0.17em}\frac{\mathit{\delta r}}{\mathit{r}}\mathrm{·}$(21)For instance, fixing δr ~ L, the injected turbulence at smaller radii can create relatively larger contrasts, inducing a flattening in the spectral cascade. The magnitude of this effect depends, however, on the relation r ← → l. Since chaotic motions are 3D, the turbulent streamlines intersect different projections of the radial gradients, hence the dependence shall be weak. Neglecting the previous KOC departures, the flattening of the simulated spectra corresponds to an average r ~ l^0.13. An opposing effect is related to the fact that the indices α_K and α_P are not constant, but both decline within the core radius (turbulent mixing also slightly lowers α_K in time; GC13). According to Eq. (21), a lower slope implies a lower contrast. This may explain why the A_v/A_ρ ratio slowly declines toward smaller scales (Fig. 2, top), though always remaining larger than unity. We will investigate other cluster atmospheres in the future to assess the impact of different K and P scale heights. In closing, we note that all these secondary effects are washed out in the presence of any significant diffusivity, which completely inverts the A_v/A_ρ downtrend (Fig. 2).

5.3. Further improvements

Finally, we discuss the limitations of the models and further improvements. We studied the evolution of the ICM, primarily in the hot regime. In future works, we plan to extend the simulated sample (as strong cool-core systems) and to test additional physics. Needless to say, 3D high-resolution 2T simulations with turbulence and diffusive terms are extremely expensive, hence small steps must be taken.

It will be interesting to include the effect of cooling, which can induce thermal instability and condensation of cold filaments (Gaspari et al. 2012a, 2013). The AGN feedback balances cooling, preserving the cluster core in global quasi-thermal equilibrium (e.g. Gaspari et al. 2012b). However, both processes just affect the inner region r< 0.1 R₅₀₀ (cf. Gaspari et al. 2014), while galaxy clusters maintain self-similarity at large radii, especially in the entropy profile (α_K ~ 1; e.g. Panagoulia et al. 2014). We thus do not expect dramatic deviations from the current ICM power spectrum (which is intrinsically volume-weighted), and we believe our results can be applied to a wide range of clusters and conditions. Strongly unrelaxed systems, such as as major mergers, might present significant variations, due to the dynamic gravitational potential for example, and requires to be further tested. We are also studying the role of very small injection scales (e.g. AGN outflows): for L< 50 kpc, entropy perturbations may be considerably weaker (as ω_BV<ω_turb), while pressure perturbations should drive the normalization of the density power spectrum, even at low Mach numbers.

In the companion work (Z14), we improve the driving, including the cosmological evolution and the turbulence generated by mergers and large-scale inflows. Albeit limited by low resolution, we find that the linear M_1D–δρ/ρ relation holds across a large sample of simulated clusters. We retrieve a relation scatter of ~30 percent. In the cosmological context, it is more difficult to disentangle the source of the velocity anisotropy, especially in unrelaxed systems. As for observational data, it is important to accurately remove the underlying radial profile and the strongly nonlinear sub-structures, which can contaminate the large-scale power. We find that the most reliable scales dominated by the turbulent cascade are l< 300 kpc, which is fortunately the optimal regime for X-ray observations (see GC13, Sect. 4.3 for a comparison with real data).

We plan to test additional physics. We currently probed the effects of (nearly maximal) Spitzer-like viscosity and electron-ion equilibration. It will be interesting to assess the role of the related magnetic suppression factors (which can be different from that of heat transport). On the other hand, diffusivities linked to the ions are roughly two orders of magnitude slower compared with electron thermal conduction, since the electron sound speed is ≃43 times that of ions. We thus expect conduction to dominate the shape of the power spectrum over a large range of scales. A viscosity lower than the present Spitzer-like value (which would require a much higher resolution) would imply a more extended inertial cascade. Compared with Fig. 2, A_v/A_ρ should thus differ only below 10 s kpc, continuing to widen in combination with high conductivity. In the presence of both low viscosity and conductivity, density and velocity spectra are expected to be tightly coupled again (KOC theory).

Fully MHD simulations are a further route of improvement, which can better model the local features (as cold fronts and filaments). However, besides the numerical complication of integrating anisotropic conduction for long times, the MHD runs can at best retrieve the geometric suppression factor (see Sect. 2). Therefore, we would still be forced to parametrize the conductivity with a factor f_∥ to include microinstabilities and line divergence below the gas mean free path.

6. Conclusions

We carried out 3D high-resolution hydrodynamic 2T simulations to study the power spectrum of the hot ICM in its various manifestations. We focused on the properties of the velocity field and the intimate relation with the driven thermodynamic fluctuations (in particular of density, the primary observable). The ICM power spectrum contains enough information to accurately constrain the dominant physics of the diffuse medium, as the strength of turbulent motions, the level of thermal diffusivity, and the thermodynamic mode, among the most notable. The spectra of v/c_s and of perturbations (e.g. δρ/ρ) are overall self-similar, varying the strength of turbulence via the 3D Mach number, M, or changing the injection scale, L.

At the large cluster scales (l ~ L), several 100 kpc, we retrieve the following:

• Weak turbulent motions in the cluster (M ≲ 0.25) mainly excite gravity waves (ω_turb<ω_BV); the leading perturbations are related to entropy variations δK/K. For stronger turbulence (M> 0.5), sound waves start to significantly contribute (ω ≳ ω_BV), passing the leading role to the compressive pressure fluctuations δP/P.

• The other thermodynamic perturbations, such as δρ/ρ and δT/T, derive from the dominant mode of the process: isobaric (for g-waves/low M), adiabatic (for p −waves/high M), or a mixed state for intermediate M. Conduction shifts the perturbations toward the isothermal mode.

• In both the regimes driven by g- or p-waves, the turbulent 1D Mach number is comparable to the variance of the leading perturbations (K or P) within order unity. For example, M_1D ~ δK/K ~ γδρ/ρ for M ≲ 0.25 flows. Quantitatively, all simulations show a linear relation given by A_v,max≃2.3 A_ρ,max (at L ~ 600 kpc) with a weak L^{1 /3} scaling. To convert between Fourier and real space, the relation to apply is instead M ≈ 4 A_ρ,max (at L~600 kpc).

• Turbulent motions with a tangential bias (β(r) < 0) mark the influence of g-waves (low M), while p-waves (high M) tend to preserve isotropy or to induce a slightly radial bias. Most clusters show intermediate Mach numbers, hence we expect a mixed regime drifting toward global isotropy (Froude ~1).

At the intermediate/small scales (l<L), 10–100 kpc, we find the following:

• The turbulent velocities develop a Kolmogorov cascade(A_v ∝ k^−1/3 or E_v ∝ k^{−5 /3}) in all subsonic runs, despite stratification (ω_turb<ω_BV). The thermodynamic perturbations, particularly entropy, act as effective “tracers” of the velocity field, developing an analogous inertial cascade, in line with the classic (Kolmogorov-Obukhov-Corrsin) advection theory of passive scalars in turbulent media.

• The cluster radial gradients, in conjunction with compressive features, conspire to moderately flatten the perturbations spectrum, slightly departing from the KOC theory and inducing a slow decrease in A_v/A_ρ.

• Thermal conduction strongly damps density/entropy perturbations (the spectral steepening occurs where Prandtl P_t< 100) but leaves the velocity cascade unaltered. This has a dramatic consequence on A_v/A_ρ, which inverts the downtrend shown in the non-diffusive model. The ratio can widen up to ~5, as a function of ${\mathit{P}}_{\mathrm{t}}^{-1}\mathrm{\propto }\mathit{f}\mathit{/}\mathit{M}$, unveiling the presence of significant conductivity in the ICM and breaking any degeneracy in the interpretation of single spectra.

The real-space and projected maps carry important information:

• The ideal or poorly diffusive flows (f ≲ 10^-2) show complex filamentary and patchy density/entropy structures (similar in the core and outskirts), which are excited by the large-scale waves and later altered by hydrodynamical instabilities. The conductive models, however, show smooth maps due to the smearing of sharp features, which does not affect the turbulent eddies. Albeit sharing similar amplitude, velocities have uncorrelated phases with the tracer, hence a filament does not necessarily imply a high local velocity.

• The thermodynamic fluctuations can be described by log-normal distributions with weak non-Gaussian deviations, strengthening the role of the power spectrum.

• Synthetic X-ray images of velocity dispersion show that the forthcoming Astro-H (and Athena) will be able to well detect subsonic ICM turbulence. Using the broadening of the Fe XXV line, the detectable turbulent broadening is ≳200 km s^-1 (M_1D ≳ 0.13 for massive clusters), which can probe density perturbations on the order of a few percent and allows a precise calibration of the linear relation. The projected velocity maps (linked to the line shift) highlight the power stored in the large-scale motions, constraining the injection scale.

The analysis presented in this work shows the wealth of information that can be extracted from the ICM power spectrum. In Coma cluster, Schuecker et al. (2004) retrieve pressure fluctuations, which are mildly adiabatic and trace a Kolmogorov spectrum, in line with a M ~ 0.4 turbulent flow. In Gaspari & Churazov (2013, see also Churazov et al. 2012), we showed that the density spectrum arising from deep Chandra data of Coma is consistent with a similar level of turbulence (several 100s km s^-1), along with highly suppressed conduction (f ~ 10^-3). Sanders & Fabian (2012) found density fluctuations (≲8 percent) having a cascade shallower than Kolmogorov in the AWM7 cluster, which implies highly suppressed conduction and M ≲ 0.18. Being able to quickly convert between thermodynamic properties and gas motions through a simple linear relation, or being able to assess the plasma diffusivity through the spectral slope or the A_v/A_ρ diagnostic, is a powerful tool for both observational and theoretical study. The same analysis can be extended to other gaseous halos, such as massive galaxies and groups. Although current constraints are in its embryonic stage, we are beginning to understand the richness of information that the ICM power spectrum can convey. Future studies and observations (not only in the X-ray band but also via SZ maps) will be able to improve and exploit the full potential of the ICM power spectrum, helping us to probe the physics of the gaseous medium with high precision.

---

Acknowledgments

The FLASH code was in part developed by the DOE NNSA-ASC OASCR Flash center at the University of Chicago. M.G. is grateful for the financial support provided by the Max Planck Fellowship. We acknowledge the MPA, RZG, and CLS center for the availability of high-performance computing resources. D.N. and E.L. acknowledge support from NSF grant AST-1009811, NASA ATP grant NNX11AE07G, NASA Chandra grants GO213004B and TM4-15007X, and the Research Corporation. M.G. thanks A. Schekochihin, R. Sunyaev, S. Borgani, F. Brighenti, F. Miniati, D. Eckert, S. Molendi, X. Shi, E. Pointecouteau for helpful comments, and the anonymous referee for a highly positive feedback.

References

All Figures