© The Authors 2026

1. Introduction

The statistical properties of supersonic compressible turbulent flows are not well predicted from first principles because the typical closure equations cannot be applied to such highly non-linear flows (Bertoglio et al. 2001; Zhou 2021). Therefore, their knowledge is mainly based on phenomenological models. For example, the 3D power spectrum (PS)¹ of the velocity field is known to be proportional to k⁻⁴ in the inertial range for supersonic compressible turbulent flows (Burgers 1948). It is interpreted as the consequence of the discontinuities created by supersonic shocks. However, the PS of the density field, though it is usually observed with a slope close to k⁻³ both in simulations of supersonic turbulence (Kim & Ryu 2005; Konstandin et al. 2016) and in the supersonic interstellar medium (Pingel et al. 2018), has not been explained so far. Simulations (Kim & Ryu 2005; Konstandin et al. 2016; Pan et al. 2022) usually show that its power-law exponent flattens as the Mach number increases. A similar evolution has also been observed in laboratory experiments (White et al. 2019). Squire & Hopkins (2017) has suggested a phenomenological model of shocks that takes into account the intermittency of the turbulent flow and that predicts the power-law exponent in the high Mach number limit. However, their model strongly depends on ill-constrained parameters. Rastegar et al. (1998) has also developed a model that neglects the pressure gradient, which could be a valid approximation in the high Mach number limit, and found a slope close to –3. None of these models predict the evolution of the power-law exponent with the Mach number.

In the astrophysical context, the PS of the density field can serve as a probe of the evolutionary stage of molecular clouds (Federrath & Klessen 2013; Burkhart et al. 2015). A shallow slope suggests that a significant amount of gas has collapsed into stars, while a steeper slope suggests that the cloud is still at the beginning of its evolution and that star formation has not started. In this paper, based on a time invariant quantity that holds in statistically homogeneous compressible turbulence (Chandrasekhar 1951; Jaupart & Chabrier 2020; Dumond et al. 2025; D25 hereafter), we derive a model that relates the slope of the density PS to the variance of the density field in the case of isothermal turbulence.

2. Model of the power spectrum of the density field

2.1. The mass invariant

Assuming ergodicity, and thus statistical homogeneity of the turbulent flow, Chandrasekhar (1951) and Jaupart & Chabrier (2021) derived a time invariant based on the continuity equation. This invariant holds if and only if the cross-correlation function, c_ρ,ρv(|q|), between the density, ρ, and the momentum, ρv, decays faster than 1/|q|² at infinity, where q denotes the spatial distance. This invariant is defined as

$$\begin{array}{c}\hfill M_{\mathrm{inv}}=\mathbb{E}(\rho )(t)\mathrm{Var}{\left(\frac{\rho }{\mathbb{E}(\rho )}\right)}_t{(l_{\mathrm{c}}^{\rho })}_t^3=\mathrm{const},\end{array}$$(1)

where $l_{\mathrm{c}}^{\rho }$ is the correlation length of the density field (Jaupart & Chabrier 2021, 2022) and $\mathbb{E}(\rho )$ and Var(ρ) are the expectation value and the variance of the density field, respectively.

This time invariance has been verified numerically in the case of decaying turbulence for various Mach numbers and equations of state in D25 (see also Appendix B). Moreover, D25 also showed that this time-invariance property implies that M_inv is Mach invariant in the supersonic regime (see discussion in D25, Sect. V.a). This Mach invariance of M_inv has been confirmed numerically for a log-density variance of ${\sigma }_s^2\lesssim 5$, corresponding to a Mach number of ℳ ≲ 10.

The invariant does not depend on the dissipation scale either because the latter will have no impact on the cross-correlation between the density and momentum at large scales. The invariant is thus not modified by variations of the dissipation scale.

To determine the value of the invariant, we relied on the model of the log-density field presented in D25 and recalled in Appendix A. Assuming that in the transonic limit (ℳ ∼ 2), the non-Gaussian features of the log-density field are small, we get

$$\begin{array}{c}\hfill M_{\mathrm{inv}}=\alpha M_{\mathrm{inj}},\end{array}$$

with α = 2 × 10⁻² and $M_{\mathrm{inj}}=\mathbb{E}(\rho )L_{\mathrm{inj}}^3$, where L_inj is the scale of maximum energy injection. This value is obtained without any fitting parameter and is consistent with the value of α measured within the 1σ uncertainties in the numerical simulations of D25. As discussed in D25, the value of α predicted by the model retains some uncertainty due to the lack of constraints on the Mach number to be used to determine M_inv from the log-density statistics. Although this will affect the prediction of α by at most 20% and the prediction of the slope of the density PS by the model by at most 5%, it will not affect the qualitative evolution of the slope with the Mach number.

