© The Authors 2025

1. Introduction

Turbulence reigns as the most ubiquitous and yet most enigmatic face of fluid movement. It drives the transport of matter at nearly all scales from micro- to astrophysical scales (Dubrulle 2019). In situations where the characteristic speed of a fluid is much less than the speed of sound in the surrounding medium, the assumption of incompressible flow offers a useful simplification. Most current knowledge pertains to this category of fluid flows. In this context, turbulence is recognised as a paradigmatic out-of-equilibrium system where the energy injected at large scales is transferred across intermediate scales via a nonlinear cascade process, before being dissipated at the smallest scales. The characterisation and prediction of velocity fluctuations at different flow scales constitutes a fundamental pursuit in the field of fluid mechanics.

In astrophysical contexts, fluid flows are however best described as compressible and/or with significant density fluctuations, a characteristic that emerges across structures of vastly different cosmological scales, from the intergalactic medium, the interstellar medium, to accretion discs, and stellar interiors. Owing to their immense scale and mass, it is then more accurate to recognise compressible and/or variable-density turbulence as the prevailing state of fluid motion. In the interstellar medium (ISM), and particularly within cold, dense molecular clouds that serve as stellar nurseries, compressible turbulence plays an ambivalent role (see for example refs. Elmegreen & Scalo 2004; Mac Low & Klessen 2004; McKee & Ostriker 2007; Hennebelle & Falgarone 2012, and references therein). It contributes to the formation of dense clumps that triggers gravitational collapse, while simultaneously injecting turbulence that counteracts this collapse and inhibits star formation. A detailed understanding of compressible turbulence is thus crucial for the characterisation of the structure and dynamics of most astrophysical systems.

Compressible turbulence is distinguished from its incompressible counterpart by the interdependence of fluid motion and thermodynamic properties, manifesting as coupled fluctuations in velocity, pressure, and density (see e.g. ref. Rabatin & Collins 2023; Scannapieco et al. 2024, and references therein). Then, the thermodynamics of the fluid matters (Passot & Vázquez-Semadeni 1998; Banerjee & Galtier 2014; Sakurai & Ishihara 2024). This coupling also yields energy conversion processes, particularly between the kinetic energy and internal energy reservoirs, thereby enriching the complexity of the mechanisms at play (Galtier & Banerjee 2011; Aluie et al. 2012; Wang et al. 2018; Schmidt & Grete 2019). For astrophysical applications, additional forms of energy such as magnetic (e.g. Banerjee & Galtier 2013; Seta & Federrath 2021) and gravitational potential energy (e.g. Banerjee & Kritsuk 2018) must also be considered, further complicating the picture of energy exchange processes.

Observations indicate that the interstellar medium (ISM) is highly inhomogeneous (Myers 1978), with regions of low and high turbulence intensity often coexisting in proximity (Linsky et al. 2022). As a result, the local turbulent conditions within a sub-region of the ISM may differ significantly from those inferred from cloud-averaged properties. Rather than introducing further complexity by explicitly modelling this inhomogeneity, we adopt a simplified approach: we treat the ISM as a collection of locally homogeneous sub-regions, each characterised by distinct turbulent properties. Our aim is to systematically characterise the statistical signatures associated with different turbulent regimes. We focus in particular on the statistics of velocity and density. For simplicity, we also assume the ISM to be isothermal, an assumption generally justified in dense molecular clouds, where radiative cooling is highly efficient (Elmegreen & Scalo 2004; Glover et al. 2010; Krumholz 2014). Despite its simplicity, this framework offers a useful first step toward capturing the statistical diversity of ISM substructure.

In practice, homogeneous turbulence can be readily realised through direct numerical simulations (DNS) of the compressible Navier-Stokes equation (Pirozzoli & Grasso 2004; Donzis & Jagannathan 2013; Jagannathan & Donzis 2016; John et al. 2021; Sakurai et al. 2021; Sakurai & Ishihara 2023, 2024; Wang et al. 2011, 2017a, 2018; Seta & Federrath 2021) or implicit large eddy simulation (ILES) of the Euler equation in a periodic box (Kritsuk et al. 2007; Federrath et al. 2008; Schmidt et al. 2009; Federrath et al. 2010; Aluie et al. 2012; Pan et al. 2019a; Federrath et al. 2021; Rabatin & Collins 2023; Scannapieco et al. 2024; Beattie et al. 2025). DNS of the Navier-Stokes equation can now achieve relatively large Re although generally restricted to relatively low Ma, from the sub- to transonic regimes. In contrast, ILES of the Euler-equation are typically used to explore the supersonic regime but with no explicit control of the viscous cut-off.

A substantial body of work has now established a relatively coherent picture of the statistical properties of compressible turbulence. For instance, the Mach-number dependence of fluctuations of the density, velocity, and pressure, and their spatial derivatives are well documented (Passot & Vázquez-Semadeni 1998; Pirozzoli & Grasso 2004; Kritsuk et al. 2007; Federrath et al. 2008; Schmidt et al. 2009; Federrath et al. 2010; Wang et al. 2011; Konstandin et al. 2012; Donzis & Jagannathan 2013; Wang et al. 2017a; Pan et al. 2019b; Rabatin & Collins 2023). However, there are still unexplored regions of the parameter space; in particular, the Re dependence of velocity and density statistics is largely unknown as most previous work in the supersonic regime are ILES.

We believe that examining the Reynolds number dependence of velocity and density statistics is important in an astrophysical context. First, in the ISM, where turbulent activity varies across subregions, the dependence of certain statistics on Reynolds numbers can reveal key information about its local structure. Second, many processes including heating and chemical reactions in the ISM, are sensitive to small-scale turbulence, particularly to kinetic energy dissipation (Godard et al. 2009). Direct Numerical Simulations (DNS) with explicitly prescribed Reynolds numbers thus offer valuable insight into small-scale dynamics and help constrain models of the ISM’s thermal and chemical energy budgets. This could not be achieved with ILES, because they lack control over the viscous cutoff scale, leading to erroneous statistics of the velocity dilatation and kinetic energy dissipation rate (Pan et al. 2019b,a). Third, the ISM exhibits internal intermittency (Falgarone & Puget 1995; Falgarone et al. 2006; Hily-Blant & Falgarone 2007; Hily-Blant et al. 2008; Falgarone et al. 2009), manifesting as extreme, localised fluctuations in gradients. For compressible, possibly supersonic turbulence, the dependence of such extreme events on Reynolds number remains an open question, which this study aims to clarify. Finally, while the ISM typically operates at Reynolds numbers exceeding those accessible in simulation, it remains unclear at what Reynolds number and at which rate, turbulence statistics converge to their asymptotic behaviour, if one were to exist (John et al. 2021). Exploring this transition in the supersonic regime helps clarify the limitations of current simulations when compared to observations.

