Kinetic and internal energy transfer in implicit large eddy simulations of forced compressible turbulence

TL;DR

This study investigates energy transfer mechanisms in compressible turbulence using simulations with cooling, revealing how Mach number influences energy flux components and scaling behaviors.

Contribution

It introduces a simulation approach with cooling to analyze energy transfer across Mach numbers, highlighting the impact of forcing modes and identifying regimes dominated by advection.

Findings

Energy transfer rates vary with Mach number and forcing composition.

High Mach number turbulence exhibits advection-dominated energy transfer.

Scaling exponents align with Burgers and Kolmogorov theories depending on variables.

Abstract

We revisit the problem of how energy transfer through the turbulent cascade operates in compressible hydrodynamic turbulence. In general, there is no conservative compressible cascade since the kinetic and internal energy reservoirs can exchange energy through pressure dilatation. Moreover, statistically stationary turbulence at high Mach number can only be maintained in nearly isothermal gas, i.e., if excess heat produced by shock compression and kinetic energy dissipation is continuously removed from the system. We mimic this process by a linear cooling term in numerical simulations of turbulence driven by stochastic forcing. This allows us to investigate turbulence statistics for a broad range of Mach numbers. We compute the rate of change of kinetic and internal energy in wavenumber shells caused by advective, compressive, and pressure dilatation effects and constrain power-law fits to compressible turbulence energy spectra to a range of wavenumbers in which the total energy transfer is close to zero. The resulting scaling exponents are significantly affected by the forcing. Depending on the root mean square Mach number, we find a nearly constant advective component of the cross-scale flux of kinetic energy at intermediate wavenumbers for particular mixtures of solenoidal and compressive modes in the forcing. This suggests the existence of a natural, Mach number dependent mixture of forcing modes. Our findings also support an advection-dominated regime at high Mach numbers with specific scaling exponents (Burgers scaling for the pure velocity fluctuation $u$ and Kolmogorov scaling for the mass-weighted variable $v=\rho^{1/3}u$).

Tables

Table 1. Table 3: Power-law fits for statistically stationary turbulence in the estimated inertial range [ k min , k max ] subscript 𝑘 min subscript 𝑘 max [k_{\rm min},k_{\rm max}] following from the criteria formulated in Section  LABEL:sec:numericalTransfer (gray shaded wavenumber ranges in Fig.s  LABEL:fig:transf_alpha , LABEL:fig:cross_scale_2to1 , and  LABEL:fig:cross_scale_vel4 ). The slopes of the energy spectra E u ​ ( k ) subscript 𝐸 𝑢 𝑘 E_{u}(k) , E v ​ ( k ) subscript 𝐸 𝑣 𝑘 E_{v}(k) , and E w ​ ( k ) subscript 𝐸 𝑤 𝑘 E_{w}(k) are given by the exponents β u subscript 𝛽 𝑢 \beta_{u} , β v subscript 𝛽 𝑣 \beta_{v} , and β w subscript 𝛽 𝑤 \beta_{w} , respectively. For details about the fitting procedure see Appendix  A .

$V$	$\zeta$	$k_{\min }$	$k_{\max }$	${\beta }_u$	${\beta }_v$	${\beta }_w$
0.25	1	9.5	45.2	-1.667	-1.663	-1.661
1.0	1	9.5	22.6	-1.766	-1.672	-1.611
1.0	2/3	9.5	19.0	-1.826	-1.721	-1.651
2.0	2/3	9.5	32.0	-1.792	-1.542	-1.391
2.0	1/2	9.5	32.0	-1.772	-1.528	-1.387
2.0	1/4	9.5	22.6	-1.867	-1.684	-1.584
4.0	2/3	22.6	38.1	-2.110	-1.587	-1.337
4.0	1/2	22.6	38.1	-2.101	-1.609	-1.363
4.0	1/4	13.5	32.0	-1.925	-1.578	-1.402
4.0	1/8	11.3	32.0	-1.936	-1.647	-1.492

Equations

$${\partial }_t\rho +{\partial }_i(\rho u_i)$$

$${\partial }_t(\rho u_i)+{\partial }_j(\rho u_i{u}_j)+{\partial }_{i}p$$

$${\partial }_t(\rho e)+{\partial }_j(\rho u_{j}e)+p{\partial }_j{u}_j$$

$$\int \underset{\equiv E_{\mathrm{kin}}(\mathbf{x})}{\underbrace{\frac{1}{2}{w}_i{w}_i}}\mathrm{d}\mathbf{x}=\frac{1}{{\left(2\pi \right)}^3}\int \underset{\equiv E_{\mathrm{kin}}(\mathbf{k})}{\underbrace{\frac{1}{2}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}w]\bigwedge 0.5ex\text{stackon}[1pt]w{\text{tmpbox}}_i\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}w]\bigwedge 0.5ex\text{stackon}[1pt]w{\text{tmpbox}}_i^{\ast }}}\mathrm{d}\mathbf{k}$$