2.2. Functional form of the power spectrum of the density field

By analogy with the Richardson cascade (Richardson 1922), we parametrized the PS of the density field with the following functional form:

$$\begin{array}{c}\hfill P_{\rho }^{3\mathrm{D}}(k)=\mathrm{Var}\left(\frac{\rho }{\mathbb{E}(\rho )}\right)\frac{A}{{(1+{\left(\frac{k{L}_{\mathrm{inj}}^{\max }}{2\pi }\right)}^{2p})}^{\frac{\eta }{2p}}}{\mathrm{e}}^{-k\frac{L_{\mathrm{diss}}}{2\pi }}.\end{array}$$(2)

Since this functional form is not derived from first principles, we introduced an additional parameter, p, which controls the sharpness of the transition between the injection range and the inertial range. Furthermore, this formulation involves $L_{\mathrm{inj}}^{\max }$, the largest scale at which forcing is applied to the medium. This functional form models the three ranges that characterize the PS of the density field. At scales larger than the injection scale ($2\pi /k>L_{\mathrm{inj}}^{\max }$), the density field is completely uncorrelated, and the PS is expected to be flat. Between the injection scale and the dissipation scale, the inertial range is characterized by a power law with a slope −η that we aim to determine. For scales smaller than L_diss, the PS decays exponentially in the dissipation range in analogy with the functional form of the dissipation range of the velocity power spectrum (Burgers 1948; Frisch 1995). We also introduced the normalization coefficient A, which ensures that the integral of the PS is equal to the variance of the density field (see Eq. (C.1)). In Fig. 1, we present the density power spectrum measured in our numerical simulations together with its fit using Eq. (2) with p = 2. The power spectrum was computed from the square of the Fourier transform (FT) of the fluctuations of the density field, $\rho -\mathbb{E}(\rho )$. From the Wiener-Khintchine theorem, this is equivalent to the FT of the autocovariance function for k ≠ 0. We find that this functional form provides a very good description of the measured power spectrum. In Appendix D, we show that our results depend only weakly on the chosen value of p for p ≥ 2. For smaller values of p, the transition between the injection and inertial ranges becomes less rapid and does not accurately reproduce the power spectrum obtained from the numerical simulations.

Fig. 1. Measured PS of the density field at ℳ = 7 and its fit with Eq. (2) with p = 2 for two different resolutions. The PS has been normalized by $8\mathbb{E}{(\rho )}^2{L}_{\mathrm{inj}}^3$ such that its asymptotic value at k → 0 corresponds to k3α. The dashed black line corresponds to k3α, with α = 2 × 10−2, which is the value predicted by our model (see also D25). The shaded area corresponds to the 1σ uncertainties estimated from the time variations of the density power spectra. The inset of the left panel shows the projected density field for the 5123 resolution.

2.3. Prediction of the slope of the density power spectrum

Due to the non-Gaussian features of the log-density field increasing with Mach number (Hopkins 2013; Squire & Hopkins 2017, D25) and the flattening of the log-density PS at high Mach numbers (Brucy et al. 2024) making Eq. (A.1) increasingly inaccurate, the model of the statistics of the log-density field (Appendix A) cannot be used to compute the PS of the density field directly. To bypass this limitation, we used the aforementioned invariant to predict the slope of the density PS in the supersonic regime. As shown in Appendix C, we can relate the density variance to the slope η:

$$\begin{array}{c}\hfill \mathrm{Var}\left(\frac{\rho }{\mathbb{E}(\rho )}\right)=32\pi \alpha {\displaystyle {\int }_0^{\infty }}\frac{u^2{\mathrm{e}}^{-u{\beta }^{-1}}}{{(1+{(\varphi u)}^{2p})}^{\frac{\eta }{2p}}}\mathrm{d}u,\end{array}$$(3)

where we have introduced β = L_inj/L_diss, which quantifies the width of the inertial range and $\varphi =L_{\mathrm{inj}}^{\max }/L_{\mathrm{inj}}$. Because α does not depend on the density variance nor on the dissipation scale and ϕ is imposed by the forcing properties of the system, the model shows that the slope η depends on both the density variance, $\mathrm{Var}\left(\frac{\rho }{\mathbb{E}(\rho )}\right)$, and the inertial width, β.

2.4. Comparison with numerical simulations