As the ISM is a highly inhomogeneous medium, it is essential to investigate not only the spatial fluctuations of field variables but also their dependence on scale. Probability density functions (PDFs) capture spatial variability, while additional statistical tools are required to characterise the flow properties across different scales. Together, these tools offer a detailed view of ISM sub-structure. A variety of scale-dependent analytical approaches has been developed, which include analysis in Fourier space (Kritsuk et al. 2007; Federrath et al. 2010; Wang et al. 2013, 2017b; Schmidt & Grete 2019; Hellinger et al. 2021b), coarse-graining (Aluie et al. 2012; Aluie 2013; Wang et al. 2013, 2018), and point-splitting methods (correlation or structure functions) (Kritsuk et al. 2007; Federrath et al. 2009; Falkovich et al. 2010; Galtier & Banerjee 2011; Konstandin et al. 2012; Wang et al. 2013; Banerjee & Galtier 2013, 2014; Banerjee & Kritsuk 2018; Ferrand et al. 2020; Hellinger et al. 2021b; Pan et al. 2022). In astrophysical contexts, because the scale at which e.g. a star is forming can differ from the scale at which turbulence is driven (Elmegreen 2008; Federrath et al. 2017), understanding how and at which scale energy is injected, transferred, converted and dissipated is crucial. This information is again valuable for developing a more representative picture of the ISM sub-structure and ultimately more accurate models for the star formation rate (Padoan & Nordlund 2011; Hennebelle & Chabrier 2011; Federrath & Klessen 2012). The influence of Ma and Re on the scale-dependent processes is characterised in some previous work (Schmidt & Grete 2019; Hellinger et al. 2021b; Wang et al. 2018). Here we significantly extend these works by considering a larger range of Ma and Re. Of particular interest is also the convergence of scale-local mechanisms at asymptotically large Re-particularly the development of an inertial range (Aluie 2013) – mirroring trends observed in high-Re incompressible turbulence (Antonia & Burattini 2006; Danaila et al. 2012). An analogous question arises in compressible flows in the limit of infinitely large Ma.

In the present study, we perform a statistical analysis of isothermal, statistically stationary, homogeneous, and isotropic compressible turbulence, using high-fidelity numerical simulations, spanning a wide range of Ma, from subsonic to supersonic regimes, and controlled Re through DNS. We begin by examining the space-local quantities by quantifying the PDFs of the density and velocity, along with their spatial derivatives. Subsequently, we analyse the scale-dependence of the kinetic energy using a point-splitting approach based on a compressible extension of the Kármán-Howarth-Monin equation.

The remainder of the paper is organised as follows. The numerical simulations are described in Section 2. Section 3.1 presents the results concerning the PDFs of the velocity and density and their derivatives, while Section 3.2 is dedicated to the scale-by-scale energy analysis. Finally, the main conclusions are summarised in Section 4.

2. Numerical simulations

We performed numerical simulations of statistically stationary homogeneous isotropic turbulence. The continuity equation together with the Navier-Stokes equation were solved on a three-dimensional Cartesian mesh, viz

$$\begin{array}{cc}& {\partial }_t\rho +\mathbf{\nabla }·\rho \mathit{u}=0,\hfill \end{array}$$(1a)

$$\begin{array}{cc}& {\partial }_t\rho \mathit{u}+\mathbf{\nabla }·\rho \mathit{u}\mathit{u}=-\mathbf{\nabla }p+\mathbf{\nabla }·\mathsf{t}+\mathit{f},\hfill \end{array}$$(1b)

where f is a forcing term that was added in order to maintain the turbulent characteristics at statistically steady state. We used the exact same forcing method as in Federrath et al. (2010, 2021). The force f acts at large scales in a band of wavenumbers between 1 × 2π/L and 3 × 2π/L, where L = 1 is the domain width. The forcing spectrum is a parabola which peaks at a wavenumber of 2 × 2π/L, and corresponds to a scale of L/2. The forcing is composed of a half solenoidal and half compressive mode, termed natural mixture by Federrath et al. (2010). The reason for this choice is that the natural mixture is believed to be an intermediate case of all possible driving mode occurring in dense molecular clouds which ranges from purely solenoidal to purely compressive (Federrath et al. 2010; Gerrard et al. 2023). The code for generating the forcing field f is available on GitHub (Federrath et al. 2022). An isothermal equation of state for a perfect gas, p = ρc_s² was used to relate the pressure p and the density ρ through the sound speed c_s, which is constant. Hence, we did not solve an equation for the internal energy. The viscous stress tensor, t, is given by

$$\begin{array}{c}\hfill \mathsf{t}=2\mu \mathsf{S}-\frac{2}{3}\mu \theta \mathsf{I}\end{array}$$(2)

where $\mathsf{S}=(\mathbf{\nabla }\mathit{u}+\mathbf{\nabla }{\mathit{u}}^T)/2$ is the strain rate tensor, $\mathsf{I}$ is the identity matrix and θ = ∇ ⋅ u is the velocity divergence (or dilatation). Here, we made the choice of setting a constant kinematic viscosity ν = μ/ρ. Indeed, it was found that the time step constraint for the simulation was driven by the viscous term, i.e. Δt ∼ ν/Δ_x² = μ/ρΔ_x² (Δ_x is the spatial resolution). Hence, using a constant dynamic viscosity μ led to much smaller time steps compared to using a constant kinematic viscosity, especially at high Mach number where density variations are large.

For solving Eqs. (1a) and (1b), we employed an optimised version of the FLASH4 code, described in Federrath et al. (2021). It is a finite volume code which makes use of a positivity-preserving MUSCL-Hancock HLL5R Riemann scheme (Waagan et al. 2011). The simulation domain is cubic of width L = 1 with triply periodic boundary conditions. The mean density in the domain is denoted as ρ₀ = ⟨ρ⟩ (the angle brackets represent spatial and time average) and was set to 1.

Our database consists of 20 different simulations where Ma and Re vary independently (see Table 1). For this purpose, the forcing amplitude was adjusted so that the velocity dispersion u′=⟨|u|²⟩^1/2 was equal to 1 (±1%) in all cases. The Reynolds number Re = u′L/ν ≡ 1/ν and Mach number Ma = u′/c_s ≡ 1/c_s were then varied independently by changing the kinematic viscosity ν and sound speed c_s, respectively. The database covers the range 263 ≤ Re ≤ 4166 and 0.25 ≤ Ma ≤ 4.00, as reported in Table 1. For each simulation, after a transient of about 4L/u′, a statistically steady state is reached. Simulations were then run for 16 additional L/u′. All statistics presented hereafter were gathered during the steady state period. For all quantities, the statistical error was estimated as the standard deviation of the mean, divided by the square root of the number of samples. As shown later, statistical uncertainty is affecting only the far tails of the PDFs, and hence only marginally the variance of the distributions.

Care has been taken to ensure an appropriate resolution of the smallest scales of the flow. For incompressible flows, a resolution η = 2Δ_x, where η is the Kolmogorov length scale, is generally assumed to be sufficient for the velocity gradients to be accurately resolved (Yeung et al. 2018). Hence, for sub- (though slightly compressible) and transonic cases (Ma ≤ 1), we adjusted the kinematic viscosity to respect this criterion.

For the supersonic cases, a good estimate of the adequacy of the numerical resolution is the value of the pressure dilatation term ⟨pθ⟩. Indeed, as per Pan et al. (2019a), this term represents the reversible transfer between the kinetic energy and internal energy reservoirs. It should therefore be zero if the flow is isothermal and at statistically steady state. When using Riemann solvers without explicitly accounting for a viscous term (an ILES of the Euler equation), Pan et al. (2019a) observe that the statistics of the post-computed dilatation are ill-estimated. They attribute this effect to the numerical schemes which artificially enforces mass conservation at the sacrifice on some errors on θ. They also show that this bias is mitigated when considering an explicit viscous term as done by Scannapieco et al. (2024) and in the present work. Here, in order to keep the pressure dilatation term ⟨pθ⟩ of the order of a few percent of the energy injection rate, it was necessary to double the mesh resolution for Ma = 2 and 4, yielding η ≈ 4Δ_x.

