Compressible Turbulence in the Interstellar Medium: New Insights from a High-resolution Supersonic Turbulence Simulation

Understanding turbulent space and astrophysical plasmas is an ongoing physical challenge that has drawn a lot of attention throughout recent years. As these plasmas are present in a wide range of astrophysical media, from the Earth's magnetosphere to distant star-forming clouds, being able to properly describe them would allow us to make significant progress in understanding the physics controlling the shape and evolution of these media. Among these recent findings, it has been shown that the introduction of compressibility in magnetohydrodynamic (MHD) and Hall-MHD models of space plasmas leads to a higher estimate of the mean dissipation rate of total energy (used as a proxy for measuring plasma heating; Banerjee et al. 2016; Hadid et al. 2017; Andrés et al. 2019) compared to the incompressible case (Sorriso-Valvo et al. 2007). For the interstellar medium (ISM) where observations indicate that turbulence is supersonic (Wilson et al. 1970; Elmegreen & Scalo 2004; Heyer & Brunt 2004; Padoan et al. 2014; Krumholz & Federrath 2019), numerical simulations performed in the framework of compressible (isothermal) HD (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Kritsuk et al. 2007; Federrath et al. 2010; Federrath 2013) have shown the presence of filaments that resemble the structures observed in the ISM (Arzoumanian et al. 2011; Federrath 2016). It was found that incompressible predictions can be restored in some cases if one considers the density-weighted fluid velocity ${\rho }^{1/3}{\boldsymbol{u}}$ instead of the simple velocity ${\boldsymbol{u}}$ (Kritsuk et al. 2007; Schmidt et al. 2008), a behavior that can be understood dimensionally with an exact law (Galtier & Banerjee 2011).

What makes this kind of study difficult is that the mechanisms governing fluid turbulence are still not fully understood. Due to its chaotic nature, the favored way of studying turbulence is to resort to a statistical approach allowing for the use of specific tools (Frisch 1995), such as exact laws. Kolmogorov was among the pioneers in this field with his so-called four-fifths law (Kolmogorov 1941), an exact relation for homogeneous incompressible isotropic hydrodynamic (HD) turbulence that paved the way to new advances in the study of nonlinear physics. This statistical law allows one to express the mean rate of kinetic energy transfer per unit volume as a function of a two-point third-order longitudinal structure function in the limit of a high Reynolds number. In the wake of Kolmogorov's work, several other exact laws were derived for HD (Monin 1959; Antonia et al. 1997), quasi-geostrophic flows (Lindborg 2007), thin elastic plates (During & Krstulovic 2018), or plasmas (Politano & Pouquet 1998; Galtier 2008; Meyrand & Galtier 2010; Yoshimatsu 2012; Ferrand et al. 2019). The influence of compressibility being important for the description of space and astrophysical plasmas, efforts have been made during recent years to derive exact laws for compressible turbulence in HD (Galtier & Banerjee 2011; Banerjee & Kritsuk 2017; Lindborg 2019) and magnetohydrodynamics (MHD; Banerjee & Galtier 2013; Galtier 2016; Andrés & Sahraoui 2017; Andrés et al. 2018; Banerjee & Kritsuk 2018).

In the quest for exact laws for compressible turbulence, the complexity may increase significantly (Banerjee & Galtier 2013; Andrés et al. 2018). It is therefore relevant to ask whether it is possible to find a compact form of these laws that reveals the most salient feature of turbulence. Ferrand et al. (2019) recently showed that several different—yet equivalent—exact laws can be derived for Hall MHD, with some being more compact and easier to compute and interpret. Following the same method, it is the first main goal of this paper to demonstrate that such a compact form—called hereafter the generalized Kolmogorov law—exists in isothermal compressible HD turbulence. In the second part, this new relation is used to study such a turbulence, at integral Mach number 4, using a massive numerical simulation with a grid resolution of $10,{048}^{3}$ points (Federrath et al. 2016a, 2020). We proceed to a global computation of the law on the whole system and to a local computation along filamentary structures. Our analysis reveals supersonic turbulence properties that can be used to better understand ISM turbulence and star formation (Heyer & Brunt 2004; McKee & Ostriker 2007; Arzoumanian et al. 2011; Federrath & Klessen 2012; Padoan et al. 2014; Orkisz et al. 2017).