$$\begin{array}{rl}{\partial }_t{E}_{\mathrm{kin}}(\mathbf{x})=& -w_i{u}_j{\partial }_j{w}_i-\frac{1}{2}{w}_i{w}_i{\partial }_j{u}_j\\ & -\frac{w_i}{\sqrt{\rho }}{\partial }_{i}p+w_i\sqrt{\rho }{a}_i.\end{array}$$

$${\partial }_t{E}_{\mathrm{kin}}(\mathbf{k})=\mathrm{Re}\left[-\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}{w}_i]\bigwedge 0.5ex\text{stackon}[1pt]w_i\text{tmpbox}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}{u}_j{\partial }_j{w}_i]\bigwedge 0.5ex\text{stackon}[1pt]u_j{\partial }_j{w}_i{\text{tmpbox}}^{\ast }-\frac{1}{2}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}{w}_i]\bigwedge 0.5ex\text{stackon}[1pt]w_i\text{tmpbox}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}{w}_i{\partial }_j{u}_j]\bigwedge 0.5ex\text{stackon}[1pt]w_i{\partial }_j{u}_j{\text{tmpbox}}^{\ast }-\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}{w}_i]\bigwedge 0.5ex\text{stackon}[1pt]w_i\text{tmpbox}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}\frac{1}{\sqrt{\rho }}{\partial }_{i}p]\bigwedge 0.5ex\text{stackon}[1pt]\frac{1}{\sqrt{\rho }}{\partial }_{i}p{\text{tmpbox}}^{\ast }+\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}{w}_i]\bigwedge 0.5ex\text{stackon}[1pt]w_i\text{tmpbox}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}\sqrt{\rho }{a}_i]\bigwedge 0.5ex\text{stackon}[1pt]\sqrt{\rho }{a}_i{\text{tmpbox}}^{\ast }\right]$$

$$E_{\mathrm{int}}=\frac{p}{\gamma -1}=\frac{\rho c_{\mathrm{s}}^2}{\gamma (\gamma -1)}$$

$$q=\sqrt{\rho }{c}_{\mathrm{s}}$$

$$E_{\mathrm{int}}=\frac{q^2}{\gamma (\gamma -1)}$$

$$E_{\mathrm{tot}}=E_{\mathrm{kin}}+E_{\mathrm{int}}=\frac{|\mathbf{w}{|}^2}{2}+\frac{q^2}{\gamma (\gamma -1)}.$$

$$\begin{array}{rl}{\partial }_t{q}^2=& -2q{u}_j{\partial }_{j}q-q^2{\partial }_j{u}_j\\ & -(\gamma -1)q^2{\partial }_j{u}_j+\gamma (\gamma -1){\Lambda }_{\mathrm{HC}}.\end{array}$$

$$E_v(k)\Delta \mathrm{K}\simeq {\int }_{\mathrm{K}}\frac{1}{2}\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}v]\bigwedge 0.5ex\text{stackon}[1pt]v{\text{tmpbox}}_i\text{savestack}\text{tmpbox}\text{stretchto}\text{scaleto}\text{scalerel}\ast [\text{widthof}v]\bigwedge 0.5ex\text{stackon}[1pt]v{\text{tmpbox}}_i^{\ast }\mathrm{d}\mathbf{k}$$

$${\int }_{k_{\mathrm{s}}}^{\infty }2E_u(k)\mathrm{d}k=c_{\mathrm{s}}^2,$$

Full text

††preprint: APS/123-QED

I Introduction

Kolmogorov’s theory of incompressible isotropic turbulence states that the flux of kinetic energy across a given wavenumber is constant for wavenumbes in the inertial subrange in which both large-scale features of the flow (stirring, boundary conditions, etc.) and viscous dissipation are negligible Frisch (1995); Alexakis and Biferale (2018). An important consequence is the two-thirds law for the scaling of the squared velocity fluctuation, corresponding to the famous $k^{-5/3}$ energy spectrum. However, it has proven to be difficult to carry over the concepts of this theory to compressible turbulence. A key problem is that, in general, kinetic energy does not cascade conservatively because it can be converted into internal energy and vice versa through pressure dilatation. An important step in addressing this problem was the demonstration of scale-locality for the transfer of kinetic energy in compressible turbulence by Aluie Aluie (2011). Under certain assumptions (sufficiently fast decline of the pressure dilatation cospectrum), it can be shown that the exchange of kinetic and internal energy becomes subdominant and the two energy budgets are decoupled on sufficiently small scales. In this case, a compressible turbulent cascade with constant kinetic energy flux can exist. At about the same time, the generalization of the relation between the energy flux and two-point-correlation functions for the fully compressible case was analytically derived by Galtier & Banerjee Galtier and Banerjee (2011). In addition to a term that reduces to the Kolmogorov result in the incompressible limit, there are contributions related to compression effects (non-zero divergence of the flow) and deviations from a conservative kinetic energy cascade.