Table 1 summarises some key quantities of the numerical database. The total number of points used in the simulation is demoted N. The energy injection ϵ_f = ⟨f ⋅ u⟩/ρ₀ equals the kinetic energy dissipation rate since the flow is at steady state. The Kolmogorov length scale is given by $\eta ={\nu }^{3/4}/{ϵ}_{\mathrm{f}}^{1/4}$ and is given in units of ${\mathrm{\Delta }}_x=L/\sqrt[3]{N}$. The Taylor-scale Reynolds number $R_{\lambda }=u{\prime }^2\sqrt{5/(3\nu {ϵ}_{\mathrm{f}})}$ varies between about 20 and 110. The pressure dilatation ⟨pθ⟩ is given in units of ϵ_f. The latter appears to be negative and thus represents a loss of kinetic energy. Its maximum value is about 5% of ϵ_f, although generally much smaller, which confirms that the resolution is adequate.

Table 1 indicates that when Ma increases, the amount of energy injected into the system increases in order to reach the same u′. This results in systematically smaller values of R_λ between sub- and supersonic conditions at the Re. To some extent, this is due to an increasing contribution of the dilatational component of the dissipation rate which writes ϵ_d = 4⟨μθ²⟩/3 (Jagannathan & Donzis 2016; John et al. 2021). Table 1 reveals that ϵ_d is negligible for Ma ≤ 0.5 but increases up to 25% of ϵ_f at Ma = 1 and roughly 60% (if not more) of ϵ_f for Ma ≥ 2. We note also that the ratio ϵ_d/ϵ_f increases substantially up to Ma = 2 where it seems to saturate at a value of roughly 60%. In the trans- and supersonic regimes, there is a slight effect of Re which manifests as a small increase of the ratio ϵ_d/ϵ_f when Re becomes larger. We note also that the amount of energy injected ϵ_f first decreases with Re before reaching an asymptotic value for Re ≥ 1886, i.e. R_λ ≥ 70. This reflects the evolution of the dissipation constant C_ϵ = ϵ_fL_i/u′³, where L_i is the integral length scale. Although it is not the scope of the present work, note that the finite value of ϵ_f when ν → 0 is known as the dissipative anomaly. For more details, the reader is referred to the work by John et al. (2021) who discussed the plausible existence of a dissipative anomaly for compressible flows.

3. Results

3.1. PDFs of density, velocity, and their derivatives

We start by investigating the effect of Re and Ma on the density fluctuations. The two-dimensional slices of s = lnρ/ρ₀ ≡ lnρ are displayed in Fig. 1. The first obvious observation is the increase of the density fluctuations when Ma increases. We further note a clear change of the density structure between sub- and supersonic conditions. Indeed, while the density field is smooth and looks rather similar to a passive scalar field for low Mach number, it becomes filamentary and reveals some zones of abrupt variations (fronts or shocks) for supersonic conditions. The effect of increasing Reynolds number is also discernible, which manifests by an increase of the density fluctuations at small scales. This reflects the smoothing effect due to viscosity which increases with viscosity.

Fig. 1. Two-dimensional slices of the density field s = lnρ for increasing Ma (from left to right) and increasing Re (from top to bottom). The colourbar spans ±3σs.

Quantitative insights can be provided by computing the volume weighted probability-density-functions (PDF) of s = lnρ/ρ₀ ≡ lnρ which are shown in Fig. 2a. It is very common (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Kritsuk et al. 2007; Federrath et al. 2008, 2010; Scannapieco et al. 2024) to assume that the fluctuations of ρ are log-normal. This prediction arises naturally by assuming a random multiplicative process and the application of the central limit theorem for the density evolution. Hopkins (2013) provides physical arguments showing that density fluctuations may not be log-normal due to intermittency of the field variable. A recent theoretical analysis of Scannapieco et al. (2024) suggests that density fluctuations might follow a stochastic differential process with time-correlated noise.

Fig. 2. (a) Probability-density-function of s = lnρ for varying Mach and Reynolds numbers. Colours from dark blue to yellow correspond to increasing Reynolds numbers as shown in the legend. The grey dashed lines are the associated Gaussian distributions (Eq. (3)). (b) Evolution of the variance σs2 with respect to Ma. The prediction (Eq. (4)) is also plotted for different type of forcing (compressive, solenoidal, mixed). The shaded regions represent the statistical uncertainty which affects only the far tails of the PDFs (a) and thus marginally the variance σs2 (b).

Log-normal fluctuations for ρ lead to a Gaussian distribution for s, i.e.

$$\begin{array}{c}\hfill \mathrm{PDF}(s)=\frac{1}{\sqrt{2\pi {\sigma }_{\mathrm{s}}^2}}exp(-\frac{1}{2}{\left(\frac{s-⟨s⟩}{{\sigma }_{\mathrm{s}}}\right)}^2),\end{array}$$(3)

where ⟨s⟩ and σ_s denote the mean and standard deviation of s, respectively.

Figure 2a indicates that increasing Ma yields larger density fluctuations. The peak of the distribution shifts towards the left as a consequence of mass conservation. Indeed, for a perfect log-normal distribution, mass conservation would have lead to ⟨s⟩= − σ_s/2 (Passot & Vázquez-Semadeni 1998; Kritsuk et al. 2007; Federrath et al. 2010; Rabatin & Collins 2023). In contrast, varying the Reynolds number while the Mach number is constant has a rather marginal effect on the PDF of s. This is in agreement with previous numerical results using either an explicit viscous term with different kinematic viscosity (Scannapieco et al. 2024), or an implicit numerical viscosity with different resolutions (Pan et al. 2019b; Federrath et al. 2010). Note however that although very little, there is a faster drop of the PDFs in the low-density tail when the Reynolds number increases.

Comparing the measured PDFs(s) to Gaussian distribution reveals that the curves are slightly skewed towards the low-density wing. This effect increases with the Mach number. As suggested by Federrath et al. (2010), Scannapieco et al. (2024), such skewed PDFs appear as a consequence of the compressive part of the forcing.

The evolution of the variance of s, i.e. σ_s², with respect to Ma is shown in Fig. 2b. Assuming a linear relationship between the density standard deviation σ_ρ and Ma, viz. σ_ρ = bMa, and a Gaussian distribution for s (i.e. a log-normal distribution for ρ), one arrives at (Federrath et al. 2008, 2010; Padoan & Nordlund 2011; Molina et al. 2012)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{s}}^2=ln(1+b^{2}M{a}^2)\end{array}$$(4)

where the coefficient b depends on the type of forcing. For instance, b = 1 for purely compressive forcing, b = 1/3 for purely solenoidal and b ≈ 0.4 for the mixed case used in the present work (Federrath et al. 2008, 2010). Fig. 2b indicates that the measured values of σ_s² gradually increase with Mach number. We further note that σ_s² tends to approach the prediction Eq. (4) at large Ma only (Konstandin et al. 2012; Mohapatra et al. 2021).