The data used to compare our predictions to numerical experiments are taken from the simulations presented in D25 and detailed in Appendix B. These are uniform grid simulations run with the hydrodynamical code RAMSES (Teyssier 2002). The numerical method is based on a second-order Godunov solver scheme, and the turbulence was forced using an Ornstein-Uhlenbeck process on the acceleration in Fourier space. The originality of this set of simulations is that turbulence is forced at scales significantly smaller than the box scale (between L_box/6 and L_box/8), which ensures that the statistical estimates of the correlation length, of the variance, and of M_inv are robust as shown in D25.

To determine the slope of the PS of the density field, we fit the measured PS of the density field with the functional form of Eq. (2), with η being the only free parameter. Because the dissipation range has been found to be universal in grid simulations, where the dissipation is induced by the grid (Federrath 2013), we used the best-fitting value of L_diss normalized to the spatial resolution Δx: L_diss = 5.5 Δx. This value is valid for the numerical solver used in the present study.

In Fig. 1, we plot our measured PS (solid) and the corresponding fit (dashed) for two different resolutions. This independence of α from the spatial resolution is illustrated by the fact that $\underset{k\to 0}{lim}{P}_{\rho }(k)$ does not vary within the 1σ uncertainties illustrated by the shaded area.

In Fig. 2, we plot the evolution of the slope of the PS of the density field predicted by the model (Eq. (3)) with the density variance for several dissipation lengths. We compare these predictions with the slope of the PS of the density field measured in the simulations of D25. The model predicts the evolution of the slope of the PS with the density variance very well.

Fig. 2. Predicted (solid lines) and measured (dots) slope −η in the inertial range of the PS of the density field with the variance at different resolutions, taking p = 2. The black curve corresponds to the predicted evolution of the slope in the absence of dissipation (Ldiss = 0). The gray shaded area corresponds to the region where the model does not predict the evolution of the slope because the invariance of Minv is not verified for such a small variance. The error bars correspond to the ±1σ uncertainty estimated from the time variations of the measured quantities. The shaded areas around the solid curves correspond to different predictions of the model considering an uncertainty of 20% on the determination of the invariant around its predicted value α = 2 × 10−2.

In the limit of zero viscosity (β → ∞), the model predicts that the slope of the PS of the density field tends to –3 in the high-variance limit. This result is consistent with the theoretical result of Rastegar et al. (1998) that predicts a slope of –2.95 in the pressureless limit (the sound speed tends to 0), which is equivalent to the high Mach number (or high variance) limit. In Appendix F we explain why we do not compare our prediction to other numerical studies.

This model highlights the fact that the variation of the slope of the density PS with the Mach number observed in numerical simulations (Kim & Ryu 2005; Konstandin et al. 2016) is a direct consequence of the dissipation scale imposed by the grid. In incompressible turbulence, for instance, the dissipation scale L_diss depends on the velocity dispersion σ_v. In realistic supersonic compressible turbulence, the dependence of the dissipation scale on the Mach number is unconstrained but has no reason to be constant with respect to the Mach number.

In Appendix E we show that adding an explicit viscosity, which reduces the inertial range without changing the resolution, modifies the slope of the density power spectrum similarly to grid viscosity. This confirms that the slope depends on viscosity. Furthermore, if the slope of the power spectrum was independent of viscosity, as is commonly assumed in turbulence, the slope should be steeper than –3. However, plasma experiments report density power spectra with slopes shallower than –3 (White et al. 2019). This supports the conclusion that the slope of the density power spectrum depends on viscosity.

3. Application to astrophysical systems: Turbulent molecular clouds

Here, we compare the slope of the observed PS in molecular clouds and our prediction. Such a comparison, however, is difficult because neither the variance of the volume density field nor the ratio of the dissipation scale over the injection scale are well constrained. It is possible to get an estimate of the variance from the Mach number using the first-order relation $\mathrm{Var}(\rho /\mathbb{E}(\rho ))=b^2{\mathcal{M}}^2$ (see, e.g., Konstandin et al. 2016). Taking b = 0.5, which corresponds to an equipartition forcing between solenoidal and compressive modes (Federrath et al. 2010), and a typical Mach number ranging between 4 and 20 for star-forming clouds (Larson 1981; Chabrier et al. 2014) leads to a variance ranging between 4 and 100 (Hennebelle & Falgarone 2012).