A further difficulty is the ambiguity in the scale-decomposition of energy density and fluxes in the compressible case. For example, the energy density in spectral space can be based on the Fourier transforms of the primitive velocity variable $u$ and momentum density $\rho u$ Graham et al. (2010). Alternatively, the energy density can be considered as square of the variable $w=\sqrt{\rho}\,u$ Kida and Orszag (1990) in analogy to the incompressible case (where the energy is simply given by $u^{2}$). For energy flux, there are even more variants Alexakis and Biferale (2018). For example, two-point correlation functions of velocity and momentum increments inspired by the derivation of the four-fifth law are used in Galtier and Banerjee (2011); Wagner et al. (2012); Kritsuk et al. (2013), while spatial coarse-graining results in expressions linked to subgrid-scale terms in large eddy simulations Eyink and Aluie (2009); Aluie (2013); Zhao and Aluie (2018). In addition, different flavors of spectral decomposition are applied to compute energy transfers between scales (basically, the transfer is the rate of change of the energy flux within a range of wave numbers). This is the most common method to analyze numerical data Alexakis et al. (2005); Graham et al. (2010); Mininni (2011); Grete et al. (2017); Verma (2019).

Numerical experiments play a crucial role in studying the phenomenology of compressible turbulence. While it remains challenging to produce strongly compressible or even supersonic turbulent flows under laboratory conditions, the regime of high Mach numbers is accessible in massively parallelized simulations on supercomputers. In astrophysics, substantial interest into supersonic turbulence was triggered by the theoretical investigation of star formation in interstellar gas clouds Mac Low and Klessen (2004); Hennebelle and Falgarone (2012). In a landmark work on scaling properties of supersonic isothermal turbulence Kritsuk et al. (2007) the hypothesis was put forward that mass-weighing of the velocity with $\rho^{1/3}$ allows for the extension of Kolmogorov scaling into the compressible regime. This idea was to some extent questioned in Schmidt et al. (2008); Federrath et al. (2010); Federrath (2013). However, the results of these studies remained inconclusive, as only particular Mach numbers were investigated and turbulence was driven by either solenoidal or compressive (dilatational) forcing modes. Further examples for numerical studies of driven hydrodynamic turbulence are Schmidt et al. (2009); Price et al. (2011); Aluie et al. (2012); Konstandin et al. (2012a); Wagner et al. (2012); Kritsuk et al. (2013).

In this article, we pick up ideas formulated in Kida and Orszag (1990); Alexakis et al. (2005) to compute the shell-to-shell transfer of kinetic energy in compressible turbulence based on the velocity variable $w=\sqrt{\rho}\,u$ (Section II). By defining $q=\sqrt{\rho}\,c_{\mathrm{s}}$, where $c_{\mathrm{s}}$ is the speed of sound, we incorporate internal energy transfers along the same lines as for magnetic energy transfers in compressible magnetohydrodynamic (MHD) turbulence with $B=\sqrt{\rho}\,v_{\rm a}$ ($v_{\rm a}$ is the Alfvénic velocity) as basic variable for the spectral decomposition Grete et al. (2017). This also guarantees positivity of the spectral energy densities. A similar approach was recently put forward by Mittal and Girimaji (2019). Our goal is to investigate whether there is a compressible inertial subrange defined by zero energy transfer (corresponding to constant flux) and how it is affected by pressure-dilatation effects for Mach numbers ranging from transonic to supersonic (Sections LABEL:sec:numericalTransfer and LABEL:sec:crossFlux). Moreover, we investigate the impact of different mixtures of solenoidal and compressive modes in the forcing. Based on our results we attempt to determine scaling exponents of turbulent energy spectra with different mass-weighing (Section LABEL:sec:numericalSpect).