We now proceed similarly for the velocity field. For this purpose, we computed the volume weighted probability-density-functions of |u|/c_s which can be interpreted as the local Mach number. Recall that in our simulation, the forcing was adjusted so that the variance of |u| is constant and independent of the chosen parameters ν, c_s. Assuming that the fluctuations of each component of the velocity vector follow a Gaussian distribution with same standard deviation, one ends up with a Maxwellian distribution for |u|/c_s (Rabatin & Collins 2023). The latter is given by

$$\begin{array}{c}\hfill \mathrm{PDF}(\left|\mathit{u}\right|/c_{\mathrm{s}})=\sqrt{\frac{2}{\pi }}\frac{{(\left|\mathit{u}\right|/c_{\mathrm{s}})}^2}{M{a}_{1\mathrm{D}}^3}exp(-\frac{1}{2}\frac{{(\left|\mathit{u}\right|/c_{\mathrm{s}})}^2}{M{a}_{1\mathrm{D}}^2})\end{array}$$(5)

where $M{a}_{1D}=Ma/\sqrt{3}$. The measured PDFs of |u|/c_s are plotted in Fig. 3. As previously shown for the density field, the PDFs of |u|/c_s appear to be independent of Re. They only depend on Ma. As expected, increasing Ma yields a shift of the different distributions towards the right and their maximum decreases in proportion of 1/c_s. That simply means that the PDFs of |u| and its components u_x, u_y, u_z are independent of Ma and Re. The PDFs of |u| would have all collapsed to a single curve. The measured PDFs follow rather nicely a Maxwellian prediction, indicating that the components of u fluctuate according to a Gaussian distribution.

Fig. 3. Probability-density-function of the local Mach number |u|/cs for varying Ma and Re. The shaded regions represent the statistical uncertainty. The grey dashed lines are the associated Maxwellian distributions Eq. (5).

The fluctuations of some derivatives of ρ and u are also worth quantifying. These appear in particular in the transport equation for ρ (or s) together with the transport equation for the kinetic energy explored later in this paper. In Figs. 4 and 5, we plot the volume-weighted PDFs of ∇ρ and θ = ∇ ⋅ u, respectively. Panels from a to e are for different Ma while the colour code indicate different Re.

Fig. 4. Probability-density-function of the density gradient ∇ρ for varying Ma and Re. Colours from dark blue to yellow correspond to increasing Re as shown in Fig. 2. The shaded regions represent the statistical uncertainty. Panels (a) to (e) correspond to Ma = 0.25 to Ma = 4.

Fig. 5. Same as Fig. 4 for the velocity divergence θ = ∇ ⋅ u. The dashed lines in panel (e) are the reproduction of the data for Ma = 2 (panel (d)).

Contrary to s and |u|, the fluctuations of the density gradient ∇ρ and dilatation θ depend on both Ma and Re. Increasing the Mach and Reynolds numbers yields larger tails in the PDFs for both ∇ρ and θ. This results in sharper density/velocity fronts when compressibility effects (i.e. the Mach number) become larger or viscous effects (i.e. the Reynolds number) get weaker. The measured PDFs deviate very significantly from Gaussian distributions revealing that ∇ρ and θ are highly intermittent. For instance, at Ma = 4 and Re = 4166, some local fluctuations of ∇ρ can exceed 100 times their standard deviation. Note that in isothermal flows, the pressure p = ρc_s² and hence the PDFs portrayed in Fig. 4 would be the exact same for ∇p.

The PDFs of dilatation θ (Fig. 5) reveal some skewed distributions that are typical of compressible flows (Passot & Vázquez-Semadeni 1998; Jagannathan & Donzis 2016; Schmidt et al. 2009; Wang et al. 2017a; Sakurai & Ishihara 2023, 2024). At large Ma and Re, the left wing of the PDFs, corresponding to compression regions where θ < 0, have a much larger extent than the right wing which drops very quickly. Strong compression events are thus more likely than strong expansion events. Despite, the total volume where θ < 0 is compensated by the volume formed by θ > 0 so that ⟨θ⟩ = 0. At the lowest Ma, we note that the PDFs(θ) are rather symmetric with respect to the θ = 0 axis and are almost Gaussian. Pirozzoli & Grasso (2004) and Sakurai et al. (2021), report similar observation but at a lower Mach number. Comparing the distributions at Ma = 1,  2 and 4, reveals that the PDFs get narrower when Ma increases from 1 to 2 (Figs. 5c and d). Wang et al. (2017a) made similar observations comparing the PDFs(θ) at Ma = 0.8 and Ma = 1 (see their Fig. 2). Note that the x-coordinate axis in Fig. 5 is normalised by σ_θ which increases between Ma = 1 and Ma = 2 (this quantity will be discussed in some forthcoming paragraphs).

Our data further indicate that the PDFs(θ) are almost unchanged between Ma = 2 and Ma = 4 (Figs. 5d and e). This is visible in Fig. 5e where the dashed lines corresponding to Ma = 2 collapse with the full lines corresponding for Ma = 4. Speculatively, one could expect that increasing further the Mach number does not have any effects on the PDFs(θ). This suggests that for Ma ≥ 2 the distribution of θ has reached a saturated state, which depends only on the Reynolds number (and probably the type of forcing) but not on the Mach number.

The variance of ∇ρ and θ, are denoted σ_∇ρ² and σ_θ², respectively. Their evolution with respect to Re and Ma are reported in Figs. 6a–b. Noticeable is the very strong increase of σ_∇ρ² when Ma varies from 0.25 to 4 (Fig. 6a). Our data suggests a power-law variation of the form σ_∇ρ² ∼ Ma⁴. It also increases with the Reynolds number although with a lower exponent, viz. σ_∇ρ² ∼ Re¹. The variance of the dilatation is presented in Fig. 6b. We first note a prompt increase of σ_θ for subsonic (though compressible) conditions before reaching a plateau for Ma ≥ 2, whose level depends only on the Reynolds number. Here again, we find that the values of σ_θ² at a given Mach number increase proportionally to Re. For comparison, we have plotted the results reported by Scannapieco et al. (2024). We selected only the cases labelled M1.3ν₀3, M2ν₀3, M3ν₀3 in Scannapieco et al. (2024) as they correspond to roughly the same effective Reynolds number (Re ∼ 3400 according to our definition). Although quantitatively different, their data are in agreement with our finding that σ_θ²/(u′/L)² may saturate for Ma ≥ 2. Since dilatation appears in the vicinity of shocks (Wang et al. 2011), the saturation of σ_θ² indicates that the shock strength at large Ma is limited by viscous effects. This saturation effect is also observed for the density-weighted variance of θ (open symbols in Fig. 6b). Note that since we have used a fluid with a constant kinematic viscosity ν, the density-weighted variance of θ, i.e. ⟨ρθ²⟩/ρ₀ = 3ϵ_d/(4ν). Therefore, the saturation of ⟨ρθ²⟩ for Ma ≥ 2 means that the dilatational component of the kinetic energy dissipation ϵ_d also saturates to a certain value that uniquely depends on Re.

