# Statistical properties of compressible isothermal turbulence from sub- to supersonic conditions

F. Thiesset, C. Federrath

► To cite this version:

F. Thiesset, C. Federrath. Statistical properties of compressible isothermal turbulence from sub-to supersonic conditions. Astronomy & Astrophysics - A&A, 2025, 703, pp.A209. 10.1051/0004-6361/202555653. hal-05386981

HAL Id: hal-05386981

https://hal.science/hal-05386981v1

Submitted on 28 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistical properties of compressible isothermal turbulence from sub- to supersonic conditions

F. Thiesset<fnref n="1" />,*, and C. Federrath<fnref n="2" />,<fnref n="3" />,★

Received 24 May 2025 / Accepted 14 October 2025

## ABSTRACT

**Context.** Turbulence is one of the key processes that control the spatial and temporal evolution of matter and energy of many astro-physical systems.

**Aims.** This paper investigates the statistical properties of isothermal turbulence in both the subsonic and supersonic regimes. The focus is on the influence of the Mach number (*Ma*) and the Reynolds number (*Re*) on both the space-local and scale-dependent fluctuations of relevant gas variables, the density, velocity, their derivatives, and the kinetic energy.

**Methods.** We carried out hydrodynamical simulations of driven turbulence with explicit viscosity and therefore controlled *Re*, at converged numerical resolutions up to 1920³ grid cells.

**Results.** We confirm previous work that the probability density functions (PDFs) of the gas density are approximately log-normal and depend on *Ma*. We provide a new detailed quantification of the dependence of the PDFs of density and velocity on *Re*, finding a relatively weak dependence, provided *Re* > 200. In contrast, derivatives of the density and velocity field are sensitive to *Re*, with the probability of extreme events (the tails of the PDFs) growing with *Re*. The PDFs of the density gradient and velocity divergence (dilatation) exhibit increasingly heavy tails with growing *Ma*, signalling enhanced internal intermittency. At sufficiently high *Ma*, the statistics of dilatation are observed to saturate at a level determined solely by *Re*, suggesting that turbulent dilatation becomes limited by viscous effects. We also examine the scale-by-scale distribution of kinetic energy through a compressible form of the Kármán-Howarth-Monin (KHM) equation that incorporates density variations. In the intermediate range of scales, a marked difference is found between subsonic and supersonic turbulence: while Kolmogorov-like scaling applies in the sub- and transonic regimes, supersonic turbulence aligns more closely with Burgers turbulence predictions. The analysis of individual terms in the KHM equation highlights the role of the pressure-velocity coupling as an additional mechanism for converting and transferring kinetic energy from large to small scales. Moreover, the contributions of the KHM terms exhibit non-monotonic behaviour with increasing *Ma*, with dilatational effects becoming more pronounced and acting to oppose the cascade of kinetic energy.

**Key words.** hydrodynamics – turbulence

## 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

<footnote n="1">¹ CNRS, CORIA, UMR 6614, Normandy Univ., UNIROUEN, INSA Rouen, 675 Avenue de l'université, BP 12, 76801 Saint Etienne du Rouvray Cedex, France</footnote>

<footnote n="2">² Research School of Astronomy and Astrophysics, Australian National University, Cotter Road, Canberra, ACT 2611, Australia</footnote>

<footnote n="3">³ Australian Research Council Centre of Excellence in All Sky Astrophysics (ASTRO3D), Cotter Road, Canberra, ACT 2611, Australia</footnote>

<footnote>* Corresponding authors: fabien.thiesset@cnrs.fr; christoph.federrath@anu.edu.au</footnote>
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

$$ \partial_t \rho + \nabla \cdot \rho \mathbf{u} = 0, \quad (1a) $$

$$ \partial_t \rho \mathbf{u} + \nabla \cdot \rho \mathbf{u} \mathbf{u} = -\nabla p + \nabla \cdot \mathbf{t} + \mathbf{f}, \quad (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 \times 2\pi/L$ and $3 \times 2\pi/L$, where $L = 1$ is the domain width. The forcing spectrum is a parabola which peaks at a wavenumber of $2 \times 2\pi/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 = \rho c_s^2$ was used to relate the pressure $p$ and the density $\rho$ through the sound speed $c_s$, which is constant. Hence, we did not solve an equation for the internal energy. The viscous stress tensor, $\mathbf{t}$, is given by

$$ \mathbf{t} = 2\mu\mathbf{S} - \frac{2}{3}\mu\theta\mathbf{l} \quad (2) $$

where $\mathbf{S} = (\nabla \mathbf{u} + \nabla \mathbf{u}^T)/2$ is the strain rate tensor, $\mathbf{l}$ is the identity matrix and $\theta = \nabla \cdot \mathbf{u}$ is the velocity divergence (or dilatation). Here, we made the choice of setting a constant kinematic viscosity $\nu = \mu/\rho$. Indeed, it was found that the time step constraint for the simulation was driven by the viscous term, i.e. $\Delta t \sim \nu/\Delta x^2 = \mu/\rho \Delta x^2$ ($\Delta x$ is the spatial resolution). Hence, using a constant dynamic viscosity $\mu$ 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 $\rho_0 = \langle\rho\rangle$ (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' = \langle|\mathbf{u}|^2\rangle^{1/2}$ was equal to 1 (±1%) in all cases. The Reynolds number $Re = u'L/\nu \equiv 1/\nu$ and Mach number $Ma = u'/c_s \equiv 1/c_s$ were then varied independently by changing the kinematic viscosity $\nu$ and sound speed $c_s$, respectively. The database covers the range $263 \le Re \le 4166$ and $0.25 \le Ma \le 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 $\eta = 2\Delta_x$, where $\eta$ 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 \le 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 $\langle p\theta \rangle$. 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 $\theta$. 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 $\langle p\theta \rangle$ 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 $\eta \approx 4\Delta_x$.

Table 1 summarises some key quantities of the numerical database. The total number of points used in the simulation is denoted $N$. The energy injection $\epsilon_f = \langle f \cdot \mathbf{u} \rangle/\rho_0$ equals the kinetic energy dissipation rate since the flow is at steady state. The Kolmogorov length scale is given by $\eta = v^{3/4}/\epsilon_f^{1/4}$ and is given in units of $\Delta_x = L/\sqrt{N}$. The Taylor-scale Reynolds number $R_\lambda = u'^2\sqrt{5/(3\nu\epsilon_f)}$ varies between about 20 and 110. The pressure dilatation $\langle p\theta \rangle$ is given in units of $\epsilon_f$. The latter appears to be negative and thus represents a loss of kinetic energy. Its maximum value is about 5% of $\epsilon_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_\lambda$ 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 $\epsilon_d = 4\langle\mu\theta^2\rangle/3$ (Jagannathan & Donzis 2016; John et al. 2021). Table 1 reveals that $\epsilon_d$ is negligible for $Ma \le 0.5$ but increases up to 25% of $\epsilon_f$ at $Ma = 1$ and roughly 60% (if not more) of $\epsilon_f$ for $Ma \ge 2$. We note also that the ratio $\epsilon_d/\epsilon_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 $\epsilon_d/\epsilon_f$ when Re
<table><caption>**Table 1.** Description of the numerical database.</caption><thead><tr><th>Re</th><th>Ma</th><th>v</th><th>c<sub>s</sub></th><th>N</th><th>ε<sub>f</sub></th><th>η/Δ<sub>x</sub></th><th>R<sub>λ</sub></th><th>−⟨pθ⟩/ε<sub>f</sub></th><th>ε<sub>d</sub>/ε<sub>f</sub></th></tr></thead><tbody><tr><td>263</td><td>0.25</td><td>3.8 × 10<sup>−3</sup></td><td>4.00</td><td>128<sup>3</sup></td><td>0.96</td><td>1.98</td><td>21.4</td><td>0.48 %</td><td>0.31 %</td></tr><tr><td>263</td><td>0.50</td><td>3.8 × 10<sup>−3</sup></td><td>2.00</td><td>128<sup>3</sup></td><td>0.92</td><td>2.00</td><td>21.9</td><td>0.08 %</td><td>2.01 %</td></tr><tr><td>263</td><td>1.00</td><td>3.8 × 10<sup>−3</sup></td><td>1.00</td><td>128<sup>3</sup></td><td>1.03</td><td>1.95</td><td>20.7</td><td>1.72 %</td><td>20.21 %</td></tr><tr><td>263</td><td>2.00</td><td>3.8 × 10<sup>−3</sup></td><td>0.50</td><td>256<sup>3</sup></td><td>1.30</td><td>3.67</td><td>18.3</td><td>1.00 %</td><td>43.27 %</td></tr><tr><td>263</td><td>4.00</td><td>3.8 × 10<sup>−3</sup></td><td>0.25</td><td>256<sup>3</sup></td><td>1.27</td><td>3.67</td><td>18.6</td><td>0.59 %</td><td>50.44 %</td></tr><tr><td>714</td><td>0.25</td><td>1.4 × 10<sup>−3</sup></td><td>4.00</td><td>256<sup>3</sup></td><td>0.70</td><td>2.02</td><td>41.2</td><td>0.48 %</td><td>0.26 %</td></tr><tr><td>714</td><td>0.50</td><td>1.4 × 10<sup>−3</sup></td><td>2.00</td><td>256<sup>3</sup></td><td>0.68</td><td>2.04</td><td>41.7</td><td>0.09 %</td><td>2.12 %</td></tr><tr><td>714</td><td>1.00</td><td>1.4 × 10<sup>−3</sup></td><td>1.00</td><td>256<sup>3</sup></td><td>0.72</td><td>2.01</td><td>40.6</td><td>1.97 %</td><td>23.74 %</td></tr><tr><td>714</td><td>2.00</td><td>1.4 × 10<sup>−3</sup></td><td>0.50</td><td>512<sup>3</sup></td><td>1.01</td><td>3.70</td><td>34.3</td><td>1.26 %</td><td>50.44 %</td></tr><tr><td>714</td><td>4.00</td><td>1.4 × 10<sup>−3</sup></td><td>0.25</td><td>512<sup>3</sup></td><td>1.00</td><td>3.70</td><td>34.5</td><td>0.87 %</td><td>54.07 %</td></tr><tr><td>1886</td><td>0.25</td><td>5.3 × 10<sup>−4</sup></td><td>4.00</td><td>512<sup>3</sup></td><td>0.59</td><td>2.04</td><td>72.7</td><td>0.67 %</td><td>0.25 %</td></tr><tr><td>1886</td><td>0.50</td><td>5.3 × 10<sup>−4</sup></td><td>2.00</td><td>512<sup>3</sup></td><td>0.59</td><td>2.04</td><td>73.0</td><td>0.20 %</td><td>2.17 %</td></tr><tr><td>1886</td><td>1.00</td><td>5.3 × 10<sup>−4</sup></td><td>1.00</td><td>512<sup>3</sup></td><td>0.64</td><td>2.00</td><td>69.9</td><td>4.06 %</td><td>25.84 %</td></tr><tr><td>1886</td><td>2.00</td><td>5.3 × 10<sup>−4</sup></td><td>0.50</td><td>1024<sup>3</sup></td><td>0.82</td><td>3.76</td><td>61.8</td><td>3.25 %</td><td>57.43 %</td></tr><tr><td>1886</td><td>4.00</td><td>5.3 × 10<sup>−4</sup></td><td>0.25</td><td>1024<sup>3</sup></td><td>0.91</td><td>3.67</td><td>58.9</td><td>1.96 %</td><td>57.62 %</td></tr><tr><td>4166</td><td>0.25</td><td>2.4 × 10<sup>−4</sup></td><td>4.00</td><td>1024<sup>3</sup></td><td>0.62</td><td>2.23</td><td>106.1</td><td>0.68 %</td><td>0.13 %</td></tr><tr><td>4166</td><td>0.50</td><td>2.4 × 10<sup>−4</sup></td><td>2.00</td><td>1024<sup>3</sup></td><td>0.59</td><td>2.25</td><td>108.6</td><td>0.22 %</td><td>2.37 %</td></tr><tr><td>4166</td><td>1.00</td><td>2.4 × 10<sup>−4</sup></td><td>1.00</td><td>1024<sup>3</sup></td><td>0.59</td><td>2.25</td><td>108.7</td><td>4.84 %</td><td>26.10 %</td></tr><tr><td>4166</td><td>2.00</td><td>2.4 × 10<sup>−4</sup></td><td>0.50</td><td>1920<sup>3</sup></td><td>0.75</td><td>3.98</td><td>96.4</td><td>4.23 %</td><td>61.35 %</td></tr><tr><td>4166</td><td>4.00</td><td>2.4 × 10<sup>−4</sup></td><td>0.25</td><td>1920<sup>3</sup></td><td>0.82</td><td>3.88</td><td>91.7</td><td>2.60 %</td><td>61.43 %</td></tr></tbody></table>

**Notes.** The number of simulation points is noted *N* and the mesh resolution is denoted Δx. The kinetic energy injection rate per unit mass (= viscous dissipation) εf = ⟨f·u⟩/ρ₀ is given in units of *u*′³/*L* (= 1). The Taylor-scale Reynolds number Rλ = *u*′²√5/(3νεf). The Kolmogorov length scale is defined by η = ν³⁴/εf¹⁴ and is given in units of Δx = 1/√N. The pressure dilation ⟨pθ⟩ and the dilatation component of the kinetic energy dissipation rate εd = 4⟨μθ²⟩/3 are reported in percentage of εf.

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ᵢ/(*u*′³, where Lᵢ is the integral length scale. Although it is not the scope of the present work, note that the finite value of εf when v → 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.

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.

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

$$PDF(s) = \frac{1}{\sqrt{2\pi\sigma_s^2}} \exp\left(-\frac{1}{2}\left(\frac{s - \langle s \rangle}{\sigma_s}\right)^2\right), \quad (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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Two-dimensional slices of the density field $s = \ln \rho$ for increasing $Ma$ (from left to right) and increasing $Re$ (from top to bottom). The colourbar spans $\pm 3\sigma_s$.</figcaption></figure>

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. $\sigma_s^2$, with respect to $Ma$ is shown in Fig. 2b. Assuming a linear relationship between the density standard deviation $\sigma_\rho$ and $Ma$, viz. $\sigma_\rho = bMa$, and a Gaussian distribution for $s$ (i.e. a log-normal distribution for $\rho$), one arrives at (Federrath et al. 2008, 2010; Padoan & Nordlund 2011; Molina et al. 2012)

$$ \sigma_s^2 = \ln(1 + b^2 Ma^2) \quad (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 \approx 0.4$ for the mixed case used in the present work (Federrath et al. 2008, 2010). Fig. 2b indicates that the measured values of $\sigma_s^2$ gradually increase with Mach number. We further note that $\sigma_s^2$ 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 $\nu, 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

$$ \text{PDF}(|u|/c_s) = \sqrt{\frac{2}{\pi} \frac{(|u|/c_s)^2}{Ma_{1D}^3}} \exp\left(-\frac{1}{2} \frac{(|u|/c_s)^2}{Ma_{1D}^2}\right) \quad (5) $$

where $Ma_{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.

The fluctuations of some derivatives of $\rho$ and $u$ are also worth quantifying. These appear in particular in the transport equation for $\rho$ (or $s$) together with the transport equation for the kinetic energy explored later in this paper. In Figs. 4 and 5, we plot the
<figure><img src="image_3.png" /></figure>

**Fig. 2.** (a) Probability-density-function of $s = \ln\rho$ 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 $\sigma_s^2$ 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 $\sigma_s^2$ (b).

volume-weighted PDFs of $\nabla\rho$ and $\theta = \nabla \cdot \mathbf{u}$, respectively. Panels from a to e are for different Ma while the colour code indicate different Re.

Contrary to $s$ and $|\mathbf{u}|$, the fluctuations of the density gradient $\nabla\rho$ and dilatation $\theta$ depend on both Ma and Re. Increasing the Mach and Reynolds numbers yields larger tails in the PDFs for both $\nabla\rho$ and $\theta$. 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 $\nabla\rho$ and $\theta$ are highly intermittent. For instance, at Ma = 4 and Re = 4166, some local fluctuations of $\nabla\rho$ can exceed 100 times their standard deviation. Note that in isothermal flows, the pressure $p = \rho c_s^2$ and hence the PDFs portrayed in Fig. 4 would be the exact same for $\nabla p$.

The PDFs of dilatation $\theta$ (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 $\theta < 0$, have a much larger extent than the right wing which drops very quickly. Strong compression events are thus

<figure><img src="image_4.png" /></figure>

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

more likely than strong expansion events. Despite, the total volume where $\theta < 0$ is compensated by the volume formed by $\theta > 0$ so that $\langle\theta\rangle = 0$. At the lowest Ma, we note that the PDFs($\theta$) are rather symmetric with respect to the $\theta = 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($\theta$) at Ma = 0.8 and Ma = 1 (see their Fig. 2). Note that the x-coordinate axis in Fig. 5 is normalised by $\sigma_\theta$ which increases between Ma = 1 and Ma = 2 (this quantity will be discussed in some forthcoming paragraphs).

Our data further indicate that the PDFs($\theta$) 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($\theta$). This suggests that for Ma $\ge$ 2 the distribution of $\theta$ 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 $\nabla\rho$ and $\theta$, are denoted $\sigma_{\nabla\rho}^2$ and $\sigma_\theta^2$, respectively. Their evolution with respect to Re and Ma are reported in Figs. 6a–b. Noticeable is the very strong increase of $\sigma_{\nabla\rho}^2$ when Ma varies from 0.25 to 4 (Fig. 6a). Our data suggests a power-law variation of the form $\sigma_{\nabla\rho}^2 \sim Ma^4$. It also increases with the Reynolds number although with a lower exponent, viz. $\sigma_{\nabla\rho}^2 \sim Re^1$. The variance of the dilatation is presented in Fig. 6b. We first note a prompt increase of $\sigma_\theta$ for subsonic (though compressible) conditions before reaching a plateau for Ma $\ge$ 2, whose level depends only on the Reynolds number. Here again, we find that the values of $\sigma_\theta^2$ 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.3v0.3, M2v0.3, M3v0.3 in Scannapieco et al. (2024) as they correspond to roughly the same effective Reynolds number (Re $\sim$ 3400 according to our definition). Although quantitatively different, their data are in agreement with our finding that $\sigma_\theta^2/(u'/L)^2$ may saturate for Ma $\ge$ 2. Since dilatation appears in the vicinity of shocks (Wang et al.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Probability-density-function of the density gradient $\nabla \rho$ 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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 4 for the velocity divergence $\theta = \nabla \cdot u$. The dashed lines in panel (e) are the reproduction of the data for Ma = 2 (panel (d)).</figcaption></figure>

2011), the saturation of $\sigma_{\theta}^2$ 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 $\theta$ (open symbols in Fig. 6b). Note that since we have used a fluid with a constant kinematic viscosity $\nu$, the density-weighted variance of $\theta$, i.e. $\langle \rho\theta^2 \rangle / \rho_0 = 3\epsilon_d/(4\nu)$. Therefore, the saturation of $\langle \rho\theta^2 \rangle$ for Ma $\ge$ 2 means that the dilatational component of the kinetic energy dissipation $\epsilon_d$ also saturates to a certain value that uniquely depends on Re.

## 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.

$$ \langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle = \langle (\delta\boldsymbol{\rho}u) \cdot (\delta\boldsymbol{u}) \rangle \quad (6) $$

where $\delta(\rho u) = (\rho u)^+ - (\rho u)^-$ and $\delta u = u^+ - u^-$ represents the increment (the difference) of $\rho 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 $\langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle$ 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

$$ \lim_{r \to \infty} \langle |\delta \boldsymbol{u}|_{\rho}^2 \rangle = 2 \langle |\boldsymbol{u}|_{\rho}^2 \rangle = 2 \langle \rho |\boldsymbol{u}|^2 \rangle, \quad (7) $$

where $\langle \rho |\boldsymbol{u}|^2 \rangle$ is twice the turbulent kinetic energy. Hence, $\langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle$ 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 $\le r$. The transport equation for $\langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle$ is known as the Kármán-Howarth-Monin equation and can formally be written as (Galtier & Banerjee 2011; Hellinger et al. 2021b):

$$ d_t \langle |\delta \boldsymbol{u}|_{\rho}^2 \rangle = \langle \mathcal{T} \rangle + \langle \mathcal{P} \rangle + \langle \mathcal{V} \rangle + \langle \mathcal{F} \rangle \quad (8) $$

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

- $d_t \langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle$ denotes the time variations of $\langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle$. In statistically stationary flows, this term is zero.

- $\langle \mathcal{T} \rangle$ represents the “transport” of $\langle |\delta\boldsymbol{u}|_{\rho}^2 \rangle$ 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:

$$ \langle \mathcal{T} \rangle = - \langle \nabla_r \cdot (\delta \boldsymbol{u}) | \delta \boldsymbol{u}|_{\rho}^2 \rangle + \langle \mathcal{R} \rangle. \quad (9) $$

The leftmost term on RHS of Eq. (9) represents the transfer (or cascade) of kinetic energy among the different scales of
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Mach-number dependence of the variance of (a) the density gradient $\nabla\rho$ (in units of $\rho_0/L \equiv 1$) and (b) dilatation $\theta$ (in units of $u'/L \equiv 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.</figcaption></figure>

the flow. It writes as the divergence in scale-space $\nabla_r \cdot$ of the flux $(\delta\mathbf{u})|\delta\mathbf{u}|_\rho^2$. The term $\langle\mathcal{R}\rangle$ is present only in compressible flows. It represents a source term due to the dilatation and writes

$$ \langle \mathcal{R} \rangle = \left\langle (\delta \mathbf{u}) \cdot \left[ (\delta(\rho \mathbf{u}))(\bar{\delta} \theta) - (\bar{\delta}(\rho \mathbf{u}))(\delta \theta) \right] \right\rangle, \quad (10) $$

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

$\langle\mathcal{P}\rangle$ 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:

$$ \langle \mathcal{P} \rangle = 2 \langle (\delta P)(\delta \theta) \rangle - \langle C(-\nabla p) \rangle, \quad (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 $-\langle C(-\nabla P) \rangle$ arises due to variations of the density. The latter was derived by Hellinger et al. (2021b) and writes:

$$ \begin{aligned} C(\mathbf{a}) &:= (\delta\mathbf{u}) \cdot (\delta\mathbf{a}) - (\delta(\rho\mathbf{u})) \cdot (\delta\nu\mathbf{a}) \\ &= (\rho^+ v^- - 1)\mathbf{u}^+ \cdot \mathbf{a}^- + (\rho^- v^+ - 1)\mathbf{u}^- \cdot \mathbf{a}^+, \end{aligned} \quad (12) $$

where $v = 1/\rho$ is the specific volume and $\mathbf{a}$ denotes any vectorial field quantity. For the pressure term, we set $\mathbf{a} \equiv -\nabla p$.

It is easy to show that in constant-density flows, $\rho^\pm v^\mp = 1$ and hence $C(\mathbf{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 $\langle\mathcal{V}\rangle$ and $\langle\mathcal{F}\rangle$ 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.

$$ \langle \mathcal{V} \rangle = 2 \langle (\delta \mathbf{u}) \cdot (\delta (\nabla \cdot \mathbf{t})) \rangle - \langle C (\nabla \cdot \mathbf{t}) \rangle \quad (13a) $$

$$ \langle \mathcal{F} \rangle = 2 \langle (\delta \mathbf{u}) \cdot (\delta \mathbf{f}) \rangle - \langle C(\mathbf{f}) \rangle \quad (13b) $$

where $\mathbf{t}$ is the viscous stress tensor defined in Eq. (2) and $\mathbf{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

$$ \lim_{r \to \infty} \langle \mathcal{T} \rangle = 0, \quad (14a) $$

$$ \lim_{r \to \infty} \langle \mathcal{P} \rangle = 4 \langle p \theta \rangle, \quad (14b) $$

$$ \lim_{r \to \infty} \langle \mathcal{V} \rangle = -4\epsilon, \quad (14c) $$

$$ \lim_{r \to \infty} \langle \mathcal{F} \rangle = 4\epsilon_f, \quad (14d) $$

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

$$ d_t \frac{1}{2} \langle |\mathbf{u}|_\rho^2 \rangle = \langle p\theta \rangle - \epsilon + \epsilon_f, \quad (15) $$

is recovered from Eq. (8) in the limit of large separations ($\rho_0 = 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 $\epsilon = \epsilon_f$. The kinetic energy dissipation rate $\epsilon = \epsilon_s + \epsilon_d$, where, for homogeneous flows, $\epsilon_s = \langle\mu|\omega|^2\rangle$ ($\omega$ is the vorticity vector) and $\epsilon_d = 4\langle\mu\theta^2\rangle/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 $\langle|\delta\mathbf{u}|_\rho^2\rangle$ is presented in Figs. 7a and b, respectively. The separation $r$ is normalised using $L$ while $\langle|\delta\mathbf{u}|_\rho^2\rangle$ is divided by its large-scale asymptotic value, viz. $2\langle|\mathbf{u}|_\rho^2\rangle$. In Figs. 7a–b, the different sets of curves corresponding to either constant-Ma or constant-Re are shifted upwards for clarity.

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 $\langle|\delta\mathbf{u}|_\rho^2\rangle$ 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 $\langle|\delta\mathbf{u}|_\rho^2\rangle$ with respect to the scale $r$, i.e. $\langle|\delta\mathbf{u}|_\rho^2\rangle \sim r^\zeta$. 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:

$$ \zeta(r) = \frac{d \ln \langle |\delta \mathbf{u}|_\rho^2 \rangle}{d \ln r} \quad (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 $\langle \mathcal{F} \rangle$ 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.

Before we start analysing the KHM equation, it is worth recalling that $\langle |\delta u|^2_p \rangle$ 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 $\langle \mathcal{T} \rangle$ and pressure term $\langle \mathcal{P} \rangle$ 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 $\langle \mathcal{T} \rangle$ or $\langle \mathcal{P} \rangle$ 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 $\langle \mathcal{T} \rangle → 0$ and $\langle \mathcal{P} \rangle → 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 $\langle \mathcal{V} \rangle$ in Fig. 8 is classically presented in the form $\langle \mathcal{V} \rangle + 4ε_f$, and hence $\langle \mathcal{V} \rangle + 4ε_f → 4ε_f$ when $r → 0$ and $\langle \mathcal{V} \rangle + 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). Irrespec-

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Distribution of kinetic energy $\langle |\delta u|^2_p \rangle / 2 \langle |u|^2_p \rangle$ 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.</figcaption></figure>

tively 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 $\langle \mathcal{P} \rangle$ 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
<figure><img src="image_9.png" /><figcaption>**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\epsilon_f$ while the separation is divided by $L \equiv 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.</figcaption></figure>

values of Re are needed to perceive compressibility effects at such Mach numbers. The analysis of Wang et al. (2018) at $R_\lambda \approx 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 $\langle T \rangle = 4\epsilon_f$ for $\eta \ll r \ll L$.

For larger $Ma$ (panels b–d), the results are rather surprising. For instance, we observe that the contribution of the transport term $\langle T \rangle$ 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 $\langle P \rangle$ may become negligible at very large $Ma$. We also note that $\langle P \rangle$ 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 $\epsilon_f$ in the limit of large Reynolds numbers. We may thus extend the Kolmogorov 4/5-law to compressible flows by writing $\langle T \rangle + \langle P \rangle = 4\epsilon_f$ for $\eta \ll r \ll 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 $\langle P \rangle$ vanishes in the limit of large $Ma$ and the generalised Kolmogorov law may reduce to $\langle T \rangle = 4\epsilon_f$ for $\eta \ll r \ll L$ when $c_s \to 0$.

We also quantified the scale-by-scale distributions of the terms where dilatation $\theta$ contributes, i.e. the pressure-dilatation term $2\langle\delta\theta\delta p\rangle$ and the source term $\langle R \rangle$. 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 \ge 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 $\langle R \rangle$. The latter term $\langle R \rangle$ appears to increase in amplitude and reaches values close to $4\epsilon_f$ at the largest $Ma - Re$ condition investigated. We do not observe any saturation of this term when $Ma$ and
<figure><img src="image_10.png" /><figcaption>Fig. 9. Scale-by-scale contribution of the terms associated to the dilatation, i.e. $\langle \mathcal{R} \rangle$ and $2\langle \delta\theta\delta p \rangle$. The colour code is same as in Fig. 7a.</figcaption></figure>

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 $\langle \delta\theta\delta p \rangle$. Again, this indicates that the dilatation is limited by viscous effects at high $Ma$.

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 $\epsilon_d$ also saturates to a certain proportion of the total kinetic energy dissipation rate $\epsilon$. 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 dif-
ferent 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 condi-
tions.

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 trans-
fer 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 pro-
posed by Falkovich et al. (2010), i.e. $\langle T \rangle + \langle P \rangle = \epsilon_f$ may hold
in the limit of very large Reynolds numbers. For infinitely large
Mach numbers ($c_s \to 0$), our data also indicates that the pres-
sure term is likely to vanish and the 4/5-law may reduce to
$\langle T \rangle = \epsilon_f$. In supersonic situations, the non-linear transport term
$\langle T \rangle$ has a strong contribution due to dilatation, i.e. $\langle R \rangle$, which
represents a loss (gain) of kinetic energy for the small (large)
scales.

Observational data of the ISM often provides some mea-
sures of the velocity structure functions. The KHM is the the-
oretical framework to interpret these measurements in terms
of the underlying physical processes, such as energy injec-
tion, transfer, conversion, and dissipation. Our results indi-
cate that *Re* and *Ma* significantly influences these pro-
cesses. We believe that this insight can be useful for bet-
ter 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 plot-
ting data, comprises approximately 55 TB. These data are avail-
able from the authors upon request, subject to availability. Please
note that certain raw datasets may no longer be accessible due to
storage capacity limitations.

**Acknowledgements.** 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**

Alam, S., Sahu, P., & Verma, M. K. 2022, Phys. Rev. Fluids, 7, 074605

Aluie, H. 2013, Phys. D, 247, 54

Aluie, H., Li, S., & Li, H. 2012, ApJ, 751, L29

Antonia, R., & Burattini, P. 2006, J. Fluid Mech., 550, 175

Banerjee, S., & Galtier, S. 2013, Phys. Rev. E, 87, 013019

Banerjee, S., & Galtier, S. 2014, J. Fluid Mech., 742, 230

Banerjee, S., & Kritsuk, A. G. 2018, Phys. Rev. E, 97, 023107P

Beattie, J. R., Federrath, C., Klessen, R. S., Cielo, S., & Bhattacharjee, A. 2025, Nat. Astron., 1

Brahami, Y. 2020, PhD Theses (Normandie Université)

Danaila, L., Antonia, R., & Burattini, P. 2012, Phys. D, 241, 224

Donzis, D. A., & Jagannathan, S. 2013, J. Fluid Mech., 733, 221

Donzis, D. A., & John, J. P. 2020, Phys. Rev. Fluids, 5, 084609

Dubrulle, B. 2019, J. Fluid Mech., 867, P1

Elmegreen, B. G. 2008, Proc. Int. Astron. Union, 4, 289

Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211

Falgarone, E., & Puget, J.-L. 1995, A&A, 293, 840

Falgarone, E., Des Forêts, G. P., Hily-Blant, P., & Schilke, P. 2006, A&A, 452, 511

Falgarone, E., Pety, J., & Hily-Blant, P. 2009, A&A, 507, 355

Falkovich, G., Fouxon, I., & Oz, Y. 2010, J. Fluid Mech., 644, 465

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Federrath, C., Klessen, R. S., Iapichino, L., & Beattie, J. R. 2021, Nat. Astron., 5, 365

Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79

Federrath, C., Klessen, R. S., & Schmidt, W. 2009, ApJ, 692, 364

Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.-M. 2010, A&A, 512, A81

Federrath, C., Rathborne, J. M., Longmore, S. N., et al. 2017, in The Multi-Messenger Astrophysics of the Galactic Centre, eds. R. M. Crocker, S. N. Longmore, & G. V. Bicknell, IAU Symp. 322, 123

Federrath, C., Roman-Duval, J., Klessen, R., Schmidt, W., & Mac Low, M. M. 2022, Astrophysics Source Code Library [record ascl:2204.0001]

Ferrand, R., Galtier, S., Sahraoui, F., & Federrath, C. 2020, ApJ, 904, 160

Galtier, S., & Banerjee, S. 2011, Phys. Rev. Lett., 107, 134501P

Gerrard, I. A., Federrath, C., Pingel, N. M., et al. 2023, MNRAS, 526, 982

Glover, S. C., Federrath, C., Mac Low, M.-M., & Klessen, R. S. 2010, MNRAS, 404, 2

Godard, B., Falgarone, E., & Des Forêts, G. P. 2009, A&A, 495, 847

Hellinger, P., Papini, E., Verdini, A., et al. 2021a, ApJ, 917, 101

Hellinger, P., Verdini, A., Landi, S., et al. 2021b, Phys. Rev. Fluids, 6, 044607P

Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29

Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 1

Hily-Blant, P., & Falgarone, E. 2007, A&A, 469, 173

Hily-Blant, P., Falgarone, E., & Pety, J. 2008, A&A, 481, 367

Hopkins, P. F. 2013, Mon. Not. R. Astron Soc., 430, 1880
Jagannathan, S., & Donzis, D. A. 2016, J. Fluid Mech., 789, 669

John, J. P., Donzis, D. A., & Sreenivasan, K. R. 2021, J. Fluid Mech., 920, A20

Kolmogorov, A. 1941, Dokl. Akad. Nauk SSSR, 30, 301

Konstandin, L., Girichidis, P., Federrath, C., & Klessen, R. 2012, ApJ, 761, 149

Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416

Krumholz, M. R. 2014, Phys. Rep., 539, 49

Lai, C., Charonko, J., & Prestridge, K. 2018, J. Fluid Mech., 843, 382

Linsky, J. L., Redfield, S., Ryder, D., & Chasan-Taber, A. 2022, AJ, 164, 106

Mac Low, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Mohapatra, R., Federrath, C., & Sharma, P. 2021, MNRAS, 500, 5072

Molina, F., Glover, S., Federrath, C., & Klessen, R. 2012, MNRAS, 423, 2680

Myers, P. C. 1978, ApJ, 225, 380

Nie, Q., & Tanveer, S. 1999, Proc. Roy. Soc., 455, 1615

Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40

Pan, L., Ju, W., & Chen, J.-H. 2022, MNRAS, 514, 105

Pan, L., Padoan, P., & Nordlund, A. 2019a, ApJ, 876, 90

Pan, L., Padoan, P., & Nordlund, A. 2019b, ApJ, 881, 155

Passot, T., & Vázquez-Semadeni, E. 1998, Phys. Rev. E, 58, 4501

Pirozzoli, S., & Grasso, F. 2004, Phys. Fluids, 16, 4386

Rabatin, B., & Collins, D. C. 2023, MNRAS, 525, 297

Saffman, P. 1968, Topics in Nonlinear Physics, 485

Sakurai, Y., & Ishihara, T. 2023, Phys. Rev. Fluids, 8, 084606

Sakurai, Y., & Ishihara, T. 2024, Phys. Fluids, 36, 085152

Sakurai, Y., Ishihara, T., Furuya, H., Umemura, M., & Shiraishi, K. 2021, ApJ, 911, 140

Scannapieco, E., Pan, L., Buie, E., & Brüggen, M. 2024, Sci. Adv., 10, eado3958

Schmidt, W., & Grete, P. 2019, Phys. Rev. E, 100, 043116

Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127

Seta, A., & Federrath, C. 2021, Phys. Rev. Fluids, 6, 103701

Thiesset, F., & Vahé, J. 2025, J. Fluid Mech., submitted [arXiv:2509.11843]

Vazquez-Semadeni, E. 1994, Astrophys. J., 423, 681

Waagan, K., Federrath, C., & Klingenberg, C. 2011, J. Comput. Phys., 230, 3331

Wang, J., Shi, Y., Wang, L.-P., et al. 2011, Phys. Fluids, 23, 125103

Wang, J., Yang, Y., Shi, Y., et al. 2013, Phys. Rev. Lett., 110, 214505

Wang, J., Gotoh, T., & Watanabe, T. 2017a, Phys. Rev. Fluids, 2, 023401

Wang, J., Gotoh, T., & Watanabe, T. 2017b, Phys. Rev. Fluids, 2, 013403

Wang, J., Wan, M., Chen, S., & Chen, S. 2018, J. Fluid Mech., 841, 581

Yeung, P., Sreenivasan, K., & Pope, S. 2018, Phys. Rev. Fluids, 3, 064603