Section 2 contains the main steps of the derivation of the new exact law for compressible HD turbulence, along with a first theoretical interpretation. In Section 3 we apply this model to our numerical simulation and expose the results obtained. These results are discussed in Section 4 and we give an overall conclusion in Section 5 in the context of ISM turbulence and star formation.

Our analysis is based on the compressible HD equations

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\partial }_{t}(\rho {\boldsymbol{u}})+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{uu}})=-{\boldsymbol{\nabla }}P+{\boldsymbol{d}}+{\boldsymbol{f}},\end{eqnarray} \tag{ 2 }$$

where ρ is the density, P is the pressure, ${\boldsymbol{d}}\equiv \mu {\rm{\Delta }}{\boldsymbol{u}}+(\mu /3){\boldsymbol{\nabla }}\theta$ is the dissipation, $\theta \equiv {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}$ is the dilatation, μ is the coefficient of viscosity, and ${\boldsymbol{f}}$ is a stationary homogeneous external force assumed to act on large scales. The system is closed with the isothermal equation of state $P={c}_{s}^{2}\rho$ with c_s being the sound speed, assumed to be constant (in practice, we only use the density in the numerical code and set c_s = 1). The energy equation takes the form

$$\begin{eqnarray}&&{\partial }_{t}\langle E\rangle =\langle {\boldsymbol{u}}\cdot {\boldsymbol{d}}\rangle +\langle {\boldsymbol{u}}\cdot {\boldsymbol{f}}\rangle ,\end{eqnarray} \tag{ 3 }$$

with $\langle \,\rangle$ an ensemble average, $E=\rho {u}^{2}/2+\rho e$ the total energy, $e={c}_{s}^{2}\mathrm{ln}(\rho /{\rho }_{0})$ the internal energy (${\rho }_{0}=\langle \rho \rangle$ is the average density), $\langle {\boldsymbol{u}}\cdot {\boldsymbol{f}}\rangle =\varepsilon$ the mean rate of total energy injection into the system, and $\langle {\boldsymbol{u}}\cdot {\boldsymbol{d}}\rangle =-\mu \langle {({\rm{\nabla }}\times {\boldsymbol{u}})}^{2}\rangle -\displaystyle \frac{4}{3}\mu \langle {\theta }^{2}\rangle$.

We define ${\boldsymbol{\ell }}$ the spatial increment connecting two points x and x' as ${\boldsymbol{x}}^{\prime} ={\boldsymbol{x}}+{\boldsymbol{\ell }}$ and, for any given field ξ, $\xi ({\boldsymbol{x}})\equiv \xi$, and $\xi ({\boldsymbol{x}}^{\prime} )\equiv \xi ^{\prime}$. We follow the same idea as Hellinger et al. (2018) and Ferrand et al. (2019) and search for a dynamical equation of a structure function for the fluctuating energy:

$$\begin{eqnarray}&&\langle { \mathcal S }\rangle =\langle \bar{\delta }\rho {(\delta {\boldsymbol{u}})}^{2}+\delta \rho \delta e\rangle ,\end{eqnarray} \tag{ 4 }$$

where for any given field ξ, $\delta \xi \equiv \xi ^{\prime} -\xi$ and $\bar{\delta }\xi \equiv (\xi +\xi ^{\prime} )/2$. The use of this structure function represents the main difference between this approach and the one of Galtier & Banerjee (2011). Developing ${ \mathcal S }$ leads to

$$\begin{eqnarray}\begin{array}{rcl}\langle { \mathcal S }\rangle & = & 2\langle E\rangle -2\langle \bar{\delta }\rho {\boldsymbol{u}}\cdot {\boldsymbol{u}}^{\prime} \rangle -\langle \rho ^{\prime} e+\rho e^{\prime} \rangle \\ & & +\,\displaystyle \frac{1}{2}\langle \rho u{{\prime} }^{2}+\rho ^{\prime} {u}^{2}\rangle .\end{array}\end{eqnarray} \tag{ 5 }$$

Therefore, finding the temporal evolution of ${ \mathcal S }$ is akin to finding the temporal evolution of every term on the RHS of Equation (5).