Fig. 6. Mach-number dependence of the variance of (a) the density gradient ∇ρ (in units of ρo/L ≡ 1) and (b) dilatation θ (in units of u′/L ≡ 1) for varying Re. In panel (b), the triangle symbols represent the data reported by Ref. Scannapieco et al. (2024) when rescaled in units of (u′/L)2. The closed (open) symbol represent the volume (density) weighted variance.

3.2. Scale-by-scale kinetic energy budgets

Light is now shed onto the Mach and Reynolds numbers effects on the kinetic energy distribution and kinetic energy transfers. For this purpose, we use a generalised version of the Kármán-Howarth-Monin (KHM) equation which accounts for the variations of density within the flow. The latter starts from the definition of the turbulent kinetic energy at a given scale proposed by e.g. Galtier & Banerjee (2011), Lai et al. (2018), Hellinger et al. (2021b), viz.

$$\begin{array}{c}\hfill ⟨|\delta \mathit{u}{|}_{\rho }^2⟩=⟨(\delta (\rho \mathit{u}))·(\delta \mathit{u})⟩\end{array}$$(6)

where δ(ρu) = (ρu)⁺ − (ρu)⁻ and δu = u⁺ − u⁻ represents the increment (the difference) of ρu and u, respectively, between two points x⁺ and x⁻ arbitrarily separated in space by a distance r, i.e. r = x⁺ − x⁻. The superscript +/− indicates that the field variable is taken at point x⁺/x⁻. The brackets in Eq. (6) denote a suitable average that depends on the symmetry of the flow. Here, the flow being statistically stationary, homogeneous, and isotropic, the average is performed over statistically independent samples, over space and over all possible orientations of the vector r using an angular average. By doing this, the quantity ⟨|δu|_ρ²⟩ depends only on the modulus of the separation vector r = |r|. It is generally interpreted as the turbulent kinetic energy at a given scale r. Note that other definitions for the scale-by-scale turbulent kinetic energy in variable-density flows can be found in the literature (see for instance Ferrand et al. 2020; Brahami 2020; Hellinger et al. 2021a). It is worth noting that for homogeneous flows, the limit of Eq. (6) at large separations yields

$$\begin{array}{c}\hfill \underset{r\to \infty }{lim}{⟨|\delta \mathit{u}|}_{\rho }^2{⟩=2⟨|\mathit{u}|}_{\rho }^2{⟩=2⟨\rho |\mathit{u}|}^2⟩,\end{array}$$(7)

where ⟨ρ|u|²⟩ is twice the turbulent kinetic energy. Hence, ⟨|δu|_ρ²⟩ is not an energy density (such as provided by spectra) but should rather be interpreted as a sort of cumulative kinetic energy distribution of all scales ≤r. The transport equation for ⟨|δu|_ρ²⟩ is known as the Kármán-Howarth-Monin equation and can formally be written as (Galtier & Banerjee 2011; Hellinger et al. 2021b):

$$\begin{array}{c}\hfill d_t{⟨|\delta \mathit{u}|}_{\rho }^2⟩=⟨\mathcal{T}⟩+⟨\mathcal{P}⟩+⟨\mathcal{V}⟩+⟨\mathcal{F}⟩\end{array}$$(8)

The different terms of Eq. (8) are explicited below:

• d_t⟨|δu|_ρ²⟩ denotes the time variations of ⟨|δu|_ρ²⟩. In statistically stationary flows, this term is zero.

• ⟨𝒯⟨ represents the “transport” of ⟨|δu|_ρ²⟩ and stems from the non-linear transport term of the Navier-Stokes equation. In statistically homogeneous flows, the latter can be decomposed into two contributions:

  $$\begin{array}{c}\hfill ⟨\mathcal{T}⟩=-⟨{\mathbf{\nabla }}_{\mathit{r}}·(\delta \mathit{u})|\delta \mathit{u}{|}_{\rho }^2⟩+⟨\mathcal{R}⟩.\end{array}$$(9)

  The leftmost term on RHS of Eq. (9) represents the transfer (or cascade) of kinetic energy among the different scales of the flow. It writes as the divergence in scale-space ∇_r⋅ of the flux (δu)|δu|_ρ². The term ⟨ℛ⟩ is present only in compressible flows. It represents a source term due to the dilatation and writes

  $$\begin{array}{c}\hfill ⟨\mathcal{R}⟩=\langle (\delta \mathit{u})·[(\delta (\rho \mathit{u}))(\overline{\delta }\theta )-(\overline{\delta }(\rho \mathit{u}))(\delta \theta )]\rangle ,\end{array}$$(10)

  where $\overline{\delta }\bullet =({\bullet }^{+}+{\bullet }^{-})/2$ denotes the arithmetic mean of any quantity •, between x⁺ and x⁻.

• ⟨𝒫⟨ represents the effect of pressure on the evolution of the scale-by-scale kinetic energy. In homogeneous compressible flows, it can be decomposed into two contributions:

  $$\begin{array}{c}\hfill ⟨\mathcal{P}⟩=2⟨(\delta P)(\delta \theta )⟩-⟨\mathcal{C}(-\mathbf{\nabla }p)⟩,\end{array}$$(11)

  The first term on RHS of Eq. (11) is the scale-by-scale contribution of the pressure-dilatation correlation. It corresponds to the scale-by-scale conversion between kinetic energy and internal energy (Aluie et al. 2012). The term −⟨𝒞(−∇P)⟩ arises due to variations of the density. The latter was derived by Hellinger et al. (2021b) and writes:

  $$\begin{array}{cc}\hfill \mathcal{C}(\mathit{a})& :=(\delta \mathit{u})·(\delta \mathit{a})-(\delta (\rho \mathit{u}))·(\delta v\mathit{a})\hfill \\ \hfill & =({\rho }^{+}{v}^{-}-1){\mathit{u}}^{+}·{\mathit{a}}^{-}+({\rho }^{-}{v}^{+}-1){\mathit{u}}^{-}·{\mathit{a}}^{+},\hfill \end{array}$$(12)

  where v = 1/ρ is the specific volume and a denotes any vectorial field quantity. For the pressure term, we set a ≡ −∇p. It is easy to show that in constant-density flows, ρ^±v^∓ = 1 and hence 𝒞(a) = 0. This term reveals a triple correlation between density, velocity, and pressure gradient. It is therefore rather similar to what Aluie et al. (2012) refer to as the baropycnal work, and thus expected to contribute as an additional process of kinetic energy transfer across scales.

• The terms ⟨𝒱⟩ and ⟨ℱ⟩ represent the scale-by-scale contribution of viscous diffusion and forcing, respectively. These terms are here written in the same formulation as Hellinger et al. (2021b), viz.

  $$\begin{array}{cc}\hfill ⟨\mathcal{V}⟩& =2⟨(\delta \mathit{u})·(\delta (\mathbf{\nabla }·\mathsf{t}))⟩-⟨\mathcal{C}(\mathbf{\nabla }·\mathsf{t})⟩\hfill \end{array}$$(13a)

  $$\begin{array}{cc}\hfill ⟨\mathcal{F}⟩& =2⟨(\delta \mathit{u})·(\delta \mathit{f})⟩-⟨\mathcal{C}(\mathit{f})⟩\hfill \end{array}$$(13b)

  where t is the viscous stress tensor defined in Eq. (2) and f is the forcing term described in Section 2. The expanded version of the viscous term has been derived by Lai et al. (2018). We keep it here in the compact form given by Eq. (13a).