Moreover, the processes that dissipate turbulence in molecular clouds are not well identified. Although it has been suggested that the main dissipation process of turbulence at small scales is ambipolar diffusion (see, e.g., Momferratos et al. 2014), this view has been challenged recently (Pineda et al. 2024). The fact that the power law of the density PS is observed over more than three decades in spatial scales means that the inertial width is β > 10³ for such systems.

In Fig. 3, we plot β against the Mach number for various slopes of the density PS. For the above-mentioned typical star-forming conditions, we expected the slope of the PS to vary from –2.6 to –3.2. This is perfectly consistent with most observations of molecular clouds (Brunt 2010; Miville-Deschênes et al. 2010; Hennebelle & Falgarone 2012; Miville-Deschênes et al. 2016; Pingel et al. 2018; Pineda et al. 2024). Shallower slopes are also observed in clouds that are actively undergoing star formation (see Federrath & Klessen 2013 for a summary). Such a behavior is also predicted from our model since the variance in such a medium is expected to be much larger than the one that would be expected in a purely turbulent medium because the self-gravity of the medium creates regions with a high density contrast.

Fig. 3. Inertial width, β, against the Mach number for various slopes, η, of the inertial range of the PS of the density field. Notably, β has been observed to be larger than 103 in the interstellar medium (Miville-Deschênes et al. 2016; Pineda et al. 2024), and thus we shade the excluded region in gray. The shaded areas around the solid curves represent the same as in Fig. 2.

Even when assuming a relation between the density variance and the Mach number, this study illustrates the fact that it is not possible to directly infer the Mach number from the slope of the power-law density field for η < 3 (see Miville-Deschênes et al. 2010, 2016; Pineda et al. 2024) because it depends on the unknown width of the inertial range. For η > 3, comparing the observed slope with the one measured in the simulations still gives a good estimate of the Mach number because the latter depends weakly on the dissipation scale.

Our model gives a criterion on the inertial width (i.e., the resolution) that is needed to simulate a turbulent flow characterized by a given density variance and a PS with a given slope. From Fig. 3, if one wants to simulate a molecular cloud at Mach 5 with a realistic slope η = 3 (Miville-Deschênes et al. 2016), one needs a resolution allowing β ≃ 200. Starting with realistic density field statistics is of prime importance because such turbulent flows are often used as initial conditions before activating gravity to study star formation in the turbulent interstellar medium (Federrath & Klessen 2013; Brucy et al. 2024). Starting with biased statistics of the density field will impact the statistics of the collapsed structures obtained at the end of the simulation. For example, a PS of the density field with a shallower slope will increase the density of small scale structures, making their collapse easier. Nevertheless, the precise quantification of this effect is beyond the scope of this paper.

4. Discussion and conclusion

We have presented a model of the slope of the PS of the density field in statistically homogeneous isotropic isothermal turbulence in numerical simulations. We demonstrated that this slope depends on both the variance of the density field and the width of the dissipation range. Unlike the conventional view that the inertial range is unaffected by small-scale viscous processes, our results suggest the opposite. Specifically, the density power spectrum does not exhibit an inertial range in the classical turbulent sense, in contrast to the velocity power spectrum, whose inertial range remains independent of viscosity. Mathematically, this distinction arises because the Chandrasekhar invariant constrains the density power spectrum at zero wavenumber to be constant in time, even in the presence of external forcing and viscosity. No analogous Chandrasekhar-type invariant exists for the velocity field when the flow is forced or viscous. These findings indicate that viscosity influences density structures across all scales in compressible turbulence. However, a clear physical interpretation of this mechanism remains to be established.

Thanks to this model, we were able to reproduce the measurements of the slope of the PS of the density field in turbulent simulations very well and quantitatively explain the observed trend for the flattening slope as the Mach number increases. The model explains the slope of the PS of the density field observed in turbulent molecular clouds. It shows that it is not possible to determine the Mach number of a molecular cloud from the observation of this slope without additional constraints on the inertial width, notably when the slope is shallower than –3. Finally, we provide a criterion for the numerical resolution that must be reached to simulate a turbulent flow characterized by a given variance and slope of the PS of the density field.

Acknowledgments

We thank the anonymous referee for constructive feedback that significantly improved the quality of this work. We thank Elliot Lynch and Guillaume Laibe for very helpful discussions that improved the clarity of the manuscript. PD, JF and NB were supported by the French national research agency grant ANR-23-CE31-0005 (BRIDGES). NB cheerfully thanks Edwin Santiago Leandro for his essential contribution to the 2D version of the viscosity implementation in Ramses.

References

Appendix A: Model for the statistics of the log-density field