Since time-averaging over many snapshots of the flow is paramount to obtain meaningful statistics, we need to maintain a statistically stationary state over several dynamical timescales. As in most simulations of compressible turbulence mentioned above, we perform implicit large eddy simulations, where dissipation results from numerical truncation terms, instead of explicitly treating viscous dissipation (see Schmidt (2015) and references cited therein). Explicit physical viscosity (bulk and shear) is neglected since Reynolds numbers in typical astrophysical systems are by far too high to allow for direct numerical simulations. Regardless of the nature of dissipation, the mean heating rate due to the dissipation of kinetic energy approximately equals the rate of energy injected by the forcing. To avoid a declining Mach number of supersonic turbulence, excess heat has to be continuously removed from the system. This is commonly achieved in numerical simulations of driven turbulence either by applying an isothermal equation of state (i.e. setting the speed of sound $c_{\mathrm{s}}$ indentical to a constant, which reduces the system of dynamic variables to mass density and momentum) or by artificially increasing the internal energy by setting the adiabatic exponent to a value that differs by a tiny fraction from unity (in this case, only the pressure is physical). Both methods are not suitable for our study because we need a physically meaningful (i.e. neither constant nor artificially increased) internal energy to compute transfers between the kinetic and internal energy reservoirs. For this reason, we define internal energy on the basis of the ideal gas equation and apply a simple toy model to mimic radiative cooling of the gas with finite cooling time. Details about the numerical simulations are outlined in Section LABEL:sec:stat and Appendix A. The outcome of our analysis is discussed in Section IV.

II Theory and methods

II.1 Dynamical equations

The equations of compressible gas dynamics in conservative form are given by

[TABLE]

The density is denoted by $\rho$, the velocities by $\mathbf{u}$, the internal energy density by $\rho e$, and the thermal pressure by $p=(\gamma-1)\rho e$, where $\gamma$ is the adiabatic exponent of the gas. The source term $\mathbf{f}=\rho\mathbf{a}$ on the right hand side of the momentum equation (2) represents an external force density and $\Lambda_{\rm HC}$ on the right hand side of the internal energy equation (3) a heating and cooling function (see Section LABEL:sec:numericalData for the modeling of these terms in our simulations).

II.1.1 Decomposition of kinetic energy

As discussed in Grete et al. (2017), there is no unambiguous definition of the spectral kinetic energy density $E_{\mathrm{kin}}$ in the compressible case. Since the symmetrized expression $\mathrm{Re}\left[\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{u}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{u}{\tmpbox}_{i}\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{\rho u}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\rho u}{\tmpbox}_{i}^{*}\right]/2$ used in Graham et al. (2010) does not guarantee positive definiteness in wavenumber space, we use for our analysis kinetic energy densities defined by

[TABLE]

where $\mathbf{w}\equiv\sqrt{\rho}\,\mathbf{u}$ Kida and Orszag (1990).

The corresponding dynamical equations are derived in Kida and Orszag (1990) for the hydrodynamic case and in Grete et al. (2017) for the general MHD case. In the hydrodynamic case, the dynamical equation for the kinetic energy density in real space reduces to

[TABLE]

An equivalent representation of the energy equation in terms of the mass-weighted variable $\mathbf{v}=\rho^{1/3}\,\mathbf{u}$ cannot be obtained, although this variable seems to be the natural choice for flux terms on the right-hand side Kritsuk et al. (2007).

In wavenumber space, the kinetic energy density is given by

[TABLE]

where wide hats signifiy Fourier transforms and an asterisk complex conjugation.

II.1.2 Decomposition of internal energy

The internal energy density $E_{\mathrm{int}}=\rho e$ can be expressed in terms of the speed of sound $c_{\mathrm{s}}$ as

[TABLE]

This suggests a decomposition of the internal energy density in terms of the variable

[TABLE]

such that

[TABLE]

and total energy can be written as

[TABLE]

The dynamical equations can be derived in much the same way as for the kinetic energy densities. Starting from equation (3), it is straightforward to show that

[TABLE]

This representation differs from the one used in (Kida and Orszag, 1990) given that we use $q$ as a building block and allows for a more straightforward interpretation of the terms. The second term on the right hand side is analogous to the corresponding term in equation (4). Although the second and third term could be combined, we will show in the following that they have different interpretations in the context of energy transfers.

The representation of the internal energy equation in Fourier space thus reads

and for the variable $\mathbf{v}=\rho^{1/3}\mathbf{u}$ introduced by Kritsuk et al. (2007):

[TABLE]