It can be proven that, for statistically homogeneous flows, the limit at large separations of each term of Eq. (8) is

$$\begin{array}{cc}\hfill \underset{r\to \infty }{lim}⟨\mathcal{T}⟩& =0,\hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill \underset{r\to \infty }{lim}⟨\mathcal{P}⟩& =4⟨p\theta ⟩,\hfill \end{array}$$(14b)

$$\begin{array}{cc}\hfill \underset{r\to \infty }{lim}⟨\mathcal{V}⟩& =-4ϵ,\hfill \end{array}$$(14c)

$$\begin{array}{cc}\hfill \underset{r\to \infty }{lim}⟨\mathcal{F}⟩& =4{ϵ}_{\mathrm{f}},\hfill \end{array}$$(14d)

where, $ϵ=⟨\mathbf{\nabla }\mathit{u}:\mathsf{t}⟩$ (the colon operator denotes the double contraction of second-order tensors) is the kinetic energy dissipation rate. Hence, the one-point kinetic energy budget

$$\begin{array}{c}\hfill d_t\frac{1}{2}{⟨|\mathit{u}|}_{\rho }^2⟩=⟨p\theta ⟩-ϵ+{ϵ}_{\mathrm{f}},\end{array}$$(15)

is recovered from Eq. (8) in the limit of large separations (ρ₀ = 1 is dropped from the notations). Recall that in statistically stationary flows, both the time derivative and pressure-dilatation terms vanish (Pan et al. 2019a) and one ends up with ϵ = ϵ_f. The kinetic energy dissipation rate ϵ = ϵ_s + ϵ_d, where, for homogeneous flows, ϵ_s = ⟨μ|ω|²⟩ (ω is the vorticity vector) and ϵ_d = 4⟨μθ²⟩/3 are the solenoidal and dilatational components of the dissipation rate, respectively.

The effect of Re and Ma on the scale-by-scale kinetic energy distribution ⟨|δu|_ρ²⟩ is presented in Figs. 7a and b, respectively. The separation r is normalised using L while ⟨|δu|_ρ²⟩ is divided by its large-scale asymptotic value, viz. 2⟨|u|_ρ²⟩. In Figs. 7a–b, the different sets of curves corresponding to either constant-Ma or constant-Re are shifted upwards for clarity.

Fig. 7. Distribution of kinetic energy ⟨|δu|ρ2⟩/2⟨|u|ρ2⟩ as a function of r/L. Panel (a) shows the effect of Re at constant Ma, while figure (b) does the opposite. The inset in (b) represents the local scaling exponent, Eq. (16), at Re = 4167. For clarity, each group of curves in (a) and (b) are shifted upwards by a factor 16 and 4, respectively.

We note that, at constant Ma (Fig. 7a), the different curves expand in the direction of smaller scales when Re increases. This results from the decreasing value of the viscous cut-off scale. A careful examination of Fig. 7b reveals that there is no noticeable effect of Ma on ⟨|δu|_ρ²⟩ at the lowest Reynolds numbers. Its effect becomes substantial only at the largest Reynolds number and in the intermediate range of scales, where one can expect a power-law behaviour of ⟨|δu|_ρ²⟩ with respect to the scale r, i.e. ⟨|δu|_ρ²⟩∼r^ζ. The curves corresponding to the largest Reynolds number in Fig. 7b indicate that the scaling exponent is larger for supersonic than subsonic conditions. This can be better quantified by looking at the local scaling exponent defined by:

$$\begin{array}{c}\hfill \zeta (r)=\frac{dln⟨|\delta \mathit{u}{|}_{\rho }^2⟩}{dlnr}\end{array}$$(16)

The latter is portrayed in the inset of Fig. 7b for Re = 4167. Note that the largest Reynolds number studied here is likely too moderate for a clear scaling range to be unambiguously perceived. We observe though that up to Ma = 1, the scaling exponent in the intermediate range of scales deflects to a value not so far from the Kolmogorov’s prediction ζ = 2/3 (Kolmogorov 1941). For Ma = 2 and Ma = 4, ζ deflects to a larger value that seems to comply with the prediction for Burgers turbulence ζ = 1 (Saffman 1968; Alam et al. 2022). An analogous observation was reported by Schmidt & Grete (2019) using a spectral decomposition of the velocity. Note that towards the smallest scales, the local scaling exponent ζ(r)→2 for all values of Ma. This agrees with the asymptotic scaling as r → 0 for both Kolmogorov and Burgers turbulence (Kolmogorov 1941; Alam et al. 2022).

In summary, in sub- up to transonic conditions, compressibility effects are relatively weak and the scale-by-scale energy distribution is not so different from what one expects for incompressible (Kolmogorov) turbulence. In the supersonic regime, compressibility effects appear to affect mainly the intermediate range of scales, where the scaling exponent is likely to comply with a Burgers kind of turbulence. More studies at higher Re and Ma are likely needed in order to draw firm conclusions about the scaling exponents in the intermediate range of scales. Note that Federrath et al. (2021) revealed that at sufficiently large Reynolds and Mach numbers, both Kolmogorov and Burgers scalings may coexist in two separate range of scales delimited by the so-called sonic scale, i.e. the scale at which the local-in-scale Mach number is 1. From the present data, this double scaling does not appear as a consequence of the moderate Reynolds number.

The scale distribution of the different terms of the KHM equation Eq. (8) are presented in Fig. 8. All terms are divided by ϵ_f. By doing this, the contribution of the forcing term ⟨ℱ⟩ appears to be independent of the conditions investigated. All other terms can then be studied in proportion of an invariant injection of kinetic energy which significantly eases the interpretation of the results. In panel a, we have combined the results for Ma = 0.25 (presented with a slight transparency) and Ma = 0.5. Panels b–d are for Ma = 1, Ma = 2, Ma = 4, respectively.

Fig. 8. Scale-by-scale kinetic energy budgets Eq. (8) from sub- to supersonic conditions. The different terms of the KHM equation are normalised by 4ϵf while the separation is divided by L ≡ 1. The colour code is same as in Fig. 7a. The shaded regions represent the statistical error. Panel (a) is for Ma = 0.25 and Ma = 0.5, the symbols are the results obtained from an incompressible code. Panel (b), (c) and (d) are for Ma = 1, Ma = 2 and Ma = 4, respectively.