The rest of the derivation is similar to the one given in Galtier & Banerjee (2011). Most of the terms have been derived in Galtier & Banerjee (2011); therefore, we only give the details for the calculation of the new term $\langle \rho u{{\prime} }^{2}\rangle$ (from which the symmetric contribution can be obtained immediately). We obtain the following expressions:

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\langle \rho {\boldsymbol{u}}\cdot {\boldsymbol{u}}^{\prime} \rangle ={{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot \langle -\rho ({\boldsymbol{u}}\cdot {\boldsymbol{u}}^{\prime} )\delta {\boldsymbol{u}}\\ \qquad +\,P{\boldsymbol{u}}^{\prime} -\rho e^{\prime} {\boldsymbol{u}}\rangle +\langle \rho \theta ^{\prime} ({\boldsymbol{u}}\cdot {\boldsymbol{u}}^{\prime} )\rangle \\ \qquad +\langle {\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}+{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}\cdot ({\boldsymbol{d}}^{\prime} +{\boldsymbol{f}}^{\prime} )\rangle ,\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\partial }_{t}\langle \rho e^{\prime} \rangle ={{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot \langle -\rho e^{\prime} \delta {\boldsymbol{u}}-P{\boldsymbol{u}}^{\prime} \rangle +\langle \rho e^{\prime} \theta ^{\prime} \rangle ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\langle \rho u{{\prime} }^{2}\rangle =\langle 2\rho {\boldsymbol{u}}^{\prime} \cdot {\partial }_{t}{\boldsymbol{u}}^{\prime} +u{{\prime} }^{2}{\partial }_{t}\rho \rangle \\ \ \ =\ \left\langle -{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{\nabla }}^{\prime} (\rho u{{\prime} }^{2})-2\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{\nabla }}^{\prime} P^{\prime} \right\rangle \\ \ \ +\,{{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot \langle \rho u{{\prime} }^{2}{\boldsymbol{u}}\rangle +\left\langle 2\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}^{\prime} +2\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}^{\prime} \right\rangle \\ \ \ =\ {{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot \langle -\rho u{{\prime} }^{2}{\boldsymbol{u}}^{\prime} +\rho u{{\prime} }^{2}{\boldsymbol{u}}\rangle +\left\langle \rho u{{\prime} }^{2}\theta ^{\prime} -2\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{\nabla }}^{\prime} P^{\prime} \right\rangle \\ \ \ +\,\left\langle 2\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}^{\prime} +2\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}^{\prime} \right\rangle ,\end{array}\end{eqnarray} \tag{ 8 }$$

where ${{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}$ is the derivative along the ${\boldsymbol{\ell }}$ direction. The combination of the different contributions gives after simplification,

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}\langle { \mathcal S }\rangle & = & 2{\partial }_{t}\langle E\rangle -{{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot \left\langle \bar{\delta }\rho {\left(\delta {\boldsymbol{u}}\right)}^{2}\delta {\boldsymbol{u}}\right\rangle \\ & & +\ \displaystyle \frac{1}{2}\langle (\rho \theta ^{\prime} +\rho ^{\prime} \theta ){\left(\delta {\boldsymbol{u}}\right)}^{2}\rangle -F-D,\end{array}\end{eqnarray} \tag{ 9 }$$

with, by definition,

$$\begin{eqnarray}\begin{array}{rcl}F & \equiv & \left\langle {\boldsymbol{u}}\cdot {\boldsymbol{f}}^{\prime} +{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}\cdot {\boldsymbol{f}}^{\prime} +\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}\right.\\ & & \left.-\,\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}\cdot {\boldsymbol{f}}-\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}^{\prime} \right\rangle ,\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}D & \equiv & \left\langle {\boldsymbol{u}}\cdot {\boldsymbol{d}}^{\prime} +{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}\cdot {\boldsymbol{d}}^{\prime} +\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}\right.\\ & & \left.-\,\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}\cdot {\boldsymbol{d}}-\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}^{\prime} \right\rangle .\end{array}\end{eqnarray} \tag{ 11 }$$