The three definitions of the turbulence energy spectrum convey different information. The spectrum of pure velocity modes is given by $E_{u}(k)\propto k^{-5/3}$ in the weakly compressible case, while $E_{u}(k)\propto k^{-2}$ is the expected scaling in shock-dominated flow at high Mach numbers. This trend is indeed observed for the exponents $\beta_{u}$ of power-law fits $E_{u}(k)\propto k^{\beta_{u}}$ in the range of wavenumbers inferred from the transfer functions (see Table 3). For example, the Kolmogorov slope $\beta_{u}\approx-1.67$ is recovered in the subsonic case ($V=0.25$) with purely solenoidal forcing and a slope $\beta_{u}\approx-1.93$ close to Burgers scaling is found for Mach numbers of about $5$ ($V=4.0$) with mainly compressive forcing ($\zeta=1/4$ and $1/8$). In Fig. 10, sample spectra are plotted for $V=1.0$ (a,b,c) and $V=4.0$ (d,e,f). In the high-Mach case, the strong impact of the forcing parameter $\zeta$ is palpable. As observed for the cross-scale fluxes, the spectrum functions for $\zeta>1/2$ are markedly different from those with a stronger solenoidal component. In the latter case the estimated inertial range shown as thin black line segments is implausible, as it falls within the bottleneck bump at high wavenumbers. It appears that the bottleneck effect becomes even more pronounced for lower Mach numbers. This might explain why the energy transfer fails to vanish at intermediate wavenumbers in some cases, such as $V=1.0$ and $\zeta=2/3$. In this case, the inferred inertial subrange shrinks to zero and no power law fits are applied (the corresponding simulations are thus omitted in Table 3).

For each spectrum, we also show the sonic wavenumber $k_{\rm s}$, which is implicitly defined by444Here we use the initial value of $c_{\rm s}$, which is close to the mean speed of sound for $\alpha=1000$.

[TABLE]

i.e. the velocity fluctuation exceeds the speed of sound for $k<k_{\rm s}$ (Schmidt et al., 2009; Federrath et al., 2010). Although the inertial-range fits for $V=1.0$ cover wavenumbers well above $k_{\rm s}\approx 4.2$, the relatively steep slopes in Table 3 reflect a significant impact of compressibility effects. The spectra plotted in Fig. 11 show that $k_{\rm s}$ is shifted with increasing $V$ toward higher wavenumbers in the inertial subrange, thus indicating the transition from moderately compressible to supersonic turbulence. Comparing to Fig. LABEL:fig:transf_alpha, however, one can see that compression-related transfers $\mathcal{T}_{\mathrm{SSc}}^{\mathrm{K}}$, $\mathcal{T}_{\mathrm{SU}}^{\mathrm{K}}$, and $\mathcal{T}_{\mathrm{US}}^{\mathrm{K}}$ have peaks at wave numbers significantly above $k_{\rm s}$. As a consequence, wavenumbers for which compression effects are relatively strong are not limited by the sonic wavenumber.

Motivated by the dimensional expression $\rho u^{3}/L=v^{3}/L$ for energy flux, it was suggested that the incompressible $k^{-5/3}$ scaling can be extended into the compressible regime for the variable $\mathbf{v}$ Kritsuk et al. (2007). The fits to $E_{v}(k)$ listed in Table 3 show that $\beta_{v}$ is close to $-5/3$ for solenoidal forcing at low to moderate Mach numbers ($V\leq 1.0$). At higher Mach numbers, it appears that $\beta_{v}=-5/3$ can be reached for increasingly compressive forcing, although the results for $V=4.0$ are not fully conclusive. The corresponding spectra are shown in Figs 10 (b,e) and 11 (b). Similar to the cross-scale energy flux discussed in the previous section, we see a transition of the shape of $E_{v}(k)$ for $V=4.0$ when the forcing becomes dominated by compressive modes. A prominent feature for $\zeta=1/4$ and $1/8$ is the shift of the kink at which the forcing peak joins into the flat part of the spectrum to wavenumbers above 10. It turns out that this is also the lower bound of the fit range (solely established on the basis of the transfer functions). For higher $\zeta$ (mainly solenodial forcing), these two points move apart and the spectra become flatter (see Fig. 10 and corresponding fit parameters in Table 3)

The spectrum function $E_{w}(k)$ corresponds to the kinetic energy density of compressible density (see Section LABEL:sec:ShellEn). It has previously been noted that $E_{w}(k)$ tends to be flatter than both $E_{u}(k)$ and $E_{v}(k)$ Kritsuk et al. (2007); Federrath et al. (2010); Grete et al. (2017), possibly approaching $k^{-4/3}$ in the strongly supersonic case. Qualitatively, the plots in Fig. 10 (c,f) indicate that $E_{w}(k)$ becomes flatter than the reference $k^{-5/3}$ spectrum for high Mach number ($V=2.0$ and $V=4.0$). As can be seen from the power-law fits listed in Table 3, $\beta_{w}$ is closer to $-4/3$ if the forcing is dominated by solenoidal modes. This trend can also be seen in the plot showing $E_{w}(k)$ for fixed $\zeta=2/3$ in Fig. 10 (c). In view of our analysis of energy transfer and cross-scale fluxes, it follows that a slope of $-4/3$ for the energy spectrum function $E_{w}(k)$ is incommensurate with an inertial-range cascade.

IV Conclusions