Before we start analysing the KHM equation, it is worth recalling that ⟨|δu|_ρ²⟩ represents a cumulative energy distribution rather than an energy density. The same applies to the different terms of its transport equation (Eq. (8)). This is important for interpreting for instance the scale contribution of the transport term ⟨𝒯⟩ and pressure term ⟨𝒫⟩ which are presented in Figs. 8a–d. The latter have a bell shape. They increase when travelling from small to large scales, reach a maximum value and then decrease again to zero at large scales. Since they are cumulative distributions, it means that scales smaller (larger) than the scale at which ⟨𝒯⟩ or ⟨𝒫⟩ is maximum receives (looses) kinetic energy. Hence, such terms are characteristics of the energy cascade or energy conversion between the different scales of the flow. The overall scale-integrated energy transfer/conversion due to the non-linear and pressure terms of the Navier-Stokes equation is zero and hence ⟨𝒯⟩→0 and ⟨𝒫⟩→0 in the limit of very large separations. The forcing and viscous term do not reveal such a trend but rather asymptote to a non-zero value at large scales which corresponds to the one-point kinetic energy budget given here by ϵ_f = ϵ. The viscous term ⟨𝒱⟩ in Fig. 8 is classically presented in the form ⟨𝒱⟩+4ϵ_f, and hence ⟨𝒱⟩+4ϵ_f → 4ϵ_f when r → 0 and ⟨𝒱⟩+4ϵ_f → 0 when r → ∞. By doing this, the viscous term appears to be confined to small scales in agreement with the theoretical reasoning of Aluie (2013). Irrespectively of Re and Ma, the terms of the KHM equation presented in Figs. 8a-d indicate that kinetic energy injection is confined to large scales (Aluie et al. 2012; Aluie 2013), then transferred/converted to smaller scales by the non-linear transport and pressure terms (when the latter contributes), and is finally dissipated by the viscous term at small scales.

Looking first to the subsonic cases Ma = 0.25 − 0.5 (Fig. 8a) reveals that the pressure term ⟨𝒫⟩ is negligible. In other words, the cases of Ma = 0.25 − 0.5 do not depart significantly from incompressible flows where only the forcing, transport and viscous term contribute to the budget. This is confirmed in Fig. 8a where the results for Ma = 0.25 − 0.5 are compared to their equivalent obtained from a fully incompressible turbulence code (Thiesset & Vahé 2025) using the exact same fluid/forcing parameters. However, it is possible that higher values of Re are needed to perceive compressibility effects at such Mach numbers. The analysis of Wang et al. (2018) at R_λ ≈ 250 indicates that the different terms of the ensemble average coarse-grained kinetic budget are not very sensitive to Ma in the subsonic regime. Figure 8a indicates that, in the intermediate range of scales, the transport term increases when Re increases. Meanwhile, the viscous term moves towards smaller scales. This behaviour suggests that in the limit of very large Re, only the transport term contributes to the budget in the intermediate range, which yields the celebrated Kolmogorov 4/5-law (Nie & Tanveer 1999; Antonia & Burattini 2006; Falkovich et al. 2010; Danaila et al. 2012), which here can be written as ⟨𝒯⟩ = 4ϵ_f for η ≪ r ≪ L.

For larger Ma (panels b–d), the results are rather surprising. For instance, we observe that the contribution of the transport term ⟨𝒯⟩ has a non-monotonic behaviour. Compared to Ma = 0.25 and Ma = 0.5, it decreases in amplitude up to Ma = 2 (Fig. 8c) before increasing again for Ma = 4 (Fig. 8d) to even overtake the incompressible case. A similar, though opposite non-monotonic behaviour applies to the pressure term. The latter increases in amplitude up to Ma = 2 before decreasing at Ma = 4. It is likely that increasing further the Mach number results in an even smaller contribution of the pressure term so that ⟨𝒫⟩ may become negligible at very large Ma. We also note that ⟨𝒫⟩ has the same shape and same sign as the transport term, although it peaks at smaller scales. Scrutinising the curves for Ma = 1 to Ma = 4 (Figs. 8b–d) indicates that the viscous term tends to contribute at smaller and smaller scales when Re and Ma increases. As a consequence, the total transfer constituted of the non-linear and pressure terms approach ϵ_f in the limit of large Reynolds numbers. We may thus extend the Kolmogorov 4/5-law to compressible flows by writing ⟨𝒯⟩+⟨𝒫⟩ = 4ϵ_f for η ≪ r ≪ L, which can be recast in the same form as Eq. (3.3) of Falkovich et al. (2010). As already stated, it is plausible that the term ⟨𝒫⟩ vanishes in the limit of large Ma and the generalised Kolmogorov law may reduce to ⟨𝒯⟩ = 4ϵ_f for η ≪ r ≪ L when c_s → 0.

We also quantified the scale-by-scale distributions of the terms where dilatation θ contributes, i.e. the pressure-dilatation term 2⟨δθδp⟩ and the source term ⟨ℛ⟩. Results are presented in Fig. 9. As observed before, these terms are negligible for Ma = 0.25 and Ma = 0.5 and start to contribute significantly only for Ma ≥ 1. They have the same bell shape as the pressure and non-linear transport term presented before but are mainly negative. This means that dilatation acts as a gain (loss) of kinetic energy for the large (small) scales. Overall, dilatation thus opposes the energy transfer of kinetic energy among the different scales. We further note that the pressure-dilatation term peaks at smaller scales than the source term ⟨ℛ⟩. The latter term ⟨ℛ⟩ appears to increase in amplitude and reaches values close to 4ϵ_f at the largest Ma − Re condition investigated. We do not observe any saturation of this term when Ma and Re increase. In contrast, the two-point pressure-dilatation correlation increases in amplitude up to Ma = 2 and finally saturates to a certain distribution for larger Ma. The saturation of the one-point statistics of dilatation was already identified before in this work. It appears that a similar behaviour applies also for two-point statistics and in particular for the scale-by-scale contribution of the pressure-dilatation term ⟨δθδp⟩. Again, this indicates that the dilatation is limited by viscous effects at high Ma.

Fig. 9. Scale-by-scale contribution of the terms associated to the dilatation, i.e. ⟨ℛ⟩ and 2⟨δθδp⟩. The colour code is same as in Fig. 7a.

Overall, our data thus suggest that the kinetic energy exchanges for compressible turbulence operates according to the following scenario. As for incompressible turbulence, part of the energy injected at large scales is transferred to smaller scales through the non-linear transport term. The specificity of compressible flows is that another part of the injected kinetic energy is transferred to smaller scales by the baropycnal work. The pressure-dilatation acts in converting the small-scale kinetic energy into the internal energy reservoir and the opposite conversion occurs at large scales.

4. Summary

We analyse data from numerical simulations of compressible isothermal turbulence covering different regimes from sub- to trans- to supersonic conditions at different Reynolds numbers. This database extends previous work on the same topic by addressing the case of larger Mach number situations, albeit with some more moderate Reynolds numbers. Care is taken to the appropriate numerical resolution of the equations by scrutinising the residual of the pressure-dilation term in the one-point kinetic energy budget.

Light is shed on some statistics of the density and velocity fields together with some of their spatial derivatives. We confirm that the statistics of density and local Mach number depend mainly on Ma but not on Re. Their spatial derivatives however show a dependence on both Re and Ma which is associated with an increased level of intermittency. Increased intermittency has significant consequences for the structure and kinematics, as well as the thermodynamics (as the most intermittent structures are believed to be the sites of the most intense dissipation Falgarone et al. (2006), Hily-Blant & Falgarone (2007), Hily-Blant et al. (2008)) of the interstellar medium. Indeed, it indicates very local fluctuations of density and velocity and therefore of the kinetic energy dissipation rate with implications for the thermal and chemical evolution of the interstellar medium.