The stationarity assumption leads to the cancellation of the term on the LHS and the energy term on the RHS of Equation (9). In this situation, all the energy injected by the forcing must necessarily be dissipated at the same rate ε, so that we have the relationship $\varepsilon =\langle {\boldsymbol{u}}\cdot {\boldsymbol{f}}\rangle =-\langle {\boldsymbol{u}}\cdot {\boldsymbol{d}}\rangle$. The content of the forcing and dissipative terms F and D can then be broken down into three parts: first, since the forcing is assumed to act on large scales only, its variations across the simulation domain should remain small. Thus, forcing cross-terms like ${\boldsymbol{u}}\cdot {\boldsymbol{f}}^{\prime}$ are expected to behave like ${\boldsymbol{u}}\cdot {\boldsymbol{f}}=\varepsilon$ so we may write (see Kritsuk et al. 2013)

$$\begin{eqnarray}&&\left\langle {\boldsymbol{u}}\cdot {\boldsymbol{f}}^{\prime} +{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}\cdot {\boldsymbol{f}}^{\prime} +\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}\right\rangle \simeq 4\varepsilon .\end{eqnarray} \tag{ 12 }$$

Second, the stationarity assumption states that mean forcing and dissipation should balance each other with ${\boldsymbol{u}}\cdot {\boldsymbol{f}}=-{\boldsymbol{u}}\cdot {\boldsymbol{d}}=\varepsilon$, leading to

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}\cdot {\boldsymbol{f}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{f}}^{\prime} +\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}\cdot {\boldsymbol{d}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}^{\prime} \right\rangle \simeq 0.\end{eqnarray} \tag{ 13 }$$

Third, as the dissipation is assumed to act on small scales only dissipative cross-terms such as ${\boldsymbol{u}}\cdot {\boldsymbol{d}}^{\prime}$ are expected to be uncorrelated and of null statistical mean, hence (limit of small μ)

$$\begin{eqnarray}&&\left\langle {\boldsymbol{u}}\cdot {\boldsymbol{d}}^{\prime} +{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}+\displaystyle \frac{\rho }{\rho ^{\prime} }{\boldsymbol{u}}\cdot {\boldsymbol{d}}^{\prime} +\displaystyle \frac{\rho ^{\prime} }{\rho }{\boldsymbol{u}}^{\prime} \cdot {\boldsymbol{d}}\right\rangle \simeq 0.\end{eqnarray} \tag{ 14 }$$

With these different estimates, the generalized Kolmogorov law for three-dimensional compressible isothermal turbulence reads