We investigated energy transfer functions and the resulting cross-scale fluxes of kinetic and internal energy in implicit large eddy simulations of forced compressible turbulence. To maintain a statistically stationary state even at high Mach numbers, the net increase of internal energy caused by the dissipation of kinetic energy is compensated by a simple linear cooling function. We varied both the amplitude of the forcing, resulting in Mach numbers ranging from subsonic to supersonic, and the mixture of solenoidal and compressive modes. Indicators of nearly isothermal turbulence are approximately constant RMS Mach numbers and log-normal density PDFs with significant skewness in the highly compressible regime Federrath et al. (2010); Federrath (2013); Hopkins (2013).

To compute transfer functions we apply a spectral decomposition into shells based on the variables $w=\sqrt{\rho}\,u$ (kinetic energy transfer) and $q=\sqrt{\rho}\,c_{\mathrm{s}}$ (internal energy transfer), where $u$ and $c_{\mathrm{s}}$ are the gas velocity and speed of sound, respectively. For a numerical resolution of $1024^{3}$, we find only limited regimes of vanishing energy transfer, which is a defining property of an inertial subrange. In a resolution study (see Appendix B) we find that extended regimes of vanishing transfers and, thus, an inertial subrange just begins to emerge for $1024^{3}$ grid cells. For the majority of simulations, the total transfer function crosses zero at intermediate wavenumbers and we consider shells to be close to inertial-range scaling if the modulus of the total transfer is less than $10\,\%$ of the peak value at higher wavenumbers.

For the comparison of the kinetic and internal energy budgets, cumulative quantities are better suited than the transfer functions. Cross-scale fluxes quantify the total amount of energy that is passed from all shells below a given wavenumber into shells above that wavenumber (we do not consider the question of interaction locality of energy transfer here). Our analysis shows that the cross-scale flux is always dominated by kinetic energy. This is of course expected in the weakly compressible case, where density fluctuations contribute only little to the energy transfer between different shells. Since the gas is nearly isothermal, changes in $q$ do not play a significant role. The flux of internal energy becomes more important at transonic Mach numbers, but decreases relative to the kinetic energy flux toward higher Mach numbers. This is simply a consequence of the large fluctuation of the velocity compared to the speed of sound. Moreover, the small transfers between the kinetic and internal energy reservoirs, which are related to pressure dilatation, support the theoretical reasoning in Aluie (2011). As a consequence, it appears that the kinetic energy cascade becomes asymptotically invariant in the limit of high Mach numbers.

We also find a strong sensitivity of the energy transfers on the composition of the forcing for supersonic turbulence similar to the previously reported sensitivity of density PDFs or energy spectra (Federrath et al., 2010). Particularly the advective component of the kinetic energy flux turns out to be nearly constant for intermediate wavenumbers if the compressive fraction of the force field is sufficiently large. At lower Mach numbers, roughly constant kinetic energy flux is only found for mostly solenoidal forcing. Together with the compressive component of the kinetic energy flux and the internal energy flux, the behaviour is less clear. A simple statistical argument suggests a ratio of two to one for the power of solenoidal and compressive large-scale modes (corresponding to $\zeta=1/2$; see (Federrath et al., 2010)) in the strongly compressible regime. Although our results indicate that there is a preferred mixture of solenoidal and compressive forcing, it appears that the ratio changes depending on the Mach number. This is probably a consequence of the increasingly strong coupling of solenoidal and compressive modes. To disentangle the coupling and its effect on the inertial range, further separation of the transfers into divergence-free and rotation-free modes is required.

The impact of the forcing mixture is corroborated by power-law fits to turbulent energy spectra in a range that is solely based on the transfer functions. Although the resolution is barely sufficient to resolve inertial-range scaling (for high Mach numbers, the fit range is shifted toward the bottleneck range in which the spectrum functions are tilted), we can nevertheless discern some trends. The slopes $\beta_{u}$ of the pure velocity spectra $E_{u}(k)$ show the expected dependence on the Mach number. While $\beta_{u}\approx-5/3$ for the lowest Mach number, the slope approaches $-2$ (Burgers turbulence) for higher Mach numbers (if the forcing is mainly solenoidal, the slope even falls below $-2$). We do not find a clear separation of scalings below and above the sonic wavenumber, as proposed in Galtier and Banerjee (2011). Based on our analysis, the inertial range crosses the sonic wavenumber and the slope gradually changes with increasing Mach number. However, as mentioned above, this also entails a stronger bias by the bottleneck effect. The compressible energy density spectrum $E_{w}(k)$ tends to become shallower with increasing Mach number, with $\beta_{w}$ between about $-4/3$ and $-5/3$ depending on the forcing. The dependence of the slope $\beta_{v}$ for the variable $v=\rho^{1/3}u$ on the weighing of solenoidal and compressive modes deserves closer attention. In Kritsuk et al. (2007), the Kolmogorov value of $-5/3$ was proposed for $\beta_{v}$ as a universal scaling exponent. With increasing Mach number, we find that $\beta_{v}$ tends to be closer to $-5/3$ for increasingly compressive forcing. In contrast to Federrath et al. (2010); Federrath (2013), we do not find a much steeper slope provided that the range of power-law fits is constrained by vanishing transfer functions (e.g. $\beta_{v}=-2.1$ for $\langle\mathrm{Ma}^{2}\rangle^{1/2}\approx 5.5$ and $\zeta=0$ in Federrath et al. (2010) was obtained for a wavenumber range that is more strongly affected by the forcing term $\mathcal{F}^{\mathrm{K}}$ in the shell-energy equation (LABEL:eq:RateEkin_shell) than in the range we obtain from our transfer analysis). Taken by themselves, the spectrum functions $E_{v}(k)$ do not exclude steeper fits. Galtier & Banerjee Galtier and Banerjee (2011) indeed favor a broader definition of compressible turbulence, with a scale-dependent compressible contribution to the energy flux in addition to a scale-invariant major component that can be identified with the advective component of the kinetic energy flux. Interpreted in this theoretical framework, our findings indicate that inertial-range scaling in the more specific sense of nearly constant energy flux is limited to a certain region in parameter space in which the net compressible contribution is relatively small. A direct comparison to observational data is difficult because the full three-dimensional structure of the velocity and density fields has to be reconstructed. For example, observations of star-forming clouds favor scaling exponents in between the Kolmogorov and Burgers exponents Federrath et al. (2010); Roman-Duval et al. (2011). By using indicators such as density PDFs, some evidence has been found that the driving mechanism resembles mixed forcing Kainulainen et al. (2013).

It has to be stressed that turbulence in a periodic box driven by stochastic forcing is an idealistic model. For turbulent flows occurring in nature, the actual mode of energy injection and realistic boundary conditions and maybe even initial conditions have to be taken into account. In the case of astrophysical turbulence, gravity, magnetic fields, and highly non-linear cooling increase the complexity even further. Moreover, microphysical effects such as deviations from local thermodynamic equilibrium are expected to play a role. However, this is beyond the scope of our simple model (as implied by the EOS and cooling function) and is left for more advanced studies in the future. Our numerical study demonstrates that even under the most ideal circumstances (statistically homogeneous and isotropic flow in a steady state), compressible turbulence exhibits a rich phenomenology. Apart from the incompressible limit ($q\gg w$), the notion of a compressible inertial-range cascade probably has to be regarded as asymptotic property in the limit of high Mach numbers ($q\ll w$). However, compressible turbulence in astrophysical systems is often found in the intermediate regime ($q\sim w$). An important example is the intracluster medium in clusters of galaxies Miniati and Beresnyak (2015); Schmidt et al. (2016).

Acknowledgments

We thank Ann Almgren (LBNL) for supporting us in the implementation of stochastic forcing into the Nyx code. PG acknowledges funding by NASA Astrophysics Theory Program grant #NNX15AP39G. The simulations presented in this article were performed on the Cray XC (TDS) system of the North-German Supercomputing Alliance (HLRN) facilities in Berlin. We also acknowledge the yt toolkit by Turk et al. (2011) that was used for our analysis of numerical data.

Appendix A Computational aspects

Nyx uses a hierarchical parallelisation scheme based on the AMReX framework555See https://amrex-codes.github.io/amrex. The grid is decomposed into a suitable number of boxes, which are distributed among MPI tasks. Each box in turn is split into several tile boxes Zhang et al. (2016), which extend over the full size of a box in $x$-direction (corresponding to the innermost loops in the Fortran core routines of the code) and covering quadratic tiles of size $16\times 16$ in the $yz$-plane. Each OMP thread computes one or more tile boxes. This scheme is also beneficial for efficient vectorization. For high performance, it is of particular importance to enable vectorization for the inverse Fourier transform of the forcing spectrum to spatial grid, which is basically a multiply-add operation over an array with more than 100 elements for each grid cell (corresponding to the number of non-zero Fourier modes; the spectrum of the forcing is exactly as specified in Schmidt et al. (2009)). Source terms (forcing and heating/cooling) are treated with a Strang splitting method which is described in detail in Almgren et al. (2013). The full code including the implementation of stochastic forcing is publicly available on Github.666AMReX repo: github.com/AMReX-Codes/amrex,

Nyx repo: github.com/AMReX-Astro/Nyx.

The simulations presented in this article were performed on Xeon Phi (Knights Landing) processors with AVX-512 vectorization. We used 32 nodes for the largest grids with $1024^{3}$ cells. On each node, 8 MPI tasks ran with 32 OMP threads per task on 64 physical cores (i.e. we employed the maximal number of hyperthreads, resulting in a speedup by a factor of about $2$ compared to no hyperthreading). The nodes were operated in cache mode with Sub-NUMA-2 clustering. With this configuration, the wallclock time per CFL timestep was about $7$ seconds. The total number of timesteps required for a sufficient number of dynamical time scales varied between roughly 20000 and 40000.

The postprocessing was carried out with Python, using parallel Fourier transforms implemented in the mpi4py-fft package777See mpi4py-fft.readthedocs.io. to compute transfer functions. To read data cubes from the simulations, we used the AMReX frontend implemented in yt888See yt-project.org/doc. Turk et al. (2011).

To fit power laws to energy spectra, we applied the scipy library function curve_fit to logarithmic data, i.e. the data model is $y=\beta x+y_{0}$ with $x=\log k$ and $y=\log E(k)$. Since our basis hypothesis is that the energy transfer vanishes in the inertial range, we assume the uncertainty sigma of the data in curve_fit to be proportional to the deviation of the energy transfer from zero. As a result, shells with energy transfer close to the upper bound of $0.1$ relative to the peak value have less weight in relation to the shell in which the transfer crosses zero.

Appendix B Resolution study

Figure 12 shows the time evolution of the mean internal (a) and kinetic (a) energies for $V=1.0$ and $\zeta=2/3$. For this particular case we performed simulations with $256^{3}$, $512^{3}$, and $1024^{3}$ grid cells. Within statistical variations due to the randomness of the forcing, the energies saturate at comparable levels. This implies that dissipation and cooling are insensitive to numerical resolution (see also Schmidt et al. (2006); Schmidt and Federrath (2011)).

The resolution dependence of the total transfer functions is shown in Fig. 13. In this case, we also computed the energy transfer due to the forcing, $\mathcal{F}^{\mathrm{K}}\equiv\mathcal{T}_{\mathrm{FU}}^{\mathrm{K}}$ in eq. (LABEL:eq:RateEkin_shell) (see also Grete et al. (2017)). Figure 13 (a,b,c) shows that the contribution of forcing is significant in comparison to non-linear kinetic energy transfer for wavenumbers lower than about $20$. For small wave numbers, the forcing dominates and $\mathcal{T}_{\mathrm{tot}}^{\mathrm{K}}\simeq-\mathcal{T}_{\mathrm{FU}}^{\mathrm{K}}$ (this can be seen in panel (a) for the highest resolution). The kinetic and internal energy transfer functions, $\mathcal{T}_{\mathrm{UU}}^{\mathrm{K}}$ (d,e,f) and $\mathcal{T}_{\mathrm{SS}}^{\mathrm{K}}$ (g,h,i), are qualitatively similar. However, an important feature that clearly changes with resolution is the location of the positive peaks. Particularly for the advective components, the peaks are shifted to the left as the range of numerically resolved wavenumbers decreases. This demonstrates that these peaks are coupled to the numerical dissipation scale, as explained in Section LABEL:sec:numericalTransfer. The compressive components, on the other hand, are less affected by the grid resolution (at least for resolutions $512^{3}$ and $1024^{3}$), which suggests a physical origin of the maximum. For resolutions below $1024^{3}$, the total transfer is positive down to low wavenumbers. As a consequence, there is no wavenumber interval qualifying as inertial subrange according to the criterion defined in Section LABEL:sec:numericalTransfer. This is also reflected by the transfers $\mathcal{T}_{\mathrm{SU}}^{\mathrm{K}}$ and $\mathcal{T}_{\mathrm{SU}}^{\mathrm{K}}$ (j,k,l), which are relatively weak in the inertial subrange for the highest resolution. For $512^{3}$, the minima (of the absolute values) at $k\approx 10$ are only marginal and they disappear completely for the lowest resolution. This supports the proposition of Aluie (2011) that pressure-dilation effects causing exchange of kinetic and internal energy should be weak in the inertial subrange.

The resulting energy spectrum functions are plotted in Fig. 14. The dissipative cutoff of the spectra shows the expected dependence on the numerical resolution. Even at intermediate wavenumbers, the spectra hardly converge as an inertial subrange only begins to emerge for $1024^{3}$ grids. This highlights the difficulty of applying power-law fits based on the appearance of the spectrum functions. One could certainly determine fits also for the lower resolutions. Depending on the chosen wavenumber range, this would result in significantly different slopes. It can also be seen that the bottleneck effect tends to become stronger with increasing resolution. These difficulties in establishing inertial range scalings were also reported in Kritsuk et al. (2007); Federrath (2013) for resolutions up to $4096^{3}$.