The results at larger Ma allow us to identify that the fluctuations of the dilatation saturate to a certain level that vary with Re but not Ma. This indicates that the shock intensity cannot grow indefinitely as it is limited by viscosity. The consequence is that the dilatational component of the kinetic energy dissipation rate ϵ_d also saturates to a certain proportion of the total kinetic energy dissipation rate ϵ. Since dilatation plays a crucial role in the dynamic evolution of the density field (Scannapieco et al. 2024), the saturation of its fluctuations may have significant consequences for the density statistics and therefore the star formation rate in astrophysical applications.

We also explore the kinetic energy exchanges among the different scales of the flow. For this purpose we use a generalised version of the KHM-equation, i.e. the transport equation for the scale-by-scale kinetic energy, which accounts for the variations of density and non-zero dilatation. We find that, at constant Reynolds numbers, the effect of Mach number on the kinetic energy distribution is perceived mainly in the intermediate range of scales. We show that, for supersonic conditions, the scaling exponent in the inertial range approaches Burgers’ predictions, while Kolmogorov scaling applies to sub- and transonic conditions.

The terms of the KHM-equation up to Ma = 0.5 are identical to their counterpart in incompressible flows. It is not excluded however that, for such values of Ma, compressibility effects become perceptible at larger Re. For transonic and supersonic conditions, the effect of Ma is rather surprising. Indeed, it is reported that the different terms of the KHM equation have a non-monotonic behaviour. For instance, the maximum transfer of kinetic energy due to the non-linear (pressure) terms of the Navier-Stokes equation appears to decrease (increase) up to Ma = 2 before increasing (decreasing) again for larger values. The two-point pressure-dilatation correlation saturates to a certain distribution above a certain Ma threshold. Even though this needs to be confirmed by some data at larger Re and Ma, our analysis suggests that the generalised 4/5-law proposed by Falkovich et al. (2010), i.e. ⟨𝒯⟩+⟨𝒫⟩=ϵ_f may hold in the limit of very large Reynolds numbers. For infinitely large Mach numbers (c_s → 0), our data also indicates that the pressure term is likely to vanish and the 4/5-law may reduce to ⟨𝒯⟩=ϵ_f. In supersonic situations, the non-linear transport term ⟨𝒯⟩ has a strong contribution due to dilatation, i.e. ⟨ℛ⟩, which represents a loss (gain) of kinetic energy for the small (large)scales.

Observational data of the ISM often provides some measures of the velocity structure functions. The KHM is the theoretical framework to interpret these measurements in terms of the underlying physical processes, such as energy injection, transfer, conversion, and dissipation. Our results indicate that Re and Ma significantly influences these processes. We believe that this insight can be useful for better interpreting observational data and refining models of ISM turbulence.

The present study opens several avenues for future research. As highlighted by Donzis & John (2020), the Mach and Reynolds numbers alone do not fully characterise the statistical behaviour of compressible homogeneous turbulence. A more comprehensive description requires the inclusion of an additional parameter: the ratio of solenoidal to dilatational velocity fluctuations. In this work, we restrict our analysis to the so-called natural mixture of forcing, which consists of an equal proportion of solenoidal and compressive modes. It would therefore be of interest to systematically vary this proportion, thereby modifying the solenoidal-to-dilatational ratio, to assess whether our conclusions will hold qualitatively under different forcing conditions. More generally, further insight could be gained by examining scale-local processes after decomposing the velocity field into its solenoidal and dilatational components. Secondly, we hypothesise that, at sufficiently high Ma, dilatational fluctuations are constrained by viscous effects. This conjecture could be tested by quantifying both the width and intensity of shock structures. Thirdly, simulations at higher Re and Ma will be required to confirm the emergence of the aforementioned asymptotic relations, analogous to the classical 4/5-law.

Data availability

The dataset supporting this study, including both raw and plotting data, comprises approximately 55 TB. These data are available from the authors upon request, subject to availability. Please note that certain raw datasets may no longer be accessible due to storage capacity limitations.

Acknowledgments

We thank the anonymous referee for their constructive feedback, which has helped us improve the manuscript. Calculations were performed using the computing resources of CRIANN (Normandy, France), under the project 2018002. We also benefited from resources from the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. C. F. acknowledges funding provided by the Australian Research Council (Discovery Project DP230102280 and DP250101526), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). 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) and the Pawsey Supercomputing Centre (project pawsey0810) 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 Flash Centre for Computational Science at the University of Chicago and the Department of Physics and Astronomy of the University of Rochester.

References

All Tables

All Figures

	Fig. 1. Two-dimensional slices of the density field s = lnρ for increasing Ma (from left to right) and increasing Re (from top to bottom). The colourbar spans ±3σs.
In the text

Fig. 2. (a) Probability-density-function of s = lnρ for varying Mach and Reynolds numbers. Colours from dark blue to yellow correspond to increasing Reynolds numbers as shown in the legend. The grey dashed lines are the associated Gaussian distributions (Eq. (3)). (b) Evolution of the variance σs2 with respect to Ma. The prediction (Eq. (4)) is also plotted for different type of forcing (compressive, solenoidal, mixed). The shaded regions represent the statistical uncertainty which affects only the far tails of the PDFs (a) and thus marginally the variance σs2 (b).

In the text

	Fig. 3. Probability-density-function of the local Mach number |u|/cs for varying Ma and Re. The shaded regions represent the statistical uncertainty. The grey dashed lines are the associated Maxwellian distributions Eq. (5).
In the text

	Fig. 4. Probability-density-function of the density gradient ∇ρ for varying Ma and Re. Colours from dark blue to yellow correspond to increasing Re as shown in Fig. 2. The shaded regions represent the statistical uncertainty. Panels (a) to (e) correspond to Ma = 0.25 to Ma = 4.
In the text

	Fig. 5. Same as Fig. 4 for the velocity divergence θ = ∇ ⋅ u. The dashed lines in panel (e) are the reproduction of the data for Ma = 2 (panel (d)).
In the text

	Fig. 6. Mach-number dependence of the variance of (a) the density gradient ∇ρ (in units of ρo/L ≡ 1) and (b) dilatation θ (in units of u′/L ≡ 1) for varying Re. In panel (b), the triangle symbols represent the data reported by Ref. Scannapieco et al. (2024) when rescaled in units of (u′/L)2. The closed (open) symbol represent the volume (density) weighted variance.
In the text

	Fig. 7. Distribution of kinetic energy ⟨|δu|ρ2⟩/2⟨|u|ρ2⟩ as a function of r/L. Panel (a) shows the effect of Re at constant Ma, while figure (b) does the opposite. The inset in (b) represents the local scaling exponent, Eq. (16), at Re = 4167. For clarity, each group of curves in (a) and (b) are shifted upwards by a factor 16 and 4, respectively.
In the text

Fig. 8. Scale-by-scale kinetic energy budgets Eq. (8) from sub- to supersonic conditions. The different terms of the KHM equation are normalised by 4ϵf while the separation is divided by L ≡ 1. The colour code is same as in Fig. 7a. The shaded regions represent the statistical error. Panel (a) is for Ma = 0.25 and Ma = 0.5, the symbols are the results obtained from an incompressible code. Panel (b), (c) and (d) are for Ma = 1, Ma = 2 and Ma = 4, respectively.

In the text

	Fig. 9. Scale-by-scale contribution of the terms associated to the dilatation, i.e. ⟨ℛ⟩ and 2⟨δθδp⟩. The colour code is same as in Fig. 7a.
In the text