In this appendix we present the statistical description of the log-density field that is used to predict the value of the invariant. Similarly to D25, the statistics of the log-density random field $s=ln(\rho /\mathbb{E}(\rho ))$ is assumed to be characterized by a PS with the following form for low Mach numbers (ℳ ≲ 5):

$$\begin{array}{c}\hfill P_s^{3\mathrm{D}}(k)=\text{Var}(ln\left(\frac{\rho }{\mathbb{E}(\rho )}\right))\frac{B{\mathrm{e}}^{-k\frac{L_{\mathrm{diss}}}{2\pi }}}{{(1+{\left(\frac{k{L}_{\mathrm{inj}}}{2\pi }\right)}^2)}^2}.\end{array}$$(A.1)

Here, B is a coefficient that ensures that the integral of the PS is equal to the variance of the log-density field and $\mathbb{E}(\rho )$ is the mean density. We also introduced the dissipation scale L_diss and the length L_inj at which most of the energy is injected in the medium².

This functional form has been justified numerically in Beresnyak et al. (2005), Brucy et al. (2024) and D25. We found in our simulations that the slope of the inertial range of the PS of the log-density field is close to -3.8. The slight difference between the present parametrization (Eq. A.1) and the measured slope affects the subsequent predictions of the slope of the PS of the density field by at most 5%. A theoretical justification for an inertial range close to -4 for the log-density PS is also under investigation (Dumond et al., in prep.).

Assuming that in the transonic limit (ℳ ∼ 2), the non-gaussian features of the log-density field are limited, the statistical properties of the field are completely determined by this power spectrum. We get M_inv = αM_inj with α = 2 × 10⁻² and $M_{\mathrm{inj}}=\mathbb{E}(\rho )L_{\mathrm{inj}}^3$, where L_inj is the scale of maximum energy injection (D25).

Appendix B: Turbulent driving and hydrodynamical solver

The numerical simulations used for this study are performed with the hydrodynamical code RAMSES (Teyssier 2002) without adaptive mesh refinement (fixed cartersian grid). The boundary conditions are periodic. The numerical method is based on a second-order Godunov solver scheme. The solver is the Harten–Lax–van Leer-Contact (HLLC) approximate Riemann solver (Toro 2009). The initial conditions consist of a fluid of atomic hydrogen at rest. The fluid follows an isothermal or polytropic equation of state. The turbulence is forced using the Ornstein-Uhlenbeck forcing on the acceleration (Eswaran & Pope 1988; Schmidt et al. 2009). The equations of conservation of mass and momentum that are solved are the following:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathbf{v})& =0,\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill \rho (\frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}·\mathbf{\nabla })\mathbf{v})& =-c_{\mathrm{s}}^2\mathbf{\nabla }\rho +\rho \mathit{f},\hfill \end{array}$$(B.2)

where ρ is the density, v the velocity, c_s the sound speed and f the turbulent driving force. The viscosity is not modeled explicitly, but acts implicitly on the flow through the numerical grid viscosity. The energy equation is not solved explicitly. It is replaced by a simple polytropic equation of state P ∝ ρ^γ with γ the polytropic index. In this study, we consider only the case γ = 1 (isothermal process). The Fourier modes $\widehat{\mathit{f}}(\mathit{k},t)$ of the turbulence driving acceleration field f follow the following stochastic differential equation:

$$\begin{array}{c}\hfill \mathrm{d}\widehat{\mathit{f}}(\mathit{k},t)=-\widehat{\mathit{f}}(\mathit{k},t)\frac{\mathrm{d}t}{T_{\mathrm{OU}}}+F_0(\mathit{k}){\mathit{P}}^{\zeta }(\mathit{k})\mathrm{d}{\mathit{W}}_t.\end{array}$$(B.3)

In this equation, dt is the integration time step and T_OU is the autocorrelation timescale. As usually done in such numerical simulations (Schmidt et al. 2009), we set T_OU to the turbulent crossing time $T_{\mathrm{cross}}=L_{\mathrm{inj}}/\sqrt{⟨{\mathrm{v}}^2⟩}$, where L_inj is the injection length of turbulence and $\sqrt{⟨{\mathrm{v}}^2⟩}$ the turbulent velocity dispersion. Scannapieco et al. (2026) recently showed that the variance of the density field depends on the ratio T_cross/T_OU. However, thoughout this study, we consider only the usual case T_OU = T_cross letting the other cases for a future study. The weighting function of the driving modes F₀ allows the turbulence to be driven only within a precise range of spatial scales. In our work, we inject the turbulence isotropically in most runs (unless stated otherwise) between L_box/6 and L_box/8:

$$\begin{array}{c}\hfill F_0(k)=\{\begin{array}{c}1-{(\frac{L_{\mathrm{box}}|\mathit{k}|}{2\pi }-7)}^2\text{if}6<\frac{L_{\mathrm{box}}|\mathit{k}|}{2\pi }<8\hfill \\ 0\text{if not.}\hfill \end{array}\end{array}$$(B.4)

In most of the turbulent simulations in the literature, turbulence is injected at larger scale, typically L_box/2. Here, we chose to inject at smaller scale to ensure that the size of the simulation box is large enough compared to the correlation length (Jaupart & Chabrier 2022). Assuming ergodicity, this ensures that spatial averages are good estimate of the statistical average.

The projection operator P^ζ is a weighted sum of the components of the Helmholtz decomposition of compressive versus solenoidal modes (Federrath et al. 2010):

$$\begin{array}{cc}\hfill P_{\mathit{ij}}^{\zeta }(k)& =\zeta P_{\mathit{ij}}^{\perp }(k)+(1-\zeta )P_{\mathit{ij}}^{\Vert }(k)\hfill \\ \hfill & =\zeta {\delta }_{\mathit{ij}}+(1-2\zeta )\frac{k_i{k}_j}{{|k|}^2},\hfill \end{array}$$(B.5)

where δ_ij is the Kronecker symbol, and $P_{\mathit{ij}}^{\perp }={\delta }_{\mathit{ij}}-k_i{k}_j/k^2$ and $P_{\mathit{ij}}^{\Vert }=k_i{k}_j/k^2$ are the fully solenoidal and fully compressive projection operators, respectively. The projection operator is used to construct a purely solenoidal force field by setting the solenoidal fraction ζ = 1, which is used to drive turbulence in an incompressible system (Eswaran & Pope 1988). In the following we have chosen ζ = 0.5, which corresponds to the energy equipartition of the velocity field: 1/3 compressive and 2/3 solenoidal.

Appendix C: Relation between the slope and the variance

We derive the relation given in Eq. 3 between the slope of the PS of the density field, −η, the variance $\text{Var}\left(\frac{\rho }{\mathbb{E}(\rho )}\right)$ and the inertial width β. Starting with the functional form of the PS given in Eq. 2, its normalization coefficient A verifies

$$\begin{array}{c}\hfill \frac{1}{{(2\pi )}^3}{\displaystyle {\int }_0^{\infty }}{P}_{\rho }^{3\mathrm{D}}(k){\mathrm{d}}^{3}k=\text{Var}(\rho ),\end{array}$$(C.1)

based on the definition of the PS that is the Fourier transform of the auto-covariance function of the density field. Thus,

$$\begin{array}{c}\hfill \frac{1}{{(2\pi )}^3}{\displaystyle {\int }_0^{\infty }}A\frac{4\pi }{{(1+{(\frac{k{L}_{\mathrm{inj}}^{\max }}{2\pi })}^{2p})}^{\frac{\eta }{2p}}}{\mathrm{e}}^{-k\frac{L_{\mathrm{diss}}}{2\pi }}{k}^2\mathrm{d}k=1.\end{array}$$(C.2)

By doing the variable change u = kL_inj/(2π), we get

$$\begin{array}{c}\hfill A=\frac{L_{\mathrm{inj}}^3}{4\pi }{\left({\displaystyle {\int }_0^{\infty }}\frac{u^2{\mathrm{e}}^{-u{\beta }^{-1}}}{{(1+{(\varphi u)}^{2p})}^{\frac{\eta }{2p}}}\mathrm{d}u\right)}^{-1},\end{array}$$(C.3)

where $\beta =\frac{L_{\mathrm{inj}}}{L_{\mathrm{diss}}}$ and $\varphi =L_{\mathrm{inj}}^{\max }/L_{\mathrm{inj}}$. In this study, ϕ = 7/6.

The correlation length $l_{\mathrm{c}}^{\rho }$ is defined as

$$\begin{array}{c}\hfill l_{\mathrm{c}}^{\rho }={\left(\frac{1}{8C_{\rho }(0)}{\displaystyle {\int }_{{\mathbb{R}}^3}}{C}_{\rho }(\mathit{q}){\mathrm{d}}^3\mathit{q}\right)}^{1/3},\end{array}$$(C.4)

where C_ρ is the auto-covariance function of the density field. It can be related to the PS by

$$\begin{array}{cc}\hfill {\mathrm{ł}}_{\mathrm{c}}^3& =\frac{1}{2^3\text{Var}(\rho )}{P}_{\rho }^{3\mathrm{D}}(0),\hfill \end{array}$$(C.5)

$$\begin{array}{cc}& =\frac{A}{8},\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& =\frac{1}{32\pi }\frac{L_{\mathrm{inj}}^3}{{\displaystyle {\int }_0^{\infty }}\frac{u^2{\mathrm{e}}^{-u{\beta }^{-1}}}{{(1+{(\varphi u)}^{2p})}^{\frac{\eta }{2p}}}\mathrm{d}u}.\hfill \end{array}$$(C.7)

Using the value of the invariant $M_{\mathrm{inv}}=\overline{\rho }\mathrm{Var}(\rho /\overline{\rho })l_{\mathrm{c}}^3$, we can relate the correlation length to the density variance. Introducing the constant α = M_inv/M_inj where $M_{\mathrm{inj}}=\mathbb{E}(\rho )L_{\mathrm{inj}}^3$, we have

$$\begin{array}{c}\hfill \text{Var}\left(\frac{\rho }{\mathbb{E}(\rho )}\right)=\alpha {\left(\frac{l_{\mathrm{c}}}{L_{\mathrm{inj}}}\right)}^{-3}.\end{array}$$(C.8)

Together with Eq. C.7, we finally find Eq. 3:

$$\begin{array}{c}\hfill \text{Var}\left(\frac{\rho }{\mathbb{E}(\rho )}\right)=32\pi \alpha {\displaystyle {\int }_0^{\infty }}\frac{u^2{\mathrm{e}}^{-u{\beta }^{-1}}}{{(1+{(\varphi u)}^{2p})}^{\frac{\eta }{2p}}}\mathrm{d}u.\end{array}$$(C.9)

Fig. C.1. Same as Fig. 2 but with the functional form of the PS of the density field given in Eq. 2 with p = 5. The shaded areas around the solid curves corresponds to different predictions of the model considering an uncertainty of 20% on the determination of the invariant around its predicted value α = 2 × 10−2.

Appendix D: Dependence on the functional form of the density power spectrum

In this appendix we test how the prediction of our model depends on the parameter p, which control the sharpness of the transition between the injection range and the inertial range.

In Fig. C.1 we plot the measurements of the slope in the simulations and the prediction of the model obtained with p = 5 in Eq. 2. Because the functional form used here has changed, the values of the slopes measured in the simulation have slightly changed compared to Fig. 2. The agreement is still good, although the discrepancy between the prediction and the measurement is slightly larger at low variance. This is due to the fact that the functional form with p = 5 reproduces slightly less accurately the measured power spectrum at low Mach number.

Appendix E: Simulations with explicit viscosity

Since this paper suggests that the slope of the inertial range depends on the ratio between the injection scale and the dissipation scale of the density power spectrum, it is important to determine whether the dissipation scale obtained in our simulations has a physical meaning, i.e., whether it is set by the viscosity of the flow (which may be numerical if no explicit viscosity is included), or whether it is purely numerical, i.e., biased and imposed by the grid discretization.

To address this issue, we run two additional 2048³ simulations at ℳ = 5 with two different values of the viscosity. In order to allow comparison with the literature on viscous compressible turbulence, we use a similar model of shear viscosity to most of the studies of compressible turbulence with viscosity (e.g., Shivakumar & Federrath 2025). The viscous stress tensor can be thus written as

$$\begin{array}{c}\hfill T_{\mathrm{vis}}^{\mathit{ij}}=\nu \rho (\frac{\partial v^i}{\partial x^j}+\frac{\partial v^j}{\partial x^i}-\frac{2}{3}\frac{\partial v^k}{\partial x^k}{\delta }^{\mathit{ij}}),\end{array}$$(E.1)

where ν is the kinematic viscosity and δ^ij the Kroneker symbol.

As shown in Fig. E.1, increasing the explicit viscosity reduces the width of the inertial range. Moreover, it modifies the slope of the density power spectrum in a manner similar to grid viscosity, within 1σ uncertainties. We therefore conclude that the effect of grid viscosity on the density power spectrum is comparable to that of an explicit physical viscosity.

Fig. E.1. Comparison between the density power spectra at resolution 10243 (left) and 5123 (right) with no explicit viscosity, and 20483 simulations with explicit viscosity such that the dissipation length of the compressible modes lη is equal respectively to Lbox/1024 and Lbox/512. The power spectrum of the 20483 simulation without explicit viscosity is shown for reference.

To relate β to the Reynolds number, we link the dissipation length l_η to the size of the cells Δx. Since the shear viscosity dissipate primarily the solenoidal modes of the velocity field, (Beattie et al. 2025), the dissipation length, defined as the scale at which the Reynolds number equals unity, is given by

$$\begin{array}{c}\hfill l_{\eta }=L_{\mathrm{inj}}/{\mathrm{Re}}^{3/4},\end{array}$$(E.2)

where Re = ℳc_sL_injν⁻¹. Here, ℳ is the Mach number, c_s the sound speed, L_inj the injection scale, and ν the viscosity.

We verify that in a simulation with spatial resolution Δx = L_box/N (with N³ the number of cells), choosing the viscosity such that the dissipation length l_η equals to half the cell size Δx′ = L_box/N′ of a lower-resolution simulation (N′³ cells, with N′ < N) reproduces the dissipation scale of that simulation, as shown in Fig. E1. We conclude that in a simulation with no explicit viscosity, the dissipation length of the compressible modes is set to half the size of the cells by the viscosity of the grid, i.e., l_η ≃ Δx/2.

We thus obtained the following relation between β and the Reynolds number in our numerical simulation:

$$\begin{array}{c}\hfill \mathrm{Re}={(11\beta )}^{4/3}={\left(\frac{2L_{\mathrm{inj}}}{\mathrm{\Delta }x}\right)}^{4/3},\end{array}$$(E.3)

defining β = L_inj/L_diss and L_diss ≃ 11l_η ≃ 5.5Δx in our numerical setup. This relation is used in Fig. 1 and Fig. 2 to link β to the Reynolds number.

Appendix F: Absence of comparison to other studies

In this appendix we explain why we restrict our quantitative comparison to our own numerical simulations. Determining the slope of P_ρ requires fitting the measured density PS with the functional form given in Eq. 2. An accurate fit necessitates resolving the scales larger than the injection scale, as well as a precise determination of the density variance. These conditions are not satisfied in existing simulations from the literature (Kim & Ryu 2005; Konstandin et al. 2016; Pan et al. 2022) because their turbulent forcing as done at half the box length. Moreover, when forcing at such scales, the statistical estimates of the variance and the correlation length are inaccurate (D25). Nevertheless, our results remain qualitatively consistent with the trends reported in those studies.

All Figures

Fig. 1. Measured PS of the density field at ℳ = 7 and its fit with Eq. (2) with p = 2 for two different resolutions. The PS has been normalized by $8\mathbb{E}{(\rho )}^2{L}_{\mathrm{inj}}^3$ such that its asymptotic value at k → 0 corresponds to k3α. The dashed black line corresponds to k3α, with α = 2 × 10−2, which is the value predicted by our model (see also D25). The shaded area corresponds to the 1σ uncertainties estimated from the time variations of the density power spectra. The inset of the left panel shows the projected density field for the 5123 resolution.

In the text

Fig. 2. Predicted (solid lines) and measured (dots) slope −η in the inertial range of the PS of the density field with the variance at different resolutions, taking p = 2. The black curve corresponds to the predicted evolution of the slope in the absence of dissipation (Ldiss = 0). The gray shaded area corresponds to the region where the model does not predict the evolution of the slope because the invariance of Minv is not verified for such a small variance. The error bars correspond to the ±1σ uncertainty estimated from the time variations of the measured quantities. The shaded areas around the solid curves correspond to different predictions of the model considering an uncertainty of 20% on the determination of the invariant around its predicted value α = 2 × 10−2.

In the text

	Fig. 3. Inertial width, β, against the Mach number for various slopes, η, of the inertial range of the PS of the density field. Notably, β has been observed to be larger than 103 in the interstellar medium (Miville-Deschênes et al. 2016; Pineda et al. 2024), and thus we shade the excluded region in gray. The shaded areas around the solid curves represent the same as in Fig. 2.
In the text

	Fig. C.1. Same as Fig. 2 but with the functional form of the PS of the density field given in Eq. 2 with p = 5. The shaded areas around the solid curves corresponds to different predictions of the model considering an uncertainty of 20% on the determination of the invariant around its predicted value α = 2 × 10−2.
In the text

	Fig. E.1. Comparison between the density power spectra at resolution 10243 (left) and 5123 (right) with no explicit viscosity, and 20483 simulations with explicit viscosity such that the dissipation length of the compressible modes lη is equal respectively to Lbox/1024 and Lbox/512. The power spectrum of the 20483 simulation without explicit viscosity is shown for reference.
In the text