$$\begin{eqnarray}&&-4\varepsilon ={{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot \left\langle \bar{\delta }\rho {(\delta {\boldsymbol{u}})}^{2}\delta {\boldsymbol{u}}\right\rangle -\displaystyle \frac{1}{2}\left\langle {(\rho \theta ^{\prime} +\rho ^{\prime} \theta )(\delta {\boldsymbol{u}})}^{2}\right\rangle .\end{eqnarray} \tag{ 15 }$$

Expression (15) is the first main result of this paper. This compact law is valid for homogeneous—but not necessarily isotropic—turbulence. As explained above, to derive this expression we have assumed the existence of an inertial range where the forcing and the dissipation are negligible (Frisch 1995; Aluie 2011). The expression found is much simpler than the one proposed in Galtier & Banerjee (2011) because (i) the flux, ${\boldsymbol{F}}\equiv \langle \bar{\delta }\rho {(\delta {\boldsymbol{u}})}^{2}\delta {\boldsymbol{u}}\rangle$, is constructed as a single term that resembles its incompressible version (that we recover by taking $\bar{\delta }\rho ={\rho }_{0}$, with ${\rho }_{0}$ a constant mass density), and (ii) the source is simply written as $S\equiv -{(1/2)\langle (\rho \theta ^{\prime} +\rho ^{\prime} \theta )(\delta {\boldsymbol{u}})}^{2}\rangle$. It is a purely compressible term, which goes to zero when the incompressible limit is taken; then, we recover the original form of the well-known four-thirds law (Antonia et al. 1997). The sign of S is directly given by the sign of the dilatation: when the flow is mainly in a phase of dilatation ($\theta \gt 0$) the source is negative, whereas in a phase of compression ($\theta \lt 0$) the source is positive. We can also see that the source tends to zero on small scales along with the ${(\delta u)}^{2}$ factor, as $(\rho \theta ^{\prime} +\rho ^{\prime} \theta )$ remains finite as ${\ell }$ goes to zero. This contrasts with the flux term that can still have a nontrivial contribution because of the ${\boldsymbol{\ell }}$–derivative introducing a $1/{\ell }$ scaling. A natural scale below which the source would be negligible is the sonic scale, i.e., the scale where the turbulence transitions from supersonic to subsonic (Federrath et al. 2010, 2020). Finally, note that expression (15) is Galilean invariant as the primitive Equations (1)–(2).

The generalized Kolmogorov law can be interpreted as if we had an effective cascade driven only by the flux term $-4{\varepsilon }_{\mathrm{eff}}\equiv {{\boldsymbol{\nabla }}}_{{\boldsymbol{\ell }}}\cdot {\boldsymbol{F}}$ (such that ${\varepsilon }_{\mathrm{eff}}=\varepsilon +S/4$) that involves the usual energy injection/dissipation rate (ε) known in incompressible theory and a new purely compressible component (source) that reflects contraction and dilatation of the turbulent structures. If we assume that the cascade driven by the flux term is direct (i.e., ${\varepsilon }_{\mathrm{eff}}\gt 0$) then a dilatation (compression) will tend to oppose (sustain) the energy cascade, preventing (enforcing) the formation of smaller structures. Furthermore, the dilatation of the structures ($S\lt 0$) can annihilate the cascade to small scales (if $\varepsilon =-S/4$) or even reverse it (if $\varepsilon +S/4\lt 0$) leading to the formation of large-scale structures via an inverse cascade. Note, however, that if S is scale dependent then compressible turbulence is not characterized by constant energy flux solutions as in incompressible theory (Kadomtsev & Petviashvili 1973; Passot et al. 1988). As we will see below, the numerical simulation will be very useful to go further in our interpretation.

In this section the exact law (15) derived above is used to investigate supersonic turbulence produced by a massive numerical simulation with a grid resolution of $10,{048}^{3}$ points and at Mach number 4. The Mach number is defined as ${ \mathcal M }={\sigma }_{v}/{c}_{{\rm{s}}}$ with ${\sigma }_{v}$ the velocity dispersion at the main forcing scale $L/2$ and L the simulation side length. The simulation was performed using a modified version of the FLASH code (Fryxell et al. 2000; Dubey et al. 2008; Federrath et al. 2020), solving the isothermal compressible HD equations in a triply periodic box. Following the methods in Federrath et al. (2010) the simulation uses a naturally mixed driving (ζ = 0.5) with an Ornstein–Uhlenbeck process acting on large scales. The forcing amplitude is a paraboloid spanning k = 1..3, peaking at k = 2 and reaching zero at both k = 1 and k = 3, where the wavenumber k is in units of $2\pi /L$. Thus, the forcing acts on scales larger than $L/3$, which are well above the ones we study in this paper. The data used here are seven snapshots of the $10,{048}^{3}$ simulation, sampled at a resolution of $2,{512}^{3}$, of the density and the three components of the velocity field, taken at 2, 3, 4, 5, 6, 7, and 8 turbulence turnover times T (downsampling the data eases its handling without affecting the results reported in this paper). Figure 1 shows through rms Mach number and minimum and maximum densities that statistics for both velocity and density have converged after two turnover times, indicating that the simulation has reached a statistically stationary state (Federrath et al. 2009), hence the use of snapshots for times $t\geqslant 2T$.

Zoom In Zoom Out Reset image size

For each snapshot the two terms of the exact law are computed along the three main axes x, y, and z, and then averaged spatially over the full box. The four signals obtained for different turnover times were eventually averaged to obtain the result displayed in Figure 2. First, we see that the mean rate of energy injection/dissipation ε (in red) is approximately constant over more than a decade. This observation indicates that the assumptions made to derive the law are satisfied on these scales of the simulation. The changes in ε observed on larger scales will be discussed in the next section in light of subsequent observations. Second, we see that the contribution of the flux term (in blue) is significantly higher than ε, which means that the source (in green) brings a correction to its contribution with an opposite sign, which is confirmed by the green dashed curve. For a better interpretation we can make a distinction between ε and ${\varepsilon }_{\mathrm{eff}}$ introduced above, i.e., the energy transferred between scales through the usual (incompressible) turbulence cascade driven by the flux term, in which case we have $\varepsilon ={\varepsilon }_{\mathrm{eff}}$. Here, we see that $\varepsilon \lt {\varepsilon }_{\mathrm{eff}}$ with a nonnegligible contribution from the source. This behavior contrasts with the one reported from direct numerical simulations of subsonic (compressible) MHD turbulence, where the overall contribution of the nonflux terms was found to be negligible with respect to the flux term (Andres et al. 2018). Note that in space plasma data, where it is not always possible to precisely measure the source (Hadid et al. 2017; Andrés et al. 2019), the variation of ${\varepsilon }^{{IHD}}$ in the inertial range may indicate the presence of nonnegligible compressible effects, especially in media where density fluctuations are high (Hadid et al. 2018).

Zoom In Zoom Out Reset image size

We further investigate the properties of this supersonic simulation in order to understand the origin of the source contribution and its influence on turbulence. As we have seen above, the source is globally positive, which reflects the dominance of compression. Therefore, for a given snapshot we searched for the grid point of minimal dilatation θ (maximum contraction of the fluid). In Figure 3 we show a slice of the data cube containing this point in the (yz) plane (other slices along (xy) or (xz) containing this point give a similar qualitative behavior). More precisely, we show the density-dilatation (top) and the modulus of the vorticity ${\boldsymbol{w}}={\rm{\nabla }}\times {\boldsymbol{u}}$ (bottom). These cuts reveal the existence of turbulent filamentary structures (elongated dark red structures for $\rho \theta$) in which both $| \theta |$ and $| {\boldsymbol{w}}|$ are up to several orders of magnitude higher than in the rest of the plane. A zoom on such structures is shown in Figure 4. These structures are typically delimited by very thin boundaries of strong contraction (dark blue lines for $\rho \theta$) in which a high turbulent activity with many vorticity tubes is observed. It is thought that the density-dilatation and the vorticity highlight the turbulence structures better than the previously used quantities θ and ρ (Kritsuk et al. 2007; Federrath et al. 2010), which are less relevant to investigate the physics involved in the generalized Kolmogorov law.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Regions with high density-dilatation are expected to drive most of the (average) source term. We therefore selected a sample of these regions (filaments) in which the mean density was high, and computed the source and flux term for increments along the main orientation of the filamentary structure. The results obtained were then averaged over the selected samples. In Figure 5 we show the result for a given filamentary structure and a "blank" region of weaker activity (top), and an average over eight filamentary structures (bottom). For all these filaments we observe a similar trend: the source is dominant, positive, and increases with spatial lag of the increments until it reaches approximately the sonic scale ${{\ell }}_{s}\simeq 0.01235L$, the scale where the scale-dependent Mach number is ${ \mathcal M }({{\ell }}_{s}/L)=1$ (see Federrath et al. 2020 for the original determination of ${{\ell }}_{s}$ in these simulations). The flux term does not have a constant sign but remains negligible with respect to the source on most scales. We note that, because of the sample size effect, the interpretation of the largest scales of the structures are subject to caution. A comparison with a blank region gives a quite different result: the values of both the flux and the source terms are up to five orders of magnitude lower than their counterparts in filamentary structures. One should note that we only evaluate here the specific contribution of small pieces of the flow. Consequently one would not expect to retrieve any form of theoretical scaling predicted by the exact law, which would only apply to the full statistical average on the whole simulation domain.

Zoom In Zoom Out Reset image size

We complement our analysis by taking the absolute value of the flux and the source before performing their statistical averages. Therefore, we define $\tilde{{\boldsymbol{F}}}=\left\langle \bar{\delta }\rho {(\delta {\boldsymbol{u}})}^{2}| \delta {\boldsymbol{u}}| \right\rangle$ and $\tilde{S}=\displaystyle \frac{1}{2}\left\langle | \rho \theta ^{\prime} +\rho ^{\prime} \theta | {(\delta {\boldsymbol{u}})}^{2}\right\rangle$ and compute them on the entire simulation domain at a given time. These two quantities represent the total activity due to the flux and the source, respectively, disregarding the sign of the local contributions and thus the direction of the resulting turbulent cascade (direct or inverse). Note that, again, these results do not have to comply to any theoretical prediction brought by the exact law, as the terms computed here are not the ones forming the exact law per se. Yet, the nonsigned quantities have the advantage of converging faster than their signed counterpart, and can lend some information about the mechanisms dominating on different scales. The results are reported in Figure 6. First, a clear power law $\tilde{S}\sim {{\ell }}^{1/2}$ emerges over two decades for the compressible source. By dimensionally equalizing the flux term and the source, we find $\delta u\sim {{\ell }}^{1/2}$ (we do not include the density that appears as a local average and not as a pure fluctuating quantity). This scaling is actually compatible with the one reported in Federrath et al. (2020) on supersonic scales using the second-order structure function (which is positive definite and therefore comparable to our calculation using absolute values), while a classical (incompressible) scaling $\delta u\sim {{\ell }}^{1/3}$ was approximately found on subsonic scales. Note that the supersonic law is dimensionally compatible with the velocity spectrum ${E}^{u}\sim {k}^{-2}$, a scaling often attributed to a purely compressible (Burgers) turbulence (Frisch 1995; Federrath 2013). In the framework of the generalized Kolmogorov law we see, however, that the change of slope reported previously can find a precise origin: it marks a transition toward a regime/scale where the absolute activity of the source becomes nonnegligible. Second, we see that the modified flux exhibits a plateau on small scales, as expected for a subsonic turbulent cascade mainly driven by the flux, which means that most of the energy transiting in either direction through these scales is transferred by a flux-driven process. A transition appears around the sonic scale above which the flux starts to drop: this behavior can be attributed to the dominance of the compressible source activity on supersonic scales.

Zoom In Zoom Out Reset image size

Based on a single simulation, realized, however, at an unprecedented spatial resolution, some conclusions can be drawn. By directly applying the new exact law derived analytically to the data we found that the amplitude of the source in the turbulent filamentary structures (Figure 5) is much higher (up to two orders of magnitude) than when it is computed on the whole simulation (Figure 2). As the overall turbulent activity is less intense in the other regions we can conclude that these filaments drive the global behavior of S. On the contrary, it has been impossible to identify a recurring behavior in specific parts of the system for the flux term as we did for the source. At intermediate, transonic scales both the flux and the source contributions reach a peak. Furthermore the sign of S in both the global and local computations is positive, leading to a value of ${\varepsilon }_{\mathrm{eff}}$ higher than ε, which is fixed externally by the forcing. We thus suggest that the energy cascade in supersonic HD turbulence reaches its maximum efficiency (i.e., ${\varepsilon }_{\mathrm{eff}}$ is maximal) around the transition from the subsonic to supersonic regimes. This efficiency decreases with scales such that ${\varepsilon }_{\mathrm{eff}}$ tends to be closer to ε on subsonic scales, where S becomes subdominant. On supersonic scales, however, the exact law shows a decrease in the energy cascade rate ε, which is no longer constant. The law being exact under several assumptions means that at least one of them is not verified. For example, we cannot exclude a nonlocal effect of the forcing that would modify the scaling law mostly on large scales. An anomalous dissipation on supersonic scales originating from the irregularities of the fields, which is not accounted for in our theory, could also contribute to the energy budget (Duchon & Robert 2000; Saw et al. 2016; Galtier 2018). On the other hand, Figure 6 shows that the total, nonsigned activity of the compressible source grows higher than that of the flux: this suggests that density-dilatation acts strongly on supersonic scales, but in such a way that local contributions of opposite signs cancel each other out in the scope of the exact law. This strong compressible activity, coupled to the flux activity becoming nonconstant at supersonic scales, suggests a transition between two regimes around the sonic scale: a subsonic turbulent cascade driven by both the compressible flux and the source, and a highly compressible regime at supersonic scales that does not feature a conservative cascade. The vorticity distribution shown in Figure 3 reinforces this interpretation since stronger vorticity is only observed inside the small-scale turbulent structures. Similar conclusions were drawn by Aluie et al. (2012) who reported for subsonic and transonic simulations that pressure-dilatation acts essentially on large scales, whereas at small scales, below a transitional "conversion" scale range, a conservative cascade appears. Here we shed a new light onto those findings using a different approach that helps better understand how various mechanisms shape supersonic turbulence.

A spectral analysis of the velocity field ${\boldsymbol{u}}$ (not shown in this paper) reveals the existence of two distinct power-law scalings, separated roughly by the sonic scale. On supersonic scales the scaling is close to ${k}^{-2}$, which is usually attributed to Burgers turbulence (Kadomtsev & Petviashvili 1973; Passot et al. 1988; Frisch 1995); on subsonic scales the scaling becomes close to ${k}^{-3/2}$, which is compatible with a theory of weak acoustic turbulence (Zakharov & Sagdeev 1970; L'vov et al. 1997). This provides other evidence of the existence of two distinct regimes at supersonic and subsonic scales: a more shock-driven compressible regime and a (possibly acoustic) turbulent regime. If we assume that the cascade rate computed here is representative of (if not identical to) the energy dissipation rate in the system, the observation that cascade/dissipation rate peaks near the sonic scale (Figure 5), where turbulence transitions from shock-like (${k}^{-2}$) to fluctuation/vortex-like (${k}^{-3/2}$) dominated regimes would be an indication of strong shock dissipation at that scale.

It is interesting to note that using a model of acoustic turbulence with weak shocks at Mach number close to unity, Lindborg (2019) reported that a scaling relation similar to the incompressible Kolmogorov law could be retrieved for a modified energy cascade rate. In the framework of said model, our scaling relation (15) still holds considering a similar modified energy cascade rate. This remark, along with the previous one, means that subsonic turbulence could be composed of a mixture of weak shocks, acoustic waves, and vortices.

A final remark can be made about the exact law: given the high resolution of the simulation one would expect the energy cascade rate to form a steady plateau over more than one decade. This small inertial range may be attributed to two possible effects not included in our exact law (15): (i) nonlocal effects due to the large-scale forcing and (ii) additional local dissipation (in the supersonic range) through shocks/discontinuities (Duchon & Robert 2000; Saw et al. 2016; Galtier 2018), since our exact law assumes smoothness of the turbulent fields. This shortcoming calls for a new theory of compressible HD turbulence where such singular fields and nonlocal effects due to large-scale forcing can be accounted for, and which would be very relevant to supersonic turbulence. In addition, it would be interesting to investigate the question of intermittency in supersonic turbulence by evaluating separately the contributions of the flux and the source. Such theories are beyond the scope of this paper and are left to future studies.

The theory developed in this paper and applied to high-resolution numerical simulations allows us to gain deep insight into supersonic ISM turbulence. The filamentary structures observed in the ISM seem to be characterized by a universal thickness of the order of the sonic scale (Arzoumanian et al. 2011; Federrath et al. 2016b). Their shape is supposed to be mainly due to HD turbulence and to be little affected by other factors such as gravity or magnetic fields (Federrath et al. 2016b; Ntormousi et al. 2016). These studies associated with our work suggest that this universality could be explained by the existence of the two distinct regimes reported here: (i) a supersonic regime dominated by shock-like structures where the energy cascade rate ε is not constant; (ii) a subsonic regime with a lower and mainly constant ε, where vortices (and acoustic waves) are produced and in which a classic conservative cascade is formed. In between, the transonic scales where turbulence reaches its peak of effective energy transfer would correspond to the size of the filaments. Our interpretation is thus that filaments are stuck on the smallest scale of the supersonic regime, which is the sonic scale, while the weaker subsonic cascade produces vorticity tubes on smaller scales.

Applications of the law to more complete simulations, featuring for instance gravitational forces or magnetic fields, would help refine this interpretation and may provide new clues on the interplay between ISM turbulence and the problem of star formation (mac Low & Klessen 2004; Hennebelle & Falgarone 2012; Padoan et al. 2014). For example, Orkisz et al. (2017) were able for the first time to observationally derive the fractions of momentum density contained in the solenoidal and compressive modes of turbulence. It was in the Orion B molecular cloud where the mean Mach number was ∼6. They showed that the compressive modes are dominant in regions with a high star formation rate (as predicted in Federrath & Klessen 2012). According to our analysis the source term is the dominant component of compressible turbulence inside filaments. Future work that would directly link the formalism of exact laws (and thus the source and flux terms) to the star formation rate could significantly advance our understanding of how turbulence controls the formation of structures on different scales in the ISM.

We thank the anonymous referee for valuable suggestions, which helped to improve this work. C.F. acknowledges funding provided by the Australian Research Council (Discovery Project DP170100603 and Future Fellowship FT180100495), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi, and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. The simulation software FLASH was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago.