key: cord-0580254-rx4q5j31
authors: Mohapatra, Rajsekhar; Federrath, Christoph; Sharma, Prateek
title: Turbulent density and pressure fluctuations in the stratified intracluster medium
date: 2020-10-23
journal: nan
DOI: nan
sha: cc5d43bb17af5a25d10f03f4e3122b63811d1603
doc_id: 580254
cord_uid: rx4q5j31

Turbulent gas motions are observed in the intracluster medium (ICM). The ICM is density-stratified, with the gas density being highest at the centre of the cluster and decreasing radially outwards. As a result of this, Kolmogorov (homogeneous, isotropic) turbulence theory does not apply to the ICM. The gas motions are instead explained by anisotropic stratified turbulence, with the stratification quantified by the perpendicular Froude number ($mathrm{Fr}_perp$). These turbulent motions are associated with density and pressure fluctuations, which manifest as perturbations in X-ray surface brightness maps of the ICM and as thermal Sunyaev-Zeldovich effect (SZ) fluctuations, respectively. In order to advance our understanding of the relations between these fluctuations and the turbulent gas velocities, we have conducted 100 high-resolution hydrodynamic simulations of stratified turbulence ($256^2times 384$---$1024^2times1536$ resolution elements), in which we scan the parameter space of subsonic rms Mach number ($mathcal{M}$), $mathrm{Fr}_perp$, and the ratio of entropy and pressure scale heights ($R_{PS}=H_P/H_S$), relevant to the ICM. We develop a new scaling relation between the standard deviation of logarithmic density fluctuations ($sigma_s$, where $s=ln(rho/left<rhoright>)$), $mathcal{M}$, and $mathrm{Fr}_{perp}$, valid till $mathrm{Fr}_perpll1$:~$sigma_s^2=lnleft(1+b^2mathcal{M}^4+0.10/(mathrm{Fr}_perp+0.25/sqrt{mathrm{Fr}_perp})^2mathcal{M}^2R_{PS}right)$, where $bsim1/3$ for solenoidal turbulence driving studied here. We further find that logarithmic pressure fluctuations $sigma_{(ln{P}/left)}$ are independent of stratification and scale according to the relation $sigma_{(ln{bar{P}})}^2=lnleft(1+b^2gamma^2mathcal{M}^4right)$, where $bar{P}=P/left$ and $gamma$ is the adiabatic index of the gas.

Turbulence and buoyancy are concurrent in several geophysical and astrophysical systems -the physics of stratified turbulence governs ocean currents and atmospheric turbulence on the earth and other planets, radiative and convective zones in the atmospheres of the sun and other stars, gas motions in hot gaseous haloes of galaxies (the circumgalactic medium or CGM), galaxy groups and clusters (intracluster medium or ICM) (Stein 1967; Goldreich & Keeley 1977; Loewenstein & Fabian 1990; Sarazin et al. 1992; Rudie et al. 2012; Parmentier et al. 2013; Skoutnev et al. 2020 ). Here we focus on stratified turbulence relevant to the ICM.

ICM refers to the gas that pervades the region between galaxies in a cluster. It is mostly composed of the hot X-ray emitting gas, ★ E-mail: rajsekhar.mohapatra@anu.edu.au † E-mail: christoph.federrath@anu.edu.au ‡ E-mail: prateek@iisc.ac.in with temperatures ranging from 10 7 -10 8 K, although a filamentary colder phase has also been detected in many clusters (Cowie et al. 1983; McDonald et al. 2010; Simionescu et al. 2018; Vantyghem et al. 2019; Olivares et al. 2019) . It is moderately stratified (with Richardson number Ri 10 or Froude number Fr ⊥ 0.1, refer to figure 1 in Mohapatra et al. 2020 ) and the gas is in rough hydrostatic equilibrium with the gravitational profile set by the dark matter halo. Turbulence in the ICM plays an important role in the gas dynamics and evolution. Turbulent energy dissipation on viscous scales and subsequent heating of the ICM, turbulent mixing of hot and cold phases of gas (Kim & Narayan 2003; Banerjee & Sharma 2014; Hillel & Soker 2020 ) may also play a key role in the gas thermodynamics, by preventing the runaway cooling of the ICM core (Zhuravleva et al. 2014a) . Anisotropy in turbulent eddies may be used to probe the orientation of ICM magnetic fields (Hu et al. 2020) . In cluster outskirts, estimating the turbulent pressure support of the gas is important to get an unbiased estimate of the halo mass that is required for cosmology with clusters (Schuecker et al. 2004 ; George et al. 2009; Bautz et al. 2009; Cavaliere et al. 2011; Nelson et al. 2014; Biffi et al. 2016; Angelinelli et al. 2020 ). However, direct measurements of the ICM turbulent gas velocities (Hitomi Collaboration 2016) are still a few years away (XRISM 1 and Athena 2 ), after the early mission end of the Hitomi satellite. Recently, observers have relied on several indirect methods to estimate gas velocities in the hot phase, such as relating X-ray surface brightness fluctuations to turbulent velocity fluctuations (Zhuravleva et al. 2014a (Zhuravleva et al. , 2015 (Zhuravleva et al. , 2019 , relating Sunyaev-Zeldovich effect (SZ) observations to turbulent pressure fluctuations (Zeldovich & Sunyaev 1969; Khatri & Gaspari 2016; Mroczkowski et al. 2019) , measuring cold-gas velocities (e.g., using the H line) and relating this to the hot phase velocity . Simionescu et al. (2019) provide a detailed review of all these different methods of measuring gas velocities in the ICM. It is important to develop accurate and robust scaling relations between the ICM hot-gas velocities and these observables.

Many recent theoretical studies have focused on the scaling of density fluctuations with the rms Mach number (M) (Gaspari et al. 2014; Zhuravleva et al. 2014b; Nolan et al. 2015; Mohapatra & Sharma 2019; Shi & Zhang 2019; Grete et al. 2020; Mohapatra et al. 2020) . Some of these studies have ignored the effect of gravitational stratification or lack a detailed parameter scan of the range of Fr ⊥ and M relevant for the ICM. The relation between density fluctuations and turbulent velocities is also seen to depend on the equation of state of the gas (Federrath & Banerjee 2015) and the adiabatic index (Nolan et al. 2015) , and whether gas cooling is included (Mohapatra & Sharma 2019; Grete et al. 2020) .

Several fluid mechanics studies (Bolgiano 1959 (Bolgiano , 1962 Carnevale et al. 2001; Lindborg 2006; Brethouwer & Lindborg 2008; Herring & Kimura 2013; Kumar et al. 2014; Rorai et al. 2014; Feraco et al. 2018; Alam et al. 2019 ) discuss the theory of stratified turbulence in the context of planetary atmospheres and oceans. In these studies, turbulence is driven perpendicular to the direction of gravity, whereas turbulence in the ICM is driven more isotropically by active galactic nuclei (AGN) jets and galaxy mergers (Churazov et al. 2002 (Churazov et al. , 2003 Omma et al. 2004) . They also mainly focus on the scaling of velocity and passive scalar spectra, intermittency and velocity anisotropy in the strong stratification limit.

In our previous study (Mohapatra et al. 2020 ) (hereafter referred to as MFS20), we performed stratified turbulence simulations with a fixed M and scanned the parameter space of weakly and moderately stratified turbulence. We also proposed a scaling relation between density fluctuations, Richardson number (Ri) and M for this regime. We found that density fluctuations also depended on a third parameter, namely the ratio between pressure and entropy scale heights ( = / ). We found that for Ri 10 (only moderately stratified), these numbers are sensitive to the turbulent driving length scale , as Ri ∝ 2 . Here we scan the parameter space of these three parameters: M, the transverse Froude number Fr ⊥ (Fr ⊥ ≈ 1/ √ Ri for Fr ⊥ 1) and through 100 simulations, extending into the strongly stratified regime (Fr ⊥ 1). The paper is organised as follows. In section 2 we briefly describe our setup and methods, present our results and their interpretations in section 3, compare our results with the literature and discuss the caveats of our work in section 4, and conclude in section 5. 1 https://global.jaxa.jp/projects/sas/xrism/ 2 https://www.the-athena-x-ray-observatory.eu/

We model the ICM as a fluid using compressible Euler equations and ideal gas equation of state. We implement gravity and turbulent forcing as additional source terms in the momentum and energy equations. We solve the following equations:

where is the gas mass density, v is the velocity, = /( ) is the pressure (we use the ideal gas equation of state), F is the turbulent acceleration that we apply, Φ is the gravitational potential, g = −∇Φ is the acceleration due to gravity, = 2 /2 + /( − 1) is the sum of kinetic and internal energy densities, is the mean particle weight, is the proton mass, is the Boltzmann constant, is the temperature, and = 5/3 is the adiabatic index.

We choose −ẑ to be the direction of the gravitational field, and pressure and density to have scale heights and , respectively. Thus, the initial pressure and density profiles are given by

We work with dimensionless units and choose 0 = 1 and 0 = 0.6, so that ,0 = √︁ 0 / 0 = 1. Since we start with the gas in hydrostatic equilibrium, the initial density, pressure, and , are related by

Hence is set as

This gives us the condition for the entropy scale height (≡ 1/[d ln /d ]), given by

This condition is satisfied for all our simulations, which locally mimic the stably stratified ICM.

When a parcel of gas in a stably stratified medium is displaced from its original position, it oscillates with a frequency defined as the Brunt-Väisälä (BV) frequency, given by

We define the turbulent time scale as ℓ ⊥ / ⊥ , where ℓ ⊥ is the integral scale, defined as

Here ⊥ is the velocity power spectrum perpendicular to the direction of g.

Here v ⊥ = ( , , 0) denotes the components of the velocity field perpendicular to the direction of gravity. The perpendicular Froude number Fr ⊥ is the ratio of these two time scales (see e.g., chapter 14 in Davidson 2013), given as

We can use the approximation ℓ ⊥ ≈ driv , where driv is the driving length scale of the turbulence. The parallel Froude number Fr is defined as

with being the integral scale parallel to the direction of gravity and is the velocity power spectrum parallel to the direction of g. Notice that the transverse velocity ( ⊥ , and not ) is used in the above expressions as peaks at scales smaller than ⊥ . In this study, we use Fr ⊥ to quantify the relative strength of stratification to turbulence 3 . The scale-dependent Froude numbersFr ⊥ (l ⊥ ) and Fr (l ), for an eddy of size (l ⊥ andl ), and perpendicular velocitỹ ⊥ are defined as

Fr ⊥ 1 indicates strong stratification (in this regime, Fr ≈ 1, see Billant & Chomaz 2001) and Fr ⊥ 1 denotes weak stratification. Strongly stratified turbulence transitions into weakly stratified turbulence at the Ozmidov length scale ℓ , which is defined as the scale on which the relative strengths of buoyancy and turbulence terms become equal. SinceFr ⊥ (ℓ ) = 1,

where = 3 ⊥ /ℓ ⊥ is the kinetic energy transfer rate and is assumed to be a constant just as in conventional turbulence. For weak and moderately stratified turbulence (Fr ⊥ 1), Fr ⊥ ≈ (1/Ri) 0.5 . Similar to MFS20, density, pressure and velocity are normalised to construct dimensionless variables, such that¯= / ( ) ,¯= / ( ) and M = / rms , where ( ) and

( ) are the average density and pressure at a − slice, respectively, is the amplitude of velocity, and ≡ ( / ) 1/2 is the local speed of sound. The potential energy per unit mass is defined as

where =¯/ is the strength of density fluctuations expressed in velocity units, and

is the ratio of pressure and entropy scale heights. The kinetic energy is defined as = 2 /2, where is the magnitude of the fluctuating velocity.

We evolve the Euler equations (1a to 1c) using the hydrodynamic version of the HLL5R Riemann solver (Bouchut et al. 2007 (Bouchut et al. , 2010 Waagan et al. 2011 ) in a modified version of the FLASH code (Fryxell et al. 2000; Dubey et al. 2008) , version 4. Our setup is the same as in MFS20 -we use a uniformly spaced 3D grid, with a box size = = 1 and = 1.5, centred at (0, 0, 0).

We use periodic boundary conditions for all variables (density, pressure and velocity) along the and directions. Along the direction, we use reflective boundary conditions for velocity and Dirichlet boundary conditions for pressure and density with the guard cells filled according to eq. (2a) and (2b), respectively. In MFS20, we used reflective boundary conditions along the direction, which led to a hydrostatic instability at the direction boundaries (due to an inverted pressure gradient at these boundaries). This instability was stronger for strongly stratified turbulence simulations and led to anomalous sound waves moving in the direction, starting from the boundaries (also seen in Shi & Zhang 2019, hereafter SZ19). The Dirichlet boundary conditions used here let us avoid this instability and allow us to extend our study to include the strongly stratified turbulence limit (down to Fr ⊥ ∼ 0.05). We restrict our analyses to a cube of size 1 centred at (0, 0, 0), with boundaries at (±0.5, ±0.5, ±0.5), to avoid any other anomalous effects near the direction boundaries. We run our simulations with shallow density profile ( > 1) on grids with 256 2 ×384 resolution elements, and the simulations with steep density profile ( ≤ 1) or simulations with high rms Mach number (M ≈ 0.4) on grids 512 2 × 768 resolution elements. We also run four strongly stratified simulations at resolution 1024 2 × 1536 for numerical convergence checks.

We follow the same spectral forcing method as in MFS20. We use the stochastic Ornstein-Uhlenbeck (OU) process to model the turbulent acceleration F with a finite autocorrelation time scale turb (Eswaran & Pope 1988; Schmidt et al. 2006; Federrath et al. 2010 ). We inject power as a parabolic function of | |, for 1 ≤ | | ≤ 3 (note that we have dropped the wavenumber unit 2 / ). The power peaks at | | inj = 2, i.e., driv = /2. For ≥ 3, turbulence develops self-consistently. We set turb = driv / , where is the standard deviation of the velocity on driv . Our driving is solenoidal (zero divergence). For further details of the forcing method, refer to section 2.7 of MFS20 and section 2.1 in Federrath et al. (2010) .

We also use the same window function, ( ), on the acceleration field, as in MFS20, such that F decays to zero near the boundaries in the direction, given by

Note that ( ) = 1 inside the analysis box (| |, | |, | | ≤ 0.5). So it only serves to exponentially decrease the acceleration amplitudes close to the boundaries.

We have conducted 100 simulations, scanning Fr ⊥ between 0.05-12.0, M between 0.01-0.4 and between 0.33-2.33, covering the parameter range relevant for the stratified ICM. We have three input parameters ( , and the acceleration field) that we vary in our various simulations to scan the range of interest in M, Fr ⊥ , and .The bin Mach number, M bin , indicates the rms Mach number that the acceleration field would produce in a homogeneous isotropic setup. We use these bins to separate our runs with different Mach numbers. The actual M can be slightly different from M bin , especially in models with steep temperature profiles ( , 0.25), as can vary by an order of magnitude with height. By definition, depends only on / (c.f., Equation (7)). Fr ⊥ depends on all three input parameters, roughly Fr ⊥ ∝ M bin / √ (as ⊥ ∝ M bin and 2 ∝ 1/[ ]). In table 1, we provide a compressed list of the simulation models, grouped under their common M bin and . We indicate the range of , the resolution, M, Fr ⊥ , and 2 ( is the standard deviation of = ln( / ) in the analysis box) for each of these groups. In table B1, we have expanded each grouping and list all of these parameters individually for each of the 100 simulations.

In order to take into account the variation in Fr ⊥ and M along due to steep temperature profiles, and to maintain uniformity during post processing among all our simulations, we divide the the analysis data cube (see definition in section 2.4.1) into four slabs along the direction (−0.5 ≤ < −0.25,−0.25 ≤ < 0, 0 ≤ < 0.25 and 0.25 ≤ ≤ 0.5). This means that for each simulation, we have four sets of data-points ( , ln(¯) , M, and Fr ⊥ ), corresponding to each of these four slabs.

In the presence of significant turbulent pressure (∝ local M 2 ), comparable to the thermal pressure, our Dirichlet+reflective boundary conditions break hydrostatic equilibrium. This leads to fluctuations in the computational domain whose amplitude increases for steeper pressure profiles, since the thermal pressure is small at the = 0.75 boundary. These fluctuations originate at the upper boundary, but they are confined close to the boundary itself, since strong stratification prevents them from travelling to lower . Therefore, as a precaution, we ignore the upper two slabs and only use the lower two slabs in our analysis of simulations that have steep pressure profiles ( < 0.2). All our simulations run for a total of 16 eddy turnover times ( eddy ≈ turb ) on the driving length scale. The simulations reach a steady state between 3-6 turb . We analyse turbulence from 6 turb to 16 turb , for a total duration of 10 turb , for statistical averaging.

Now we describe the results of our simulations and discuss their possible interpretations. We also compare our results against¯--M relations in other stratified turbulence simulations of the ICM.

In fig. 1 , we compare 6 representative simulation models (3 different models with high to low Froude number, from top to bottom, and 2 different Mach numbers, left versus right). Each panel shows the column density fluctuationsΣ = ∫¯d − 1, with the projected velocity field superimposed as vectors, where denotes the line of sight (LOS) direction (the absolute column density Σ = ∫ d is shown in the insets to provide a sense of the strength of the stratification). We have chosen as the LOS, which is perpendicular to the direction of stratification. For small density fluctuations (¯< 1), theΣ plots provide a sense of comparison (refer to section 3.2 of MFS20) to X-ray surface brightness fluctuations in Zhuravleva et al. (2014a) , which have been used to reconstruct turbulent velocities of ICM gas.

For this figure, we have chosen 6 representative simulations, with two different M and three different Fr ⊥ , such that we roughly present the extremes of these two parameters. We have plotted them such that M is approximately constant along a column, Fr ⊥ is approximately constant along a row and is a constant for all. For the Σ plots shown in the insets, we have used a separate log-scale colourbar for each row. The steepness of the Σ profile is inversely proportional to . For these plots, Fr ⊥ ∝ M , as decreases between the left and right panels, and also decreases from the top to the bottom panels. The column density fluctuations (Σ ) increase with both M and Fr ⊥ , but are more sensitive to the changes in M.

In the top and middle rows (weakly and moderately stratified turbulence), we observe that the eddies are roughly circular and hence the velocity field is roughly isotropic (similar to figure 3 of MFS20). However, for strongly stratified turbulence shown in the bottom row, the eddies become flatter in the direction and turbulence forms layered stratified structures. We also observe the correlation betweenΣ and in the middle row, where upward velocity arrows ( > 0) are associated with regions whereΣ > 0 (shown in red), and downward velocity arrows ( < 0) are associated with regions whereΣ < 0 (shown in blue). This positive correlation implies that some of the kinetic energy is converted into buoyancy potential energy.

In homogeneous idealised turbulence, the velocity field is expected to be isotropic and to follow a nearly Gaussian distribution. However, external fields such as gravity (lower panel of fig. 4 in MFS20; SZ19) and magnetic fields (Federrath 2016; Beattie et al. 2020) can induce strong anisotropies in velocity fields. In our strongly stratified simulations, we expect the component of velocity to be the most affected by the strength of the stratification.

Understanding the variation in as a function of stratification is of key significance to our study, since in MFS20 we showed that the additional density fluctuations introduced by buoyancy effects (¯b uoy ) are correlated to . We also derived a scaling relation for¯b uoy , assuming the velocity field to be isotropic, i.e., 2 ≈ 2 /3. However, this assumption breaks down for strongly stratified . The Froude number in each row decreases from ∼ 10 (top panels), to ∼ 1 (middle panels), to ∼ 0.2 (bottom panels). The colourbar for theΣ panels (shown below the sub-panels) has a linear scale. Note that each row of the Σ insets has its own colourbar, which is in log scale. A movie of the time evolution of these representative simulations is available at this youtube link. 

Notes: Column 1 shows the simulation name. The numbers following 'M' and ' ' are the binned Mach number (M bin ) and the ratio of pressure to entropy scale heights / in the simulations, respectively. In columns 2 and 3 we list the range and the resolution range. These parameters are defined in eqs. (2a) and (2b) and section 2.5. The default resolution of all the weak stratification runs ( > 1) is 256 2 × 384. The runs with stronger stratification ( < 1) or M bin = 0.4 are run with 512 2 × 768 resolution elements. We also run four simulations at a resolution of 1024 2 × 1536, for convergence checks and four strongly stratified simulations at M = 0.01 and resolution 512 2 × 768. Column 4 lists the actual rms M range, which can be different from the targeted M in strongly stratified simulations. In column 5, Fr ⊥ refers to the mean perpendicular Froude number of the simulations (see equation 6d). Column 6 shows 2 , the standard deviation of = ln¯squared. All quantities (M, Fr ⊥ , ) were averaged over 10 turbulent turnover times, for 6 ≤ / turb ≤ 16. A more detailed version of this table, which lists each of the total of 96 simulations used here, is available in the appendix table B1. turbulence (Fr ⊥ 1), since the velocity field becomes strongly anisotropic and motion is confined to layers perpendicular to the direction of gravity, as we see in the lower panels of fig. 1 .

In order to investigate this further and understand the scaling of with Fr ⊥ , we show the ratio of rms to ⊥ (defined in eq. (6c)) in fig. 2 . For weak stratification (Fr ⊥ 1), the velocity is roughly isotropic ( / ⊥ ≈ 1) 4 , perpendicular and parallel velocities being roughly the same. As we move from right to left in fig. 2 , in the moderately stratified turbulence regime (Fr ⊥ ≈ 1), the ratio starts decreasing slowly with decreasing Fr ⊥ . This is expected, as more of the direction kinetic energy gets converted into buoyancy potential energy with increasing strength of stratification. As we move further left in the plot to the strong stratification limit (Fr ⊥ 1), the ratio shows a sharp decrease with decreasing Fr ⊥ , with / ⊥ ∝ Fr 0.6 ⊥ . For strongly stratified turbulence, buoyancy dominates on large scales, for ℓ <l < driv (defined in eq. 6j), till the relative strengths of buoyancy and turbulence terms become equal. For ℓ < ℓ , the turbulence transitions to weakly stratified turbulence. Thus, unlike isotropic turbulence, is set by˜ (ℓ ) instead of the integral scale velocity . Forl > ℓ , we may use the Boussinesq approximation (changes in density are small relative to the mean density), for which eq. (1a) reduces to ∇ · v = 0. Then the ratiõ

using eqs. (6h) and (6i), where we assumeFr ∼ 1 for strong stratification and that the vertical velocity fluctuations peak at the Ozmidov scale and not the driving scale. This gives us

using eqs. (6j) and (8) (see example 14.2 in Davidson 2013). We observe / ⊥ ∝ Fr 0.6±0.1 ⊥ in our strongly stratified simulations, which roughly agrees with the theoretical prediction. Some deviations from the theoretical prediction may arise because the low-Froude simulations are not in the limit of Fr ⊥ 1, but around Fr ⊥ ∼ 0.1.

In this subsection, we discuss the variation of density and pressure fluctuations as a function of our simulation parameters -M, Fr ⊥ 10 −1 10 0 10 1 Fr ⊥ 10 0 

. The scaling relations between these quantities are important for obtaining ICM gas velocity estimates from X-ray and SZ observations. The X-ray surface brightness fluctuations and SZ effect fluctuations are used to calculate the amplitude of density and pressure fluctuations, respectively, which are further used to calculate velocity fluctuations using these relations (Zhuravleva et al. 2013 (Zhuravleva et al. , 2014b Khatri & Gaspari 2016) .

In fig. 3 , we show the density and pressure fluctuations squared ( 2 and 2 ln¯) versus M for all our simulations. We also present approximate fits based on the scaling relations that we propose in the next section 3.3.2. Clearly, 2 increases with increasing M, stratification strength (decreasing Fr ⊥ ) and

. In comparison, 2 ln¯o nly increases with M and is independent of Fr ⊥ and . As we discussed in section 3.5 of MFS20, the density fluctuations are comprised of two components, which can be written as the sum of an un-stratified and a stratified turbulence component:

where¯2 turb scales as 2 M 4 (Mohapatra & Sharma 2019; MFS20), and is the turbulence driving parameter ). The parameter = 1/3 for solenoidal turbulence, as we use for the present set of simulations. Our fits in fig. 3 show this dependence on M for Fr fit 1 (light bronze coloured fits). We also showed that¯t urb corresponds to adiabatic density fluctuations, so the corresponding pressure fluctuations¯2 scale as 2¯2 turb ∝ M 4 for For moderate stratification (Fr ⊥ ∼ 1) or low M turbulence, thē buoy term dominates, but it corresponds to the isobaric motions of isotropic gas parcels, which have zero contribution to the net pressure fluctuations. Hence¯still scales as¯t urb and shows a M 4 variation throughout as seen in the lower panel of fig. 3 . This scaling seems to hold even in the strongly stratified turbulence limit, which means that the nature of¯b uoy is isobaric even for Fr 1. This behaviour is expected, since a rising/falling parcel of gas with subsonic velocity is always in pressure equilibrium with its immediate surroundings (such that¯b uoy = 0). Thus, the expression for logarithmic pressure fluctuations becomes 2 ln¯= ln 1 + 2 2 M 4 .

In the lower panel of fig. 3 , we fit the data to eq. (11) using the fitting tool LMfit (Newville et al. 2016 ). The results are in good agreement with our expectations.

In MFS20, we also showed that¯2 buoy increases with M as approximately M 2 , for constant Fr ⊥ 1 (or Ri 1). The motions in the direction associated with¯b uoy are strongly constrained for Fr ⊥ 1 and because of the large energy cost, BV oscillations (with small displacement in the direction) dominate over turbulence in the vertical direction. But the scaling with M still holds in this limit, as is seen in the dark bronze fits in the upper panel of fig. 3 .

In this subsection, we discuss the scaling of density and pressure fluctuations with the stratification parameters Fr ⊥ and . We know that the¯2 turb and¯2 buoy terms scale as M 4 and M 2 , respectively. In order to make comparisons between different Fr ⊥ and easier, we normalise 2 and 2 ln¯t o 2 ,bin and 2 ln¯,bin , respectively, given by 

In the weak and moderately stratified turbulence limit (Fr ⊥ 1), we simplify this expression further and showed in MFS20 that

where 2 is a fitting parameter. However, this simplification involved assuming the velocity field to be roughly isotropic, which is true for Fr ⊥ 0.5 (see fig. 2 ), not for strongly stratified turbulence. In the limit of strongly stratified turbulence (Fr ⊥ 1), the direction motions are heavily suppressed by buoyancy and the velocity field is no longer isotropic, as we showed in section 3.2. In this limit, 

using eqs. (3a), (6a) and (9). Interpolating between the two asymptotic expressions for¯b uoy in eq. (13b) and eq. (13c), we get

where 1 and 2 are fitting parameters. The combined expression for the net density fluctuations is then given bȳ

Here we substituted 2 = ln 1 + 2 . This is valid for log-normal distributions of¯but it also holds approximately for non-log-normal distributions with small 2 (see Appendix 6). Thus, we find that introducing stratification has two main effects. For Fr ⊥ 0.5, due to the existing density gradient, the turbulent motions in the direction produce higher density fluctuations. This happens simply because when a gas parcel moves along the direction, it only attains pressure equilibrium at that height, but still has a density contrast with respect to its surroundings.On assuming isotropic gas velocities, which is roughly valid for Fr 0.5 (see fig. 2 ), we obtain¯2 buoy = 2 M 2 /Fr 2 ⊥ (eq. 13b). However, on further increasing the stratification (Fr ⊥ 0.5), the turbulence becomes anisotropic as gas motions along the direction are suppressed, with the kinetic energy along the direction being converted into buoyancy potential energy. The turbulent eddies flatten along the direction and become pancake-like (see lower panels of fig. 1 ). In this limit, the motion along the direction is best described by BV oscillations. The amplitude of these oscillations is proportional to , which decreases with decreasing Fr ⊥ for constant M (see fig. 6 ). Substituting the Fr ⊥ dependence of the ratio / ⊥ , we obtain¯2 buoy = 2 M 2 Fr ⊥ (eq. 13c). The general expression for the dependence of¯2 buoy (eq. 13d) is thus an interpolation between these two forms. For even stronger stratification (Fr ⊥ 0.001), we expect the amplitude of¯2 buoy to decrease below¯2 turb , and the net density fluctuations would again be given by¯2 ∼ 2 M 4 , similar to the unstratified subsonic turbulence scaling. One can interpret this limit as 2D subsonic turbulence, with M representing the 2D Mach number, since ⊥ .

Here we derive the fitting parameters 1 and 2 , such that they describe the variation in density fluctuations for all of our 96 simulations simultaneously. We use the fitting tool LMfit with eq. (13f) as the fitting function, M, Fr ⊥ and as independent variables, and 2 as the dependent variable. We obtain 2 1 = 0.10 ± 0.01 and 2 = 0.25 ± 0.02.

In order to present all our data and the fitting function together, we scale all density and pressure fluctuations to M scaled = 0.25 and = 2.33. These are the same parameter values we used in all simulations of MFS20 and hence they provide a direct comparison between both of our studies. We construct two new scaled quantities, given by 

In fig. 5 , we show these two quantities as a function of Fr ⊥ . We also show the fitting function eq. (13f) for M = 0.25 and = 2.33. Thus, our scaling relation for 2 works well for the entire range of parameter space we have scanned using the 96 simulations. The value of 2 ln¯,scaled is roughly constant, as expected from eq. (11). We expect the variation to be mostly due to two reasons both leading to inaccurate estimates of M, and these variations are amplified in the plot since 2 ln¯,scaled ∝ M 4 . Firstly, the box heats up in larger M bin simulations (especially M bin = 0.25 and 0.4), which leads to an increased speed of sound and a decreased Mach number. Secondly, due to the steep temperature profiles for simulations with , < 0.2, ( ≠ ), the speed of sound significantly varies along the direction, whereas we drive an isotropic velocity field. This leads to a strong variation in M along the direction, which may cause the differences from the scaling relations for Fr ⊥ 0.5. 

In this subsection we discuss the importance of pressure fluctuations for estimating turbulent gas velocities in the ICM. In the previous subsections, we showed that density fluctuations are sensitive to the parameters Fr ⊥ and , whereas pressure fluctuations depend only on M. We found that the relation 2¯= 2 2 M 4 describes all our simulations very well (see lower panels of fig. 3 and fig. 5 ). This happens because the stratified turbulence component of density fluctuations is isobaric (see figure 8 of MFS20). In figure 14 of Mohapatra & Sharma (2019) , we also showed that pressure fluctuations do not significantly depend on whether radiative cooling is included or not. While density fluctuations are sensitive to the thermodynamic and stratification parameters, pressure fluctuations are independent of these. Thus, compared to density fluctuations, we expect pressure fluctuations to provide a robust estimate of turbulent velocities.

Hot electrons in the ICM, on average, up-scatter the cosmic Notes: All these simulations have grid resolution 512 2 × 768. The columns (2) -(5) have their usual meanings. For these runs, / = 1. microwave background (CMB) photons through inverse Compton scattering, and the change in the CMB brightness temperature is proportional to the ICM electron pressure integrated along the line of sight. This effect is known as the thermal Sunyaev-Zeldovich (tSZ) effect (see Mroczkowski et al. 2019 , for a review). Although it is possible to recover turbulent gas velocities from the resolved observations of the tSZ effect (Khatri & Gaspari 2016) , the present observations suffer from a lack of angular resolution (e.g., up to a few 100 kpc for the Coma cluster). The tSZ observations also suffer from projection effects (more than X-ray surface brightness which is proportional to 2 ), which makes the measured pressure fluctuations even smoother than the density fluctuations. Many of these limitations are expected to be addressed by future facilities (see section 6.2 of Mroczkowski et al. 2019) . The relative robustness of the relation between pressure fluctuations and turbulent velocities (unlike density fluctuations which depend on stratification and thermodynamic parameters) provides a strong motivation for high angular resolution SZ observations. simulation parameters. Since these simulations have = , the sound speed is uniform throughout the domain and the turbulent pressure is negligible compared to the thermal pressure (turbulent pressure ∼ 10 −4 thermal pressure). Since , > 1.0, the density and pressure gradients are not very steep. Thus, these simulations do not suffer from fluctuations near the boundaries mentioned in sections 2.4.1 and 2.5.1.

In fig. 6 we show the ratio between and ⊥ for these simulations. The ratio follows the Fr 0.6 scaling. In fig. 7 we show 2 ,scaled and 2 ln(¯),scaled , and we find that they also follow the scaling relation in the low-Fr ⊥ regime-2 ,scaled decreases with decreasing Fr ⊥ and 2 ln(¯),scaled is independent of Fr ⊥ .

In this subsection, we test our scaling relation against data from SZ19, who study decaying turbulence in a stratified medium, using a setup similar to ours. To construct 2 , we use the approximation 2 = ln 1 + 2 , and we scale these quantities to 2 ,scaled with M = 0.25 and = 2.33, using eq. (14a). We also found that our definitions of the integral scale ⊥ differ by a factor of 2 (see equation 3 in SZ19) which is the conversion factor from wavenumber space to real space. We have corrected for this by dividing Fr ⊥ in SZ19 by the same factor (2 ). ,scaled ) vs Fr ⊥ with the coloured (cyan, green and orange) datapoints taken from SZ19 and scaled (see the main text). We follow the same colour and marker scheme as in figure 8 of SZ19, where different colours represent different values of (which correspond to simulations with different levels of stratification). The open symbols are for ≤ 2/ and the filled symbols are for > 2/ . We show our data points from the top panel of fig. 5 in grey. We have also plotted data points from our turbulence decay simulations in red and purple. The dashed line shows the best fit to SZ19's data and the solid line shows the best fit to our data. There seems to be a small scaling discrepancy which becomes larger at smaller Fr ⊥ . We discuss more about this discrepancy and its possible sources in the main text.

In fig. 8, we show 2 ,scaled as a function of Fr ⊥ , using the same colour scheme (cyan, green and orange) and markers as figure 8 in SZ19. We have greyed out our data points in the background. SZ19 study decaying turbulence, and therefore all variables are time ( ) dependent in their study. The open symbols in fig. 8 are for ≤ 2/ and the filled symbols are for > 2/ . The dashed line is the best fit to the scaled SZ19 data and the solid line is the best fit to our data.

We find that 2 ,scaled for SZ19's data shows a similar dependence on Fr ⊥ as our scaling relation predicts. However, the exact dependence of 2 ,scaled on Fr ⊥ is somewhat different, i.e., we see that 2

,scaled is higher by almost a factor of 2-4 for Fr ⊥ 0.5.

We find that the fitting parameters 2 1 and 2 (see eq. (13f)) take the values 0.25 ± 0.02 and 0.19 ± 0.02 respectively. Despite the difference in the fit parameter, the functional form of our relation provides a good fit to both our and SZ19's data. The open data points for ≤ 2/ in SZ19 correspond to the moderate-weakly stratified turbulence regime (Fr ⊥ 0.5). In this limit,

For a given stratification profile (fixed , and corresponding ) and ≤ 2/ ,¯2 is independent of M, as seen in figure 8 of SZ19. This has significant implications for observational studies, which infer turbulent velocities from surface brightness fluctuations, which are in turn caused by density fluctuations. For > 2/ , Fr ⊥ 0.5, where the 2 dependence on Fr ⊥ is weak near the peak of the proposed scaling relation and hence¯2 ∝ M 2 .

In order to directly compare against SZ19, we conduct two . The ratio of density fluctuations to velocity power spectra, 2 , defined in eq. (16), as a function of the wavenumber , for our high-resolution, strongly stratified simulations. In this strongly stratified limit (which does not necessarily apply to the ICM in general; see Fig. 1 in MFS20), 2 depends strongly on and weakly on , M and Fr ⊥ .

sets of turbulence decay simulations with M bin = 0.05 and 0.25, = = 1.0 with grid resolution 512 2 × 768. Once turbulence reaches a roughly steady state (after ≈ 3 eddy ), we switch off external driving. We simulate five statistically similar instances of decaying turbulence by changing the random number seed for the turbulent driving.

We calculate the average value of 2 over time intervals of 2 eddy in the decay phase. We further average these values across the five runs. The red and purple circles in fig. 8 represent 2 ,scaled for these decaying turbulence simulations. Both the purple and red data points agree with our proposed scaling relation. It is worth noting that 2 ,scaled in SZ19's data becomes independent of Fr ⊥ for Fr ⊥ 0.3, whereas 2 ,scaled in both our steady-state turbulence and decaying turbulence simulations decreases with decreasing Fr ⊥ .

The differences in the fitting parameters may arise due to several factors-it is possible that some of the discrepancy in 2 ,scaled arises due to slight inaccuracies of the scaling that we have to apply in order to bring SZ19's data into our analysis parameter space, in particular because of the large (two orders of magnitude) difference in M, and the scaling involving terms proportional to M 2 (for the Fr ⊥ dependence) and proportional to M 4 (for the turbulence term in the relation). Finally, the very low Fr ⊥ simulations in SZ19 may be affected by spurious sound waves as discussed in section 2.4.1. Our turbulence decay simulations do not suffer from this because of the new Dirichlet boundary conditions that we apply to the density and pressure. These simulations are at low rms Mach number and use an initially isothermal stratification profile ( = = 1.0). Hence they also do not suffer from the numerical fluctuations at the upper boundary discussed in section 2.5.1.

In the previous subsections we have investigated the scaling of rms density fluctuations with the rms Mach number, Fr ⊥ , and . Another quantity relevant for observations is the amplitude and spectral scaling of density fluctuations and velocity power spectra. This is ( 2 ) is averaged over the inertial range of turbulence (6 ≤ ≤ 15). The fit (solid line) shown is plotted using M = 0.25 and = 2.33, corresponding to the parameters used in MFS20. The ratio ( 2 ) roughly scales as 2 /M 2 and is independent of M in the strongly stratified turbulence limit. We also show the values of 2 scaled and Fr clusters from Zhuravleva et al. (2014b) in the olive coloured shaded region. The values seem to be slightly larger than predicted by our scaling relation at low Fr ⊥ . We discuss about this discrepancy in the main text because turbulent velocities are more difficult to measure directly in observations, and often the surface brightness fluctuations are used to infer properties of the turbulence (Zhuravleva et al. 2014a , 2015 and Simionescu et al. 2019 for a review).

The ratio between Fourier one-component density fluctuation amplitudes¯and Fourier one-component velocity amplitudes 1, is denoted by and it is related to the density fluctuations and velocity power spectra by 2 = (¯) 2 /( 1, / ) 2 ≈ 3 (¯)/ (M ),

where (¯) and (M ) are the density fluctuations and velocity power spectra. Density and velocity are normalised with respect to the stratification profile and the speed of sound, respectively. The ratio given by is an important parameter for the observational studies listed above. The Fourier transform of the normalised Xray SB fluctuationsSB are used to calculate¯and then 1, is calculated using eq. (16).

All the observational studies we mentioned above assume 2 ≈ 1, independent of , M and Fr ⊥ . A value of 2 ≈ 1 was calibrated using cosmological simulations in Zhuravleva et al. (2014b) and large cluster scale (1000 kpc) simulations in Gaspari et al. (2014) . In Mohapatra & Sharma (2019) , we found 2 ∝ M 2 through idealised box simulations of unstratified turbulence. We also noted that the amplitude of 2 increases by two orders of magnitude on including cooling. In subsection 3.9.3 of MFS20, we studied the variation of 2 for weakly and moderately stratified turbulence (Fr ⊥ 0.5, M ≈ 0.25 and = 2.33). We showed that 2 increases with increasing stratification and reaches 2 ≈ 1 for Fr ⊥ ≈ 0.7. Here we have extended this study to the strongly stratified turbulence regime (up to Fr ⊥ ≈ 0.1).

In fig. 9 , we show 2 as a function of wavenumber for our high-resolution, strongly stratified simulations. For reference, we show 2 for the unstratified turbulence run (Ri0, Fr ⊥ − → ∞) from MFS20. The quantity 2 is relatively flat with respect to , for all of these simulations, which implies that (¯) and (M ) follow similar spectral scaling. We also show 2 , the amplitude of 2 averaged over 6 ≤ ≤ 15 as a function of Fr ⊥ in fig. 10 for the high-resolution simulations (1024 2 × 1536) of this study and MFS20 5 . For the two simulations with ≠ 2.33, we show 2 scaled = 2 × (2.33/ ). We fit the data using the expression for 2 /M 2 from eq. (13e) and M = 0.25 and = 2.33, with 1 and 2 as fit parameters. The fitted factor of (3.4 ± 0.1) is close to the expected factor of 3 from the 3D to 1D conversion; see eq. (16). We have also plotted 2 scaled for cluster simulations from Zhuravleva et al. (2014b) in the olive coloured shaded region in fig. 10 , which are slightly higher than our values and scaling relation. The differences in values of 2 scaled may arise from the choice of averaging regions, in this case 500 kpc regions in the clusters of Zhuravleva et al. 2014b . In addition to this, since turbulence in these cosmological simulations is driven more naturally by galaxy mergers and in-falls, they may include compressive components (which are excluded in our simulations), increasing the value of the turbulence driving parameter and generating larger density fluctuations. The cosmological simulations may also suffer from a lack of resolution on small scales and have a higher numerical viscosity compared to our idealised box simulations.

The amplitude of 2 has the same dependence on M, Fr ⊥ , and , as the ratio of the rms fluctuations¯2/M 2 . Since¯2 buoy is the main component of¯2 for Fr ⊥ ≈ 0.1, 2 scales similar to¯2 buoy /M 2 . Hence 2 is roughly independent of M, weakly dependent on Fr ⊥ (since versus Fr ⊥ peaks at around Fr ⊥ ≈ 0.1), and linearly dependent on , in agreement with eq. (13d). Combining the results of this work with those of Mohapatra & Sharma (2019) and MFS20, we have studied the variation of 2 with M, Fr ⊥ , and , scanning the parameter space relevant for the ICM. We find that the ratio 2 is almost invariant with the wavenumber 6 . The amplitude of 2 varies as¯2/M 2 . Based on eq. (13e) and the fit in fig. 10 , we propose a new scaling relation for 2 : 2 = (3.4 ± 0.1) 2 M 2 + (0.11 ± 0.01)

For Fr ⊥ 1, 2 approaches 1 and is independent of M, which was the limit studied in Zhuravleva et al. (2014b) and Gaspari et al. (2014) . However, as we showed in figure 1 of MFS20 (also in figure 2 of SZ19), the ICM can have Fr ⊥ between 0.1-100. Thus, we suggest that observational studies use eq. (17) to obtain a more accurate estimate of turbulent velocities from density (surface brightness) fluctuations. Since Fr ⊥ 1 also marks the onset of large-scale density anisotropy (see lower panels of fig. 1 and also  fig. 10 in MFS20), one could use the peak value of 2 , where it varies slowly with Fr ⊥ for relating the ICM density and velocity fluctuations

In this section we discuss possible shortcomings of our work and possible methods to address some of them. 5 All simulations in MFS20 had = 2.33 and M ∼ 0.25. 6 Except for the heating and cooling simulations in Mohapatra & Sharma (2019) , where switching on radiative cooling led to a steepening.

The buoyancy Reynolds number Re buoy for stratified turbulence is given by Re ⊥ Fr 2 ⊥ , where Re ⊥ is the turbulent Reynolds number in the transverse direction. When Re buoy 1, viscous forces can be ignored on the integral scale ℓ ⊥ (Davidson 2013). The ratio Re ⊥ scales as 4/3 in numerical simulations where is the number of resolution elements (Frisch 1995; Haugen & Brandenburg 2004; Benzi et al. 2008; Federrath et al. 2011) . Hence for low Fr ⊥ simulations, Re buoy ≈ 4096 × Fr 2 ⊥ ≈ 10 for Fr ⊥ 0.05, using = 512. Thus, we approach the limit Re buoy ≈ 1 for our strongly stratified simulations, implying that viscous forces are almost of the same order as the inertial forces on the integral scale for Fr ⊥ ≈ 0.05. Using higher grid resolution is one of the ways to avoid this issue.

In section 2.5.1, we mentioned that the positive boundaries (at = 0.75) are unstable for steep pressure profiles ( < 0.25). This happens because for the steepest pressure profiles with ≈ 0.1, the initial pressure varies by more than a factor of 10 6 across the entire domain. These steep gradients could possibly make the code numerically unstable near the low-pressure positive boundary. Using smaller box sizes (so that pressure values at the positive box boundary are higher) and/or higher resolution along the direction is a possible solution to this problem. This method has been used by many fluid mechanics studies to study strongly stratified turbulence (Lindborg 2006; Brethouwer & Lindborg 2008) .

In this paper, we have only studied solenoidally-driven stratified turbulence. The driving parameter is also part of the scaling relation (Federrath et al. , 2010 . Compressively-driven turbulence is expected to generate larger density fluctuations and its effect is supposed to be captured by the driving parameter , which we have not varied in any of our simulations. In nature, turbulence is supposed to be a mixture of compressive and solenoidal mode. Galaxy infall and mergers could possibly drive some compressive modes (Churazov et al. 2003; Federrath et al. 2017) in the ICM. Effects of different levels of compressively-driven turbulence would be an interesting follow-up study.

We have ignored the effects of radiative cooling, thermal conduction, and magnetic fields in this work. Radiative cooling (Mohapatra & Sharma 2019; Grete et al. 2020 ) and thermal conduction (Gaspari & Churazov 2013; Gaspari et al. 2014) have been shown to affect the¯-M scaling relation. However, this additional physics is beyond the scope of the current paper and will be addressed in a follow-up study.

We have studied the parameter space relevant to subsonic, stratified turbulence in the ICM through idealised high-resolution hydrodynamic simulations. We have covered the parameter regime 0.05 < Fr ⊥ < 12.0, 0.05 ≤ M ≤ 0.4 and 0.33 ≤ ≤ 2.33 through 96 simulations. The main results of this study are as follows:

(i) We have extended the scaling relation between the rms density fluctuations (denoted in log-scale by ), the rms Mach number (M), the perpendicular Froude number (Fr ⊥ ), and the ratio between pressure and entropy scale heights ( ) to the strong stratification limit (Fr ⊥ 1). The new scaling relation is

2 . The density fluctuations increase with decreasing Fr ⊥ for Fr ⊥ 0.2, saturate and then decrease slowly for Fr ⊥ . We have shown that the density fluctuations in all of our 96 simulations follow this scaling relation. Our results also qualitatively agree with the turbulence decay study in SZ19.

(ii) We have also extended the scaling of pressure fluctuations to the limit Fr ⊥ 1. We find that ln¯i s independent of the stratification parameters and depends only on M. The scaling relation remains ln¯= ln 1 + 2 2 M 4 . Since pressure fluctuations are unaffected by stratification, they can be used to obtain fairly accurate estimates of turbulent velocities through tSZ observations.

(iii) The ratio between normalised one-dimensional Fourier density (¯) and velocity amplitudes ( 1, / ) approximately scales as¯/M and saturates to ≈ 1 at Fr ⊥ ≈ 0.2. The square of its amplitude is given by 2 = 3.4 2 M 2 + 0.11 (Fr⊥+0.19/ √ Fr ⊥)

2 . (iv) The ratio between velocity components parallel and perpendicular to gravity ( / ⊥ ) is roughly constant with respect to Fr ⊥ for Fr ⊥ 1, as expected for weak stratification. For 0.5 Fr ⊥ 2, the ratio weakly decreases with decreasing Fr ⊥ , and for Fr ⊥ 0.5, it varies as Fr 0.6 ⊥ .

This work was carried out during the ongoing COVID-19 pandemic. The authors would like to acknowledge the health workers all over the world for their role in fighting in the frontline of this crisis. We thank the anonymous referee for their useful comments, which helped improve this work. RM acknowledges helpful discussions with Piyush Sharda. RM thanks Xun Shi and Irina Zhuravleva for providing data for subsection 3.3.6 and section 3.4, respectively and useful discussions. CF acknowledges funding provided by the Australian Research Council (Discovery Project DP170100603 and Future Fellowship FT180100495), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). PS acknowledges a Swarnajayanti Fellowship from the Department of Science and Technology, India (DST/SJF/PSA-03/2016-17). We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. The simulation software FLASH was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago.

All the relevant data associated with this article is available upon request to the corresponding author.

Movies of projected density and density fluctuations of different simulations are available at the following links on youtube: ( 2 ) for all our runs. The symbol shape corresponds to , the colour to the Mach number, and the shading to Fr ⊥ . The data is fit well by 2 = ln 1 + 2 , as discussed in the last paragraph of section 3. This paper has been typeset from a T E X/L A T E X file prepared by the author.

Turbulence in Rotating, Stratified and Electrically Conducting Fluids

Astronomical Society of the Pacific Conference Series

The Multi-Messenger Astrophysics of the Galactic Centre

Lmfit: Non-Linear Least-Square Minimization and Curve-Fitting for Python

 Figure A1 shows the relation between 2 and 2 . The relation 2 = ln 1 + 2 holds for log-normal distributions of density and/or small density fluctuations (¯ 1) (Price et al. 2011) . We use this relation to explain the scaling of 2 based on the turbulent and buoyant components of 2 in our study. Table B1 provides a long version of table 1 , listing all the simulations, including all input parameters and relevant calculated/output variables.