# Mass invariant in a compressible turbulent medium

Pierre Dumond, Jérémy Fensch, Gilles Chabrier, Etienne Jaupart

► To cite this version:

Pierre Dumond, Jérémy Fensch, Gilles Chabrier, Etienne Jaupart. Mass invariant in a compressible turbulent medium. Physical Review Research, 2025, 7 (3), pp.033035. 10.1103/PhysRevResearch.7.033035. hal-05154806

HAL Id: hal-05154806

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

Submitted on 16 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Mass invariant in a compressible turbulent medium

Pierre Dumond ◦, 1, ∗ Jérémy Fensch ◦, 1 Gilles Chabrier ◦, 1, 2 and Etienne Jaupart ◦, 3

¹CRAL, Ecole Normale Supérieure de Lyon, Université de Lyon, UMR CNRS 5574, F-69364 Lyon Cedex 07, France

²School of Physics, University of Exeter, Exeter EX4 4QL, United Kingdom

³CEA-DAM-DIF, F-91297 Arpajon, France

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

(Received 20 December 2024; accepted 6 April 2025; published 8 July 2025)

Predicting the measurable statistical properties of density fluctuations in a supersonic compressible turbulent flow is a major challenge in physics. In 1951, Chandrasekhar derived an invariant under the assumption of the statistical homogeneity and isotropy of the turbulent density field and stationarity of the background density. Recently, Jaupart and Chabrier [Astrophys. J. Lett. 922, L36 (2021)] extended this invariant to nonisotropic flows in a time-evolving background and showed that it has the dimension of a mass. This invariant $M_{inv}$ is defined by $M_{inv} = \mathbb{E}(\rho)\text{Var}(\frac{\rho}{\mathbb{E}(\rho)})(l_c^\rho)^3$, where $\rho$ is the density field and $l_c^\rho$ is the correlation length. In this article, we perform numerical simulations of homogeneous and isotropic compressible turbulence to test the validity of this invariant in a medium subject to isotropic decaying turbulence. We study several input configurations, namely different Mach numbers, injection lengths of turbulence, and equations of state. We confirm that $M_{inv}$ remains constant during the decaying phase of turbulence. Furthermore, we develop a theoretical model of the density field statistics which predicts without any free parameter the evolution of the correlation length with the variance of the log-density field beyond the assumption of the Gaussian field for the log density. Noting that $M_{inv}$ is independent of the Mach number, we show that this invariant can be used to relate the non-Gaussian evolution of the log-density probability distribution function to its variance with no free parameters.

DOI: 10.1103/PhysRevResearch.7.033035

I. INTRODUCTION

Properly understanding the statistical properties of com-
pressible turbulence remains a major problem in physics.
Most attempts have been limited to the study of weakly
compressible turbulence, i.e., three-dimensional (3D) Mach
number $M = \sqrt{\langle v^2 \rangle}/c_s \le 1$, where $\sqrt{\langle v^2 \rangle}$ is the 3D velocity
dispersion and $c_s$ is the sound speed, both on the theoreti-
cal side (see [1-3] for a review) and on the numerical side
(e.g., Refs. [4,5]). Little attention has been paid to supersonic
compressible turbulence, mainly because in this case none of
the usual closure equations (e.g., eddy damped quasinormal
Markovian or direct interaction approximation) can be used to
derive the statistical properties of interest for the flow. There-
fore, our current understanding of compressible turbulence is
based essentially on phenomenology, notably the description
of the field of density fluctuations as being the result of the
interaction of a large number of random shocks, allowing the
use of the central limit theorem for multiplicative independent
random processes [6,7]. However, there is no exact result to
describe the commonly measured statistics of the flow (prob-
ability density function [PDF], power spectra). Some exact

relations do exist [8], but it is difficult to extract from them
a prediction of the statistical evolution of the flow properties.

Assuming ergodicity, stationarity of the background den-
sity and statistical homogeneity [9] and isotropy of the
turbulent flow, Ref. [10] derived a temporal invariant based
on the continuity equation. He also had to assume that the
cross-correlation C<sub>ρ,ρv</sub>(|q|) between the density ρ and the
momentum ρv decays faster than 1/|q|<sup>2</sup> at infinity, where
q denotes the spatial distance. Reference [11] has recently
extended this invariant to nonisotropic turbulence and evolv-
ing statistical mean density of the flow E(ρ), yielding the
following form:

$$
M_{\text{inv}} = \mathbb{E}(\rho)(t)\text{Var}\left(\frac{\rho}{\mathbb{E}(\rho)}\right)_t (l_c^\rho)_t^3 = \text{const}, \quad (1)
$$

where $l_c^\rho$ is the correlation length, defined as

$$
l_c^\rho = \left( \frac{1}{8C_\rho(0)} \int_{\mathbb{R}^3} C_\rho(q) d^3 q \right)^{1/3}, \quad (2)
$$

where $C_\rho$ is the autocovariance function (ACF) of the density
field. We recall that, by definition, $C_\rho(0)$ is the variance of the
density field. Note that the basic mathematical expressions of
the invariant derived by [10] and by [11] are the same. The
only difference between the two is the context within which
they are applied. The invariant derived by Ref. [10] applies to
homegeneous isotropic turbulence with constant background
average density, while the one derived by Ref. [11] applies
to homogeneous nonisotropic background with time-evolving
background average density.

<footnote>*Contact author: pierre.dumond@ens-lyon.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.</footnote>
In this paper we mainly verify the invariant in the context studied by [10], e.g., for isotropic turbulence with no time variation of the statistical mean density of the flow and no gravity. We will briefly consider the case of a time-varying background density also and leave the study of the invariant in anisotropic turbulent flows and in the presence of gravity for a future study.

In Sec. II, we first present the setup we use to numerically verify this invariant, making sure that the hypotheses it relies on are verified in the simulated turbulent flow. We show in Sec. III that this quantity is indeed constant over time in decaying turbulence. In Sec. IV we derive an analytical model that describes well the statistical properties of the density field and allows one to predict the variation of the correlation length with the Mach number without any fitting parameters. Finally, in Sec. V, we show that this invariant can be used to determine analytically the evolution with the Mach number of the parameter T introduced by Ref. [12] to describe the non-Gaussianity of the PDF of the logarithm of the density field (called log-density field in the following).

## II. SIMULATIONS

### A. Numerical setup

We use the hydrodynamical code RAMSES [13] without adaptive mesh refinement (fixed Cartesian 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 [14]. 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 [15–17]. The equations of conservation of mass and momentum that are solved are the following:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \quad (3) $$

$$ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right) = -c_s^2 \nabla \rho + \rho \mathbf{f}. \quad (4) $$

The viscosity is not modeled explicitly, but acts implicitly on the flow through the numerical grid viscosity. The equation energy is not solved explicitly. It is replaced by a simple polytropic equation of state $P \propto \rho^\gamma$ with $\gamma$ the polytropic index. In this study, we will consider two cases: $\gamma = 1$ (isothermal process) and $\gamma = 5/3$ (isentropic process for monoatomic gas). The Fourier modes $\hat{f}(k, t)$ of the turbulence driving acceleration field $f$ follow the following stochastic differential equation:

$$ d\hat{f}(k, t) = -\hat{f}(k, t) \frac{dt}{T_{\text{OU}}} + F_0(k) \mathbf{P}^{\zeta}(k) d\mathbf{W}_t. \quad (5) $$

In this equation, $dt$ is the integration time step and $T_{\text{OU}}$ is the autocorrelation timescale. As usually done in such numerical simulations (e.g., Ref. [16]), we set $T_{\text{OU}}$ to the turbulent crossing time $T_{\text{cross}} = L_{\text{inj}}/\sqrt{\langle v^2 \rangle}$, where $L_{\text{inj}}$ is the injection length of turbulence and $\sqrt{\langle v^2 \rangle}$ the turbulent velocity dispersion. The weighting function of the driving modes $F_0$ 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_{\text{box}}/6$ and $L_{\text{box}}/8$:

$$ F_0(k) = \begin{cases} 1 - \left( \frac{L_{\text{box}}|k|}{2\pi} - 7 \right)^2 & \text{if } 6 < \frac{L_{\text{box}}|k|}{2\pi} < 8, \\ 0 & \text{if not.} \end{cases} \quad (6) $$

In most of the turbulent simulations in the literature, turbulence is injected at larger scale, typically $L_{\text{box}}/2$ (see discussion in Sec. II B). 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 (see Sec. II B). Assuming ergodicity, this ensures that spatial averages are good estimates of the statistical average (see discussion in Sec. III).

The projection operator $\mathbf{P}^{\zeta}$ is a weighted sum of the components of the Helmholtz decomposition of compressive versus solenoidal modes [18]:

$$ \begin{aligned} P_{ij}^{\zeta}(k) &= \zeta P_{ij}^{\perp}(k) + (1 - \zeta) P_{ij}^{\parallel}(k) \\ &= \zeta \delta_{ij} + (1 - 2\zeta) \frac{k_i k_j}{|k|^2}, \end{aligned} \quad (7) $$

where $\delta_{ij}$ is the Kronecker symbol and $P_{ij}^{\perp} = \delta_{ij} - k_i k_j / k^2$ and $P_{ij}^{\parallel} = 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 $\zeta = 1$, which is used to drive turbulence in an incompressible system [15]. In the following we have chosen $\zeta = 0.5$, which corresponds to the energy equipartition of the velocity field: 1/3 compressive and 2/3 solenoidal.

We carried out simulations at different large scale Mach numbers: $M = 2.5, 3.5, 5, 7$, and 10 to test the invariant under different conditions. The projected density of the simulation at $M = 3.5$ is shown in Fig. 1. We also consider several spatial resolutions and different injection lengths to test the convergence of the results. The list of the various simulations is given in Table I. We also give the measured $\sigma_s^2$, which is a more primitive quantity of interest, and which depends on Mach and the type of forcing and equation of state.

To test the temporal invariance of the quantity $M_{\text{inv}}$ [see Eq. (1)] suggested by Refs. [10,11], we first drive the turbulence for at least two turbulent crossing times $T_{\text{cross}}$ before stopping the driving. This ensures that the turbulence is fully developed when the driving is stopped [19]. During the turbulence decay phase, we measure $M_{\text{inv}}$ and verify whether it remains constant or not.

### B. Subboxes

Since the boundary conditions of the simulation domain are periodic, the periodicity could create artificial large scale correlations which could prevent the decay of the correlation between $\rho$ and $\rho v$, which is one of the main assumptions for the quantity $M_{\text{inv}}$ to be time invariant. Furthermore, it can be shown that the correlation length measured in a periodic box with a periodic estimator is 0 [20]. This suggests that any estimate of the correlation length in a periodic box will be inaccurate. To overcome these problems, we study only a subsystem of the whole simulation domain. The size of the subbox should be (i) small enough compared with the whole
<figure><img src="image_3.png" /><figcaption>FIG. 1. Projected density Σ of the simulation with the turbulence driven at Mach 3.5 and injected at $L_{box}/7$. The region within the white square is the one in which the study of the invariant is preformed in order to avoid the effect of the boundary conditions (see Sec. II B). The horizontal double arrow shows the size of the injection scale $L_{inj}$.</figcaption></figure>

box to be unaffected by the boundary conditions and (ii) large enough compared with the correlation length to resolve it properly. We choose a subbox with a volume of 1/8 of the whole box. We verified that changing the volume from 1/8 to 1/64 of the whole box does not affect qualitatively the results. All further analysis will be carried out in such a system.

<table><caption>TABLE I. List of simulations. From left to right: the 3D Mach number $M$, the polytropic index $\gamma$, the injection length $L_{inj}$, the number of cells $N_{cells}$, and the variance of the log-density field $\sigma_s^2$. The * denotes the fiducial run used in this study.</caption><thead><tr><th>M</th><th>γ</th><th>L<sub>inj</sub></th><th>N<sub>cells</sub></th><th>σ<sub>s</sub><sup>2</sup></th><th>Fiducial</th></tr></thead><tbody><tr><td rowspan="3">2.5</td><td rowspan="3">1</td><td>L<sub>box</sub>/7</td><td>1024<sup>3</sup><br>2048<sup>3</sup></td><td>1.37<br>1.34</td><td></td></tr><tr><td>L<sub>box</sub>/14</td><td>1024<sup>3</sup><br>2048<sup>3</sup></td><td>1.37<br>1.36</td><td></td></tr><tr><td>L<sub>box</sub>/7</td><td>512<sup>3</sup><br>1024<sup>3</sup><br>2048<sup>3</sup></td><td>2.17<br>2.07<br>2.04</td><td>*</td></tr><tr><td rowspan="3">3.5</td><td rowspan="3">1</td><td>L<sub>box</sub>/7</td><td>1024<sup>3</sup><br>2048<sup>3</sup></td><td>2.01<br>2.01</td><td></td></tr><tr><td>L<sub>box</sub>/14</td><td>512<sup>3</sup><br>1024<sup>3</sup><br>2048<sup>3</sup></td><td>3.23<br>2.94<br>2.83</td><td>*</td></tr><tr><td>L<sub>box</sub>/7</td><td>512<sup>3</sup><br>1024<sup>3</sup><br>2048<sup>3</sup></td><td>4.39<br>4.19<br>4.01</td><td>*</td></tr><tr><td rowspan="3">5</td><td rowspan="3">1</td><td>L<sub>box</sub>/7</td><td>512<sup>3</sup><br>1024<sup>3</sup><br>2048<sup>3</sup></td><td>4.39<br>4.39<br>4.01</td><td></td></tr><tr><td>L<sub>box</sub>/7</td><td>512<sup>3</sup><br>1024<sup>3</sup></td><td>2.70<br>2.70</td><td></td></tr><tr><td>L<sub>box</sub>/7</td><td>512<sup>3</sup><br>1024<sup>3</sup><br>2048<sup>3</sup></td><td>5.70<br>5.44<br>4.95</td><td>*</td></tr><tr><td>10</td><td>1</td><td>L<sub>box</sub>/7</td><td>512<sup>3</sup><br>1024<sup>3</sup><br>2048<sup>3</sup></td><td>5.70<br>5.44<br>4.95</td><td>*</td></tr></tbody></table>

In this subbox, while the volumic mean density $\bar{\rho}$ is not constant but exhibits small variations (of the order of 5–10%), the statistical mean density $\mathbb{E}(\rho)$ is constant. It is equal to the mass contained in the periodic box divided by its volume. Furthermore, the anisotropy introduced by the cubic geometry of the subbox does not affect the calculation of the correlation length. In fact, this anisotropy has an effect only at the boundary of the subbox and thus on scales much larger than the correlation length. The invariant that is tested is thus the one derived by Ref. [10] for which the statistical mean density $\mathbb{E}(\rho)(t)$ does not change with time and the flow is statistically isotropic, in contrast to the more general derivation of Ref. [11].

### C. Computation of the correlation length

To verify the time invariance of $M_{inv}$, a major difficulty is to compute the correlation length of the density field. Under the assumption of ergodicity and statistical homogeneity, the ACF of a random field $X$ can be calculated from its power spectrum $|\mathcal{F}[X - \mathbb{E}(X)]|^2$ using the Wiener-Khintchine theorem:

$$ \begin{aligned} C_X(q) &= \mathbb{E}[X(r)X(r+q)] - \mathbb{E}(X)^2 \\ &= \mathcal{F}^{-1}(|\mathcal{F}[X - \mathbb{E}(X)]|^2)(q), \end{aligned} \quad (8) $$

where $\mathcal{F}$ is the Fourier transform operator and $\mathbb{E}$ the statistical average. Numerically, however, the correlation function exhibits numerical noise at large distances which prevents its integral over the volume [see Eq. (2)] from converging. These spurious oscillations are probably due to the limited statistics at large distance.

To avoid this problem, various estimates of the correlation length have been suggested [11,20–22]. In Appendix B, we detail these possible estimates. However, none of them can be applied to the present case, either because they are too impractical or because they need too much computational resources. Here, to compute the correlation length, we choose to stop the volumic integration of the ACF at its first zero, assuming that the ACF at larger distance contains mostly noise. Due to the rapid decay of the ACF at large distance, the cutoff of ACF at its first zero does not impact too much the obtained value of the autocorrelation length.

To check that the value of the correlation length obtained by this method is accurate, we confront it to another method based on the value at 0 of the power spectrum [20]. However, we do not have access to this value because of the finite size of the box within which the measurements are done. To obtain this value, we thus have to impose a functional form of the power spectrum. We chose the following one with the slope of the inertial range $2\eta$ and the dissipation scale $\beta$:

$$ P_{\rho}^{1D}(k) \propto k^2 P_{\rho}^{3D}(k) \propto \frac{e^{-k\beta}}{\left(1 + \left(\frac{kL_{inj}}{2\pi}\right)^2\right)^{\eta}} k^2. \quad (9) $$

In Fig. 2, we compare the value of the correlation length computed from the fit of the power spectrum and the integration to the first zeros of the ACF. These computations are performed during the stationary regime of the $2048^3$ simulations. Both estimates are in good agreement with each other. A comparison between our measurement based on the integration to the first zeros of the ACF and on the ergodic
<figure><img src="image_4.png" /><figcaption>FIG. 2. Left: 1D density power spectrum (blue) and fit using Eq. (9) as a function of the wave number $\lambda = 2\pi/k$ for one of the output in 2048³ simulation at Mach 5 during the stationary regime. Right: correlation length measured in the 2048³ simulations by the integration ACF to its first zeros (red) and by the value at 0 of the power spectrum (black). The error bars correspond to the $\pm 1\sigma$ estimated from the time variations of the measured correlation length.</figcaption></figure>

theory is presented in Appendix B for simulations at 1024³. We have not made this comparison for our 2048³ because of the high computational cost needed to use the estimation of the correlation length with the ergodic method.

A key condition to be sure that the correlation length is well resolved is that it must be significantly smaller than the size of the system under study. From our numerical tests, we have verified that the correlation length is well resolved when it is at most 1/15 of the size of the system within which the measures are done (see Sec. II B). In addition, we find that the correlation length must be resolved by at least 10 cells in order to be properly converged (see Appendix A).

### III. TEMPORAL INVARIANCE OF $M_{inv}$

#### A. Steady state background

In this section we present the results of our numerical test of the time invariance of $M_{inv}$ in a steady state background for two different equation of states. We first consider an isothermal equation of state ($\gamma = 1$) before studying a polytropic one ($\gamma = 5/3$).

In Fig. 3, for the isothermal equation of state, we show that, during the turbulence decay phase, $M_{inv}$ remains relatively constant (less than 30% variation) compared to the decrease of the variance and the increase of the cube of correlation length, which vary between factors 2 and 10. We stop the integration of the simulation when the latter is no longer well resolved.

Fluctuations in the evaluations of the invariant reflect fluctuations in the estimates of statistical quantities. The biggest fluctuations come from the evaluations of the correlation length, which is notoriously delicate to estimate (see Sec. IV C and Ref. [23]). These fluctuations are the result of averaging over a finite volume. This is equivalent to saying that we have access to a finite number of independent samples from which to estimate the various statistics. The greater the ratio $l_c/L_{box}$, the greater the fluctuations of the estimates around the true statistical average.

Since the true statistical quantity $M_{inv}$ is known to be constant by definition during the stationary phase, the fluctuations of the estimate of $M_{inv}$ in Fig. 3 reflect the typical fluctuations generated by our estimator around this constant average and can be estimated. They are of the order of a factor 2 to 3. Since the fluctuations of the measured invariant during the decaying phase are smaller or of the order of those measured in the stationary phase, we can conclude that the true statistical quantity $M_{inv}$ is constant during the decaying phase.

In our numerical setup, the injection of turbulence at a scale smaller than the size of the box ensures that the volumic average over the measurement box of the quantities of interest is a good estimate of the true statistical value. This is confirmed by the small fluctuations of the mean density and variance during the statistically stationary regime before the turbulence driving is switched off. The relative fluctuations of these quantities are of the order of 5–10%. We compared this with a simulation in which the turbulence is driven at $L_{box}/2$ and $M = 7$. In that case, the relative fluctuations of the mean density and the density variance are of the order of 30% (see Fig. 4) and the variance $\sigma_{\rho}^2$ can vary by up to a factor of 2 during the simulation.

This demonstrates that turbulence injection at a scale significantly smaller than the box size is necessary to obtain good estimates of the statistical quantities of interest when performing volumic averaging. If the turbulence is injected at the box scale, an accurate estimate of the statistical mean can still be made by considering the time average over several turbulent crossing times. However, this method is limited to a statistically stationary turbulence regime. To study the statistical properties of nonstationary turbulent flows, e.g., decaying or self-gravitating and thus fragmenting, it is therefore of paramount importance to have sufficient statistical realizations in volume, i.e., $L_{inj} \ll L_{box}$, contrary to what is usually done in the literature where turbulence is usually injected at scales of the order of $L_{box}/2$. Their results may not be statistically accurate.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Evolution of mass $M_{inv}$ (blue), correlation length $l_c^3$ (orange), density variance $\text{Var}(\rho/\mathbb{E}(\rho))$ (red), and the Mach number (black dashed) for three Mach numbers in the 2048³ simulations. The equation of state is isothermal in these three simulations. The black dotted line corresponds to the time $T_{dec}$ at which the turbulence forcing is switched off. The integration is stopped when the correlation length is no longer well resolved and only the times after which the turbulence is fully developed are shown. This corresponds to $t > 2T_{cross}$ (see text). All the four quantities shown are normalized to their value at the time the turbulence forcing is switched off. The time is normalized to the time $T_{dec}$.</figcaption></figure>

Moreover, the invariant $M_{inv}$ is a dynamics invariant, independent of the equation of state, as it is simply based on the conservation of mass equation. For the sake of completeness, and for reasons that will become apparent in Sec. III B, we have also chosen to verify the invariance of $M_{inv}$ using a polytropic equation of state with $\gamma = 5/3$. It corresponds to the usual polytropic index for a monoatomic gas. In Fig. 3, and as in Fig. 3, we plot the evolution of the variance $\text{Var}(\rho/\mathbb{E}(\rho))$, the cube of the correlation length $l_c^3$, and the invariant $M_{inv}$. During the decay, each term exhibits statistical fluctuations but the invariant $M_{inv}$ remains constant with a standard deviation of 10%. This shows that the invariant is indeed conserved.

## B. Expanding or contracting medium

The derivation of Ref. [11] is more general than the one of [10] because it shows that $M_{inv}$ is still time invariant in an evolving background. The invariant should be conserved in a globally expanding or compressing medium in which $\mathbb{E}(\rho)$ changes with time.

To show that this is indeed the case, our consideration on the polytropic equation of state with $\gamma = 5/3$ is interesting

because it allows one to deduce the behavior of the invariant in an isotropically expanding or contracting box with velocity $V = H(t)r$. This type of flow is called Hubble flow in the astrophysical context and widely used in cosmology to describe the expansion of the universe. Reference [24] shows that the equations of conservation of mass, momentum, and energy have the same form in a nonexpanding medium with a polytropic equation of state with $\gamma = 5/3$ and in the comoving frame of an expanding medium if we choose the following transformation between the physical frame and the comoving frame:

$$ \tilde{r} = a^{-1}r, \quad d\tilde{t} = a^{-2}dt, \quad \tilde{\rho} = a^3\rho, \qquad (10) $$

where the tilde variables are those in the comoving frame and $a$ is the expansion factor. It can be related to the Hubble constant $H(t)$ by $H(t) = \frac{1}{a}\frac{da}{dt}$. In the comoving frame, $\mathbb{E}(\tilde{\rho})$ is constant in time and the problem in this frame is thus equivalent to the one already solved in Sec. III A with the polytropic equation of state.

Since the equations of conservation are the same for a polytropic gas with $\gamma = 5/3$ and in the comoving frame of an

<figure><img src="image_6.png" /><figcaption>FIG. 4. Evolution of the relative fluctuation of the mean density of the gas (left) and its variance (right) in the sub-box for a turbulence driving injected at $L_{inj} = L_{box}/7$ (solid) and $L_{inj} = L_{box}/2$ (dash) from the simulations at $M = 7$ and $512^3$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 5. Same as Fig. 3 but for a polytropic equation of state with $\gamma = 5/3$. The simulation is run at $M = 7$ with a 1024³ resolution.</figcaption></figure>

expanding/contracting flow, we also show from the polytropic simulation that the invariant is also verified in the comoving frame of an expanding/contracting medium. Thus we have

$$ \tilde{M}_{\text{inv}} = \mathbb{E}(\tilde{\rho})\text{Var}\left(\frac{\tilde{\rho}}{\mathbb{E}(\tilde{\rho})}\right) (\tilde{I}_c^\rho)^3 = \text{const.} \quad (11) $$

From the transformation given by Eq. (10) we have that $M_{\text{inv}} = \tilde{M}_{\text{inv}}$. Thus the invariant is also verified in the physical frame.

## IV. ANALYTICAL MODELING OF THE CORRELATION LENGTH AND THE INVARIANT

### A. Gaussian model for the log-density field: Predicting $I_c^\rho$

We now present a model for predicting the evolution of the correlation length with the Mach number. We start with the first order assumption that the log-density $s = \ln[\rho/\mathbb{E}(\rho)]$ is a 3D Gaussian random field. Whereas, strickly speaking, there is no exact physical or mathematical justification for the assumption of Gaussianity for $s$, a reasonable justification has been given in Ref. [25] for a gravoturbulent medium. The $s$ field is completely characterized by the following power spectrum:

$$ P_s^{\text{1D}}(k) \approx A \frac{8\pi \sigma_s^2 L_{\text{inj}}^3 e^{-k\beta/(2\pi)}}{\left(1 + \left(\frac{kL_{\text{inj}}}{2\pi}\right)^2\right)^2} k^2, \quad (12) $$

where $k$ is the wave vector, $\sigma_s^2$ is the variance of the log-density $s$, and $A$ is a coefficient close to unity

<figure><img src="image_8.png" /><figcaption>FIG. 6. In red, the 1D power spectrum of the log-density field measured during the stationary regime of the $M = 3.5 \times 1024^3$ simulation and the suggested model [see Eq. (12)] with the dissipation factor $\beta = 4\Delta x$ and without dissipation in blue ($\beta = 0\Delta x$). They are both similar to each other, especially in the inertial range. The gray vertical dotted line shows the injection length.</figcaption></figure>

that ensures that the integral of the power spectrum is equal to the variance of the log-density field. The coefficient $\beta$ corresponds to the dissipation scale. The dissipation range is found to be universal in the grid simulations where the dissipation is induced by the grid [26]. As a consequence, we chose a unique value of $\beta$ normalized to the spatial resolution $\Delta x$: $\beta = 4\Delta x$. The power spectrum given in Eq. (12) is plotted in red in Fig. 6.

This form reproduces particularly well the one measured in our numerical simulations as shown in Fig. 6. The factor $2\pi$ in Eq. (12) can be understood as follows. In the code, the forcing is performed at the wave vector $k_{\text{inj}} = 2\pi/L_{\text{inj}}$ with $L_{\text{inj}} = L_{\text{box}}/7$ in most of our simulations [see Eq. (6)]. It is quite natural that this particular scale is imprinted in some of the fields of the system. However, it is unclear why it affects in particular the log-density field. This point should be investigated in more detail in dedicated studies.

With the two above assumptions, i.e., that the log-density is a Gaussian random field and that its power spectrum has a form given by Eq. (12), the statistics of the $s$ field is specified. While the ACF of the log-density defined by Eq. (8) has no simple analytical form, we can still derive the correlation length of the log-density field:

$$ I_c^s = \frac{L_{\text{inj}}}{2\pi} \left( \frac{\pi^2}{2 \text{Ci}(\tilde{\beta})[\sin(\tilde{\beta}) + \tilde{\beta} \cos(\tilde{\beta})] + [\pi - 2 \text{Si}(\tilde{\beta})][\cos(\tilde{\beta}) - \tilde{\beta} \sin(\tilde{\beta})]} \right)^{1/3}, \quad (13) $$

where Ci and Si are respectively the cosinus and sinus integral function [27] and $\tilde{\beta} = \beta/L_{\text{inj}}$.

From the Gaussian assumption on the $s$ field, we get the correlation of the normalized density field $\tilde{\rho} = \rho/\mathbb{E}(\rho)$:

$$ C_{\tilde{\rho}}(\boldsymbol{q}) = \exp[C_s(\boldsymbol{q})] - 1. \quad (14) $$

This expression enables us to compute the correlation length of the density field:

$$ I_c^\rho = \left( \frac{\pi}{2} \frac{1}{C_{\tilde{\rho}}(0)} \int_0^\infty C_{\tilde{\rho}}(q) q^2 dq \right)^{1/3}. \quad (15) $$

This integral has no analytical expression. The variation of the density correlation length with $\sigma_s$ predicted by the model
<figure><img src="image_9.png" /><figcaption>FIG. 7. Evolution of the correlation length of the log-density field (left) and the density field (right) with the variance $\sigma_s^2$ of the log-density field. The model for the log-density correlation length is plotted as a dashed line [see Eq. (13)] in the left panel. The model for the density correlation length is plotted as a solid line [see Eq. (15)] in the right panel. The measurements of our 2048³ fiducial simulations are shown in red. The error bars correspond to the $\pm 1\sigma$ estimated from the time variations of the measured quantities.</figcaption></figure>

is plotted as a solid line in the right panel of Fig. 7. This model only holds when the turbulence is forced and fully developed. It does not give a prediction for the evolution of the correlation length of the density field when the turbulence is decaying, as the power spectrum of the log-density field would not have the functional form given by Eq. (12) anymore. Therefore, the evolution of the correlation length of the log-density field is not known.

We now compare the prediction of our model with the measurements of the correlation lengths made in the simulations after reaching the statistically stationary regime and before turbulence starts to decay. In Fig. 7 we have plotted the average correlation length of the density and log-density fields, respectively, measured in the four 2048³ runs with turbulence forced at scale $L_{box}/7$ and Mach numbers $M = 3.5, 5, 7$, and 10. We did not show the correlation length measured at $M = 2.5$ because it is not converged relative to the size of the box (see Appendix A). The error bars correspond to the uncertainties at $1\sigma$ estimated from the time variations of the measured quantities. As can be seen in the figure, our model agrees well (less than one $\sigma$) with the measured correlation lengths of the density and the log-density fields without any fitting parameters. Some discrepancies (around 1.5 to $2\sigma$) between the measurements and the model are observed at high Mach, especially for the log-density correlation length. We found that this is due to the flattening of the log-density power spectrum with Mach number, which is not taken into account in the model, but also to the fact that non-Gaussian features of the log-density field become non-negligible at high Mach. The convergence studies of the density and log-density correlation lengths are presented in Appendix A.

## B. Modeling the evolution of $M_{inv}$ with the Mach number

From the model of the correlation length presented in the previous section, we can model the dependence of $M_{inv}$ on the Mach number during the statistically stationary regime. Within the Gaussian approach, its expression can be derived

as

$$M_{inv} = \mathbb{E}(\rho)L_{inj}^3 \frac{\pi}{2C_{\bar{\rho}}(0)} \int_0^{\infty} C_{\bar{\rho}}(\tilde{q}L_{inj})\tilde{q}^2 d\tilde{q}(e^{\sigma_s^2} - 1). \quad (16)$$

This model is plotted in purple dotted line in Fig. 8. We see that this model does not reproduce at all the measurements of $M_{inv}$ made in the simulations after reaching the statistically stationary regime. This reflects the fact that the variance $\sigma_{\bar{\rho}}^2 = b^2 M^2$ of the density field predicted from the variance of $s$ under a Gaussian field assumption is significantly

<figure><img src="image_10.png" /><figcaption>FIG. 8. Evolution of the invariant $M_{inv}$ with the variance $\sigma_s^2$ of the log-density field. The model for the mass invariant based on the Gaussian assumption only is plotted as a purple dotted line. The model with a non-Gaussian correction for the density variance is plotted as a purple dashed line. Both rely on the model of the correlation length described in Sec. IV A. The model describing a constant mass invariant is plotted in black solid line (see Sec. V). The measurements from a 2048³ simulation with turbulence injected at $L_{box}/7$ are shown in red. The error bars correspond to the $\pm 1\sigma$ estimated from the time variations of the measured quantities.</figcaption></figure>
overestimated. Consequently, Eq. (16) overestimates $M_{inv}$, especially at high Mach.

The model proposed above has the major drawback of relying on the Gaussian assumption for the statistics of the log-density random field. It is well known that, although the PDF of the log-density field can be reasonably well modeled by a Gaussian function, the non-Gaussianities of turbulence reduce its high-density tail and thus its variance [12,18]. To better predict the evolution of the mass invariant with Mach number, we use the log-Poisson model to describe the log-density PDF [7,12,28,29] (see also Ref. [30] for another model). The log-density PDF has thus the following form:

$$ \mathcal{P}(s) = I_1(2\sqrt{\lambda}u)e^{-(\lambda+u)\sqrt{\frac{\lambda}{T^2u}}}, \quad (17) $$

defining the quantities

$$ u = \frac{\lambda}{1+T} - \frac{s}{T}, \quad \lambda = \frac{\sigma_s^2}{2T^2}. \quad (18) $$

This introduces a Mach-dependent parameter $T$ to reproduce the negative skewness of the density PDF observed in compressible turbulence simulations: the larger $T$, the larger the negative skewness. In this model, the variance of the density field becomes

$$ \sigma_{\bar{\rho}}^2 = \exp\left(\frac{\sigma_s^2}{1 + 3T + 2T^2}\right) - 1. \quad (19) $$

For $T=0$, we recover the Gaussian prediction for $\sigma_{\bar{\rho}}$. From a phenomenological model of the structure of supersonic shocks, a parametrization of the dependence of $T$ on the Mach number as $T = \kappa[1 - (bM)^{-2}]$ has been proposed [7]. We find that a value $\kappa = 0.16$ best fits our data, slightly different from $\kappa = 0.2$ found by the aforementioned authors. The Mach number can be related to the log-density variance by the following equation:

$$ \sigma_s^2 = \ln(1 + b^2 M^2), \quad (20) $$

with a forcing parameter $b$ whose value depends on the nature of the turbulent forcing [18,31,32]. These studies suggest that for mixed driving (solenoidal fraction $\zeta = 0.5$), $b$ should be close to 0.4. This value differs, however, from other studies [33]. Here we find that $b = 0.8$ reproduces best the relation we measure between $\sigma_s^2$ and $M$.

The expression of the invariant then becomes

$$ M_{inv} = \mathbb{E}(\rho)L_{inj}^3 \frac{\pi}{2C_{\bar{\rho}}(0)} \times \int_0^{\infty} C_{\bar{\rho}}(\tilde{q}L_{inj})\tilde{q}^2 d\tilde{q} \left[ \exp\left(\frac{\sigma_s^2}{1 + 3T + 2T^2}\right) - 1 \right]. \quad (21) $$

This non-Gaussian correction concerns only the relation between the density and the log-density variance, not the one between the density and the log-density correlation lengths, which seems to be well described by the Gaussian model proposed in Sec. IV A. Taking into account the non-Gaussianities in the relationship between the ACFs of the log-density and density field from a theoretical point of view is much more difficult because it would require the knowledge of all the orders of the two points correlation of the log-density field.

In Fig. 8, we plot the mass invariant against $\sigma_s^2$. Our model, which includes the above correction for the non-Gaussianity of the log-density PDF, reproduces the measurements reasonably well at high Mach, while it significantly underestimates the value at low Mach.

### C. Discussion of the hypotheses

As described above, our basic Gaussian model relies on a first order assumption about the statistics of the log-density field, namely that it is a Gaussian random field. This assumption is inaccurate because the supersonic shocks generate high order correlations, which induce phase coupling within the field (see Ref. [25] for a discussion) stronger at high Mach. Therefore, the statistics of $s$ field cannot be described only by its power spectrum and PDF. However, the fact that our Gaussian model successfully reproduces the evolution of the density correlation length, $l_c^\rho$, with the variance of the log-density field, $\sigma_s^2$, suggests that the underlying assumptions describe the statistics of the log-density field sufficiently well.

In Fig. 9 we compare the predicted ACFs of the density and log-density fields, respectively, with the ones measured in the three $2048^3$ runs at $M = 2.5, 3.5$, and $7$. As seen, there is good agreement, without any fitting parameters. As the Mach number increases, Eq. (14) becomes increasingly inaccurate to describe the two-point statistics of the density field. It must be kept in mind that the non-Gaussian correction to the model (see Sec. IV B) modifies only the variance of the density field but not the relation between the density and log-density ACFs [see Eq. (14)]. Thus, as shown in Fig. 9, the functional form of the density field ACF predicted by Eq. (14) deviates from the one obtained in the simulations at high $M$ because of non-Gaussian processes. However, this deviation has small consequences on the accuracy of the model to predict the correlation length (see Fig. 7) because the value of the integral of the ACF normalized to the density variance is not very sensitive to these small fluctuations of its functional form, at least for $M \lesssim 10$. More precisely, we can note that the deviation of the model compared to the measured one is large at short distance, but much smaller at large one. This is expected since the non-Gaussianities will play an important role at short distance when the structures are highly correlated. They will be less important at large distance, where the fields are almost independent. Because the correlation length involves the integral of $q^2C_\rho(q)$, the discrepancies observed between the model and the data play a minor role in the result of the integral. This explains the success of the model in describing the evolution of the correlation lengths as the function of $\sigma_s$ and $M$.

### V. DERIVATION OF THE LAW $T(\sigma_s)$

As an example of what we can learn from this invariant, we show that we can compute the variation of the parameter $T$ with the Mach number and thus predict some of the non-Gaussian properties of the density field.

#### A. Dependence of the mass invariant $M_{inv}$ with the Mach number

In this subsection we first give arguments showing that the quantity $M_{inv}$ should be independent of the Mach number.
<figure><img src="image_11.png" /><figcaption>FIG. 9. Autocovariance function of log-density s (top) and density ρ (bottom) normalized to their respective variance measured in the $M = 2.5, 3.5$, and 7 runs during the statistically stationary phase. Our model for the log-density and density autocorrelation functions is shown by the red dashed line. The measured value of the correlation length of the log-density and density fields is shown by the vertical black dashed line. The y scale is log above $e^{-1}$ and linear below.</figcaption></figure>

Based on this, in the next subsection, we will derive the evolution of the log-density PDF with $\sigma_s$ and $M$ and compare it with the existing parametrizations of this relation in the literature.

Let us consider a medium that is initially at rest, in which turbulence is driven. Once the correlation length is small enough compared to the box size to be properly resolved, the mass $M_{inv}$ is a constant quantity. In particular, when we stop driving the turbulence and let it decay, this quantity remains constant. It only breaks down when the correlation length is too large compared to the size of the system, i.e., when volu- mic averages are bad estimates of statistical averages. Instead of letting the turbulence decay, we can go the other way and increase the driving force until a new target Mach number is reached. During this process, none of the hypotheses on which the invariant is based are violated. The quantity $M_{inv}$ should therefore remain constant during this process. Since this can be done for any Mach number, we conclude that the invariant should be constant for any Mach number.

However, by considering the limit $M \to 0$ we can see that this reasoning is limited. In fact, at very low Mach, the variance of the density field tends to 0, while the cor- relation length cannot go to infinity because it is limited by the injection length. Indeed, regions of the fluid sepa- rated by distances larger than injection length are completely uncorrelated. This can be seen from the shape of the 3D power spectrum of the density which is flat for $k < k_{inj}$ with $k_{inj} = 2\pi/L_{inj}$. Thus, at least in forced subsonic turbulence, the invariant should depend on the Mach number. This de- pendence on the Mach number in the subsonic regime is the consequence of the fact that the cross-correlation between the density and the momentum is not negligible due to the large

scale correlations imposed on the flow by the driving. Once we reach the supersonic regime, the cross-correlation becomes negligible.

This behavior can be seen by considering the accretion timescale $\tau_{acc}$ of the invariant, corresponding to the time needed to modify it substantially. From the computation of the invariant derived by [10,11], it is defined as

$$ \tau_{acc} = \left| \frac{\int_{\Omega} C_{\rho}(q) d^3 q}{\int_{\partial \Omega} \mathbb{E}(\rho \rho' \mathbf{v}') d^2 \mathbf{S}} \right|, \quad (22) $$

where $d^2\mathbf{S}$ is an infinitesimal surface element and where $\Omega$ is the measurement domain, typically the sub-box of size $L_{box}/2$. To ensure that the invariant is constant, the accretion time $\tau_{acc}$ should be greater than the turbulent crossing time $T_{cross} = L_{inj}/\sqrt{\mathbb{E}(\mathbf{v}^2)}$. In Fig. 10, we plot the measurement of $\tau_{acc}$ done in the $1024^3$ run divided by the crossing time during the transitory regime. We do not perform this calculation on a $2048^3$ simulation because of the high computational cost of this calculation, but to study a $1024^3$ simulation is sufficient given the convergence study (see Appendix A). $\Omega$ is taken to be a ball whose radius is the lag corresponding to the first zeros of the ACF. We see that these measurements cross the horizontal line $\tau_{acc}/T_{cross} = 1$ at $\sigma_s^2 \simeq 1$. $M_{inv}$ should thus be constant soon after the time such that $\sigma_s^2 = 1$, i.e., in the supersonic regime, as observed in the simulation.

To compute the value of the invariant, we focus on the transonic regime where the statistical properties of the density field are reasonably well known. For a Mach number close to 1, the density PDF is well approximated by a log normal (e.g., Ref. [12]), i.e., $T \ll 1$. Basically, we want to choose the largest supersonic Mach number (or largest $\sigma_s^2$) at which $T$ is
<figure><img src="image_12.png" /><figcaption>FIG. 10. Black dots show the measurement of $\tau_{acc}/T_{cross}$ in the 1024³ run during the transitory regime for different target Mach numbers: $M = 3.5$ (black), $M = 5$ (green), and $M = 7$ (blue). It is evaluated at the first zeros of the ACF $C_{\rho}$. The ratio is small in the subsonic regime ($\sigma_s^2 \lesssim 1$), suggesting that the invariant in this regime is strongly dependent on the Mach number. In the supersonic regime ($\sigma_s^2 \gtrsim 1$) it increases significantly, confirming that the invariant will be independent of the Mach number in this regime.</figcaption></figure>

still small, typically below a few percent, such that the s-PDF is still well approximated by a Gaussian and not to be affected by the variation of the invariant in the subsonic regime. In Eq. (16), we choose $\sigma_s^2 = 1.3$ because it corresponds to a Mach number of 2. This choice is motivated by the fact that, at such Mach number, the measurements of $T$ available in the literature report a value between 1% and 5%. Moreover, $\sigma_s^2 > 1$ suggests from Fig. 10 that $M_{inv}$ should already be constant. With this value, we find that the mass of the invariant is equal to $\alpha M_{inv}$ with $\alpha = 2.0 \times 10^{-2}$ and $M_{inv} = \mathbb{E}(\rho)L_{inj}^3$. This corresponds to the solid black horizontal line in Fig. 8. The fact that the invariant is constant for different Mach numbers is consistent with the numerical data plotted in Fig. 8.

Our analytical model slightly underestimates its value because the model of the correlation length also slightly underestimates the measured one and $M_{inv}$ depends steeply on $l_c^\rho$. This is a consequence of the non-Gaussianities not taken into account by the model.

## B. Computation of the parameter $T$

In the framework described above, we can determine the evolution of the parameter $T$ introduced by [12] for compressible turbulence versus $\sigma_s$ by solving the following equation derived from Eq. (21) using $M_{inv} = \alpha \rho L_{inj}^3$

$$ \alpha = \frac{\pi}{2C_{\rho}(0)} \int_{0}^{\infty} (\exp(C_s(\bar{q}L_{inj})) - 1) \bar{q}^2 d\bar{q} \\ \times \left[ \exp\left(\frac{\sigma_s^2}{1 + 3T + 2T^2}\right) - 1 \right], \quad (23) $$

where $\alpha$ is the constant determined in the previous section. The determination of the parameter $T$ by this method assumes that the previously derived model for the density correlation length is accurate even at high Mach numbers. This is confirmed by our numerical simulations, at least for Mach numbers below 8 (see Sec. IV).

To predict the variation of $T$ with the Mach number instead of the log-density variance $\sigma_s^2$, we use Eq. (20). In the left panel of Fig. 11 we plot the evolution of $T$ with the Mach number. They are obtained by fitting the s-PDF measured in our simulation with the log-Poisson PDF [see Eq. (17)] with $T$ being the only fit parameter. As seen in most of the numerical simulations, we see that this parameter increases with the Mach number. We compare this with other laws for $T$ derived from numerical simulations and/or phenomenological models [7,12,33]. All these models are parametrized with the compressive Mach number $M_c = bM$.

The model derived with the value of the invariant computed at Mach 2 is plotted in blue solid line. Because this value is subject to uncertainties, the shaded area corresponds to the predicted of the model with an invariant calculated at a Mach between 1.6 and 2.5, corresponding to $\sigma_s^2$ varying between 1 and 1.6. The best fit to our simulations, including that of F21 and excluding $M = 10$ (see below), is obtained for $M = 2$. The model reproduces well the measurements of $T$ in the high resolution simulation performed by Ref. [34] and the measurements in our 2048³ simulations for Mach numbers below 7, while the existing parametrizations in the literature overestimate it [7,12,33]. This may suggest that the measurements of $T$ in some of the low resolution numerical simulations are not well converged even at small Mach number because of a lack of statistics, the turbulence being usually injected at $L_{box}/2$. The parametrizations given by [12,33] are also approximate, given the large scatter in the measurements of $T$ used in these studies to determine these parametrizations. For Mach number larger than 7, the underestimation of $T$ predicted by the model compared to our measurement in the Mach 10 simulation and to the existing parametrizations is likely to be a consequence of the non-Gaussianities that become dominant at such high Mach numbers. These non-Gaussian features result in an underestimation of the correlation length by the model described in Sec. IV A. The range of validity of our model is therefore limited to Mach numbers between 2 and 7. Below Mach 2, $M_{inv}$ is not constant while, above Mach 7, non-Gaussian features are not negligible.

On the right panel, we show the s-PDF measured in our 2048³ simulation at $M = 5$ and plot the PDF predicted by our model. The model reproduces very well the data. For comparison, we also plot the PDF that would be given by the Gaussian model, i.e., with $T = 0$. It significantly overestimates the amount of dense gas, highlighting the importance of considering a model able to describe the negative skewness of the s-PDF. To predict accurately the amount of dense gas in a supersonic turbulent medium is indeed of prime importance in an astrophysical context, namely to model the star formation process [33,35].

## VI. CONCLUSION

In this article, we have studied numerically the accuracy of the mass invariant introduced by Ref. [10] for compressible isotropic turbulence and extended by Ref. [11]. We show that the predicted time invariance of the quantity $M_{inv}$ is well verified during the decaying phase of isotropic turbulence. A correlation length small enough compared to the box size is crucial for a correct computation of the invariant. In
<figure><img src="image_13.png" /><figcaption>FIG. 11. Left: evolution of the parameter *T* with the compressible Mach number *M*<sub>c</sub>. Different parametrizations of *T* found in the literature are shown in dotted green [12], orange [7], and red [33]. The red points correspond to the *T* parameter extracted by fitting the measured PDF with the functional form suggested by [12] [his Eq. (5)] in our four 2048³ simulations with turbulence forced at *L*<sub>box</sub>/7. The gray point is the measurement from the high resolution turbulence simulation from [34]. The error bars correspond to the ±1σ estimated from the time variations of the measured quantities. The model is shown in blue. The shaded area corresponds to different predictions of the model considering various values of the invariant. Top right: in green, measured *s*-PDF in the 2048³ simulation at *M* = 5. The shaded area corresponds to the 1σ uncertainties estimated from the time variations of the PDF. The blue solid line is the PDF predicted by the model presented in Sec. V. The dotted black line corresponds to the log normal mode of the PDF (*T* = 0) and the dashed line corresponds to the PDF predicted by [7]. Bottom right: relative error between each model plotted on the top panel and the measured PDF.</figcaption></figure>

order to fulfill this condition, turbulence must be injected at
a scale significantly smaller than the size of the system under
study. Although the calculations of Ref. [11] for the invariant
are valid for both isotropic and nonisotropic turbulence, we
have studied here only the isotropic case which is represen-
tative of a vast variety of turbulent flows both in physics
and astrophysics. Anisotropic turbulence does not present any
particular difficulties, but is outside the scope of this paper and
will be considered in a future study.

Even though we use the Gaussian random field assumption
which is a first order assumption to describe the statistical
properties of the supersonic flow, the derived model repro-
duces well the evolution of the correlation length l_c^ρ with
M. This suggests that the Gaussian random field assumption
describes the statistics of the log-density field sufficiently
well. However, in order to correctly describe the evolution of
M inv, we have derived analytically a non-Gaussian correction
to better capture the evolution of the density variance σ^2
with M. This improved model better agrees with the results
of the simulations. Finally, we argue why M inv should not
depend on the Mach number from the trans to supersonic
regime (M ~ 7) and show that this invariant can be a pow-
erful tool to predict the evolution of the log-density PDF
with the Mach number with no free parameters. This is quite
a fundamental result for our understanding of compressible
turbulence as it provides an exact way to infer the evolution

of the main statistical properties of a compressible turbulent
flow.

In the astrophysical context, gravity plays a major role in
the formation of structures in gravoturbulent molecular clouds
[36,37]. It modifies the statistical properties of the flow, as
seen in the form of the density PDF, which develops a power-
law tail at high density [38,39]. In that context, the invariant
may represent a fundamental quantity for the formation of
structures in molecular clouds, as it would determine the
characteristic mass of the most correlated structures induced
by turbulence in a gravitational system. On the theoretical
point of view, this has been studied by [11]. It will be thus
of prime importance to verify the validity of this invariant and
its physical interpretation in such a context, which will be the
purpose of a forthcoming study.

ACKNOWLEDGMENTS

The authors thank the referees for helpful comments that
improved the clarity of the manuscript. The authors thank L.
Chevillard for reading the manuscript and his helpful com-
ments. They also thank Q. Vigneron, E. Lynch, T. Gillet,
N. Brucy, M. Bourgoin, B. Commerçon, and L. Cherry for
helpful discussions and comments. Simulations were pro-
duced using GENCI allocations (Grants No. A0150411111
and No. A0170411111). We gratefully acknowledge support
from the CBPSmn (PSMN, Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. The platform operates the SIDUS solution [40] developed by E. Quemener. P.D. and J.F. acknowledge support from Agence Nationale pour la Recherche through the BRIDGES project (ANR-23-CE31-0005).

## DATA AVAILABILITY

The data that support the findings of this article are not publicly available upon publication because it is not technically feasible and/or the cost of preparing, depositing, and hosting the data would be prohibitive within the terms of this research project. The data are available from the authors upon reasonable request.

## APPENDIX A: NUMERICAL CONVERGENCE

In this Appendix we present our convergence test for the different quantities of interest, namely the density and log-density variances and the correlation lengths. First, we have increased the resolution of the simulation from 512³ to 2048³ for the simulation where turbulence is forced at $L_{box}/7$. As shown in Fig. 12, the variances and correlation length are well converged for the 2048³ simulation. Doubling the number of cells from 1024³ to 2048³ changes the correlation length by less than 5%. This corresponds to the case where the correlation is resolved by more than 10 cells. At high Mach (Mach 7 and 10) the log-density variance shows a slightly larger variation, until 15% when the resolution is doubled. To properly confirm the convergence of the log-density variance at these Mach numbers, a 4096³ should be performed, but this requires too much computational resources.

To test the effect of the finite volume of the studied system, we also increased the size of the subbox relative to the injection length. Therefore, we have performed 1024³ and 2048³

simulations where the turbulence is forced at $L_{box}/14$. The physical spatial resolution of this simulation is the same as in the 512³ and 1024³ run with turbulence forced at $L_{box}/7$, respectively. When we force the turbulence at $L_{box}/14$ with a resolution of 1024³ at $M = 3.5$, the log-density variance and the density correlation length are similar to those measured in the 512³ simulation with turbulence forced at Mach $L_{box}/7$. This confirms that the size of the system has little effect on these quantities, provided it is large enough. However, at $M = 2.5$ the density correlation length is not converged at all with respect to the injection length. This suggests that at such a small Mach number one needs a very large box to properly resolve the correlation length of the fields of interest. To do this, one should perform a 4096³ simulation with an injection of turbulence at $L_{box}/28$ to confirm the convergence of the simulations run at $L_{box}/14$, but this also requires too much computational resource. This problem highlights the difficulty of this study, which requires both a very good spatial resolution to properly resolve small scale turbulent structures and a turbulent box large enough compared to the injection scale of turbulence to have sufficient statistics.

## APPENDIX B: USUAL ESTIMATES OF THE CORRELATION LENGTH

To estimate the correlation length, several other methods have been suggested in the literature. The first one consists in considering the integral length [21] as a proxy of the correlation length. The integral length is defined as

$$l_i^\rho = \frac{1}{C_\rho(0)} \int_0^\infty C_\rho(q) dq. \quad (B1)$$

As justified by Ref. [11], for appropriate functional forms of correlation functions, $l_i^\rho \simeq l_c^\rho$. However, this estimate is only relevant for isotropic turbulence because it is based on the azimuthally averaged ACF $C_\rho(q)$.

<figure><img src="image_14.png" /><figcaption>FIG. 12. Evolution of the correlation length of the log-density field (left) and the density field (right) with the variance $\sigma_s^2$ of the log-density field. The model for the log-density correlation length is plotted as a dashed line [see Eq. (13)] in the left panel. The model for the density correlation length is plotted as a solid line [see Eq. (5)] in the right panel. The measurements from the 512³ runs are shown in green, those from the 1024³ runs in blue, and those from the 2048³ runs in red. The measurements from the 1024³ and 2048³ with turbulence injected at $L_{box}/14$ are shown in purple and brown, respectively. The error bars correspond to the $\pm 1\sigma$ estimated from the time variations of the measured quantities.</figcaption></figure>
To avoid the computation of the ACF of the density field, which can be numerically expensive in high resolution simulations, we can also use the following estimate of the integral length [11,22]:

$$l_i^\rho = \frac{L_{\text{box}}}{4} \frac{\text{Var}(\frac{\Sigma}{\bar{\Sigma}})}{\text{Var}(\frac{\rho}{\bar{\rho}(\rho)})}, \quad (B2)$$

where $\Sigma$ is the column density field along a given line of sight and $\bar{\Sigma}$ is its spatially averaged value. In a subbox of size $L_{\text{box}}/2$ with Cartesian coordinates, it is defined by

$$\Sigma(x, y) = \int_{0}^{\frac{L_{\text{box}}}{2}} \rho(x, y, z) dz. \quad (B3)$$

Again, this estimate is only relevant for isotropic turbulence. For anisotropic turbulence, the result will depend on the choice of the line of sight. We have checked that the two estimates given by Eqs. (B1) and (B2) give the same result.

However, the ratio between the correlation length and the integral length can be as large as 2, which means that the quantity $M_{\text{inv}}$ can be over- or underestimated by almost an order of magnitude. Furthermore, the ratio $l_c^\rho/l_i^\rho$ can vary with the Mach number. To calculate $M_{\text{inv}}$ with the quantity $l_i^\rho$ would lead to an apparent violation of the time invariance of this quantity. The study of the properties of $M_{\text{inv}}$ calculated from the integral length is therefore useless.

The consistency of the estimate of the correlation length based on the integral of the ACF until its first zero can also be checked by comparing it with other estimates of the correlation length based on ergodic theory. Reference [20] shows

that the correlation length can be estimated by

$$l_c^\rho = \frac{1}{2} \frac{L_{\text{box}}}{2} \left( \frac{\text{Var}(\bar{\rho}(t))}{\text{Var}(\rho)} \right)^{1/3}, \quad (B4)$$

where $\text{Var}(\bar{\rho}(t))$ is the variance of the temporal variation of the volumic estimate of the mean density over time and $\text{Var}(\rho)$ is the statistical variance of the density field.

This ergodic estimate of the correlation length can only be used for stationary systems. It cannot therefore be used to check the invariance of $M_{\text{inv}}$ in decaying turbulent flows. Furthermore, the convergence of the variance of the statistical variation of the mean density within the box is difficult to achieve. It is necessary to integrate the simulation over a very long time (more than several hundred turbulent crossing times) to have an accurate estimate of this quantity, which is not feasible for the high resolution simulations used in this study.

We also made a comparison in our 1024³ simulations between the correlation length of the density field computed from the integration of the ACF to its first zeros and that computed from the ergodic theory [Eq. (B4)] for Mach 3.5 and Mach 10. To obtain a converged result we had to run the simulations (with turbulent forcing) during 70 crossing times. At Mach 3.5 we get $l_c^{\text{1st 0}}/L_{\text{inj}} = 0.168 \pm 0.027$ and $l_c^{\text{ergo}}/L_{\text{inj}} = 0.166$. At Mach 10 we get $l_c^{\text{1st 0}}/L_{\text{inj}} = 0.126 \pm 0.020$ and $l_c^{\text{ergo}}/L_{\text{inj}} = 0.130$. Note that we have no easy way of estimating the uncertainties on $l_c^{\text{ergo}}$. These estimates are in very good agreement with each other, confirming the robustness of our calculation of the correlation length using the first zeros of the ACF.

[1] B. J. Bayly, C. D. Levermore, and T. Passot, Density variations in weakly compressible flows, Phys. Fluids A 4, 945 (1992).

[2] J.-P. Bertoglio, F. Bataille, and J.-D. Marion, Two-point closures for weakly compressible turbulence, Phys. Fluids 13, 290 (2001).

[3] Y. Zhou, Turbulence theories and statistical closure approaches, Phys. Rep. 935, 1 (2021).

[4] J. Wang, Y. Yang, Y. Shi, Z. Xiao, X. T. He, and S. Chen, Statistics and structures of pressure and density in compressible isotropic turbulence, J. Turbul. 14, 21 (2013).

[5] J. Wang, M. Wan, S. Chen, and S. Chen, Kinetic energy transfer in compressible isotropic turbulence, J. Fluid Mech. 841, 581 (2018).

[6] E. Vazquez-Semadeni, Hierarchical structure in nearly pressureless flows as a consequence of self-similar statistics, Astrophys. J. 423, 681 (1994).

[7] J. Squire and P. F. Hopkins, The distribution of density in supersonic turbulence, Mon. Not. R. Astron. Soc. 471, 3753 (2017).

[8] S. Galtier and S. Banerjee, Exact relation for correlation functions in compressible isothermal turbulence, Phys. Rev. Lett. 107, 134501 (2011).

[9] We recall that a stochastic field is statistically homogeneous if all its statistical properties are invariant under space translations. Not to be confused with spatial homogeneity.

[10] S. Chandrasekhar, The fluctuations of density in isotropic turbulence, Proc. R. Soc. London A 210, 18 (1951).

[11] E. Jaupart and G. Chabrier, Generalized transport equation for the autocovariance function of the density field and mass invariant in star-forming clouds, Astrophys. J. Lett. 922, L36 (2021).

[12] P. F. Hopkins, A model for (non-lognormal) density distributions in isothermal turbulence, Mon. Not. R. Astron. Soc. 430, 1880 (2013).

[13] R. Teyssier, Cosmological hydrodynamics with adaptive mesh refinement: A new high resolution code called RAMSES, Astron. Astrophys. 385, 337 (2002).

[14] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Springer, Berlin, Heidelberg, 2009).

[15] V. Eswaran and S. Pope, An examination of forcing in direct numerical simulations of turbulence, Comput. Fluids 16, 257 (1988).

[16] W. Schmidt, W. Hillebrandt, and J. C. Niemeyer, Numerical dissipation and the bottleneck effect in simulations of compressible isotropic turbulence, Comput. Fluids 35, 353 (2006).

[17] W. Schmidt, C. Federrath, M. Hupp, S. Kern, and J. C. Niemeyer, Numerical simulations of compressively driven interstellar turbulence: I. Isothermal gas, Astron. Astrophys. 494, 127 (2009).

[18] C. Federrath, J. Roman-Duval, R. S. Klessen, W. Schmidt, and M.-M. Mac Low, Comparing the statistics of interstellar
turbulence in simulations and observations: Solenoidal versus
compressive turbulence forcing, *Astron. Astrophys.* **512**, A81
(2010).

[19] C. Federrath, On the universality of supersonic turbulence, Mon. Not. R. Astron. Soc. 436, 1245 (2013).

[20] E. Jaupart and G. Chabrier, Statistical properties and correlation length in star-forming molecular clouds: I. Formalism and application to observations, *Astron. Astrophys.* **663**, A113 (2022).

[21] G. K. Batchelor, *The Theory of Homogeneous Turbulence* (Cambridge University Press, Cambridge, UK, 1953).

[22] S. Bialy and B. Burkhart, The driving scale–density decorrelation scale relation in a turbulent medium, Astrophys. J. Lett. 894, L2 (2020).

[23] A. Papoulis and S. U. Pillai, *Probability, Random Variables, and Stochastic Processes*, 4th ed. (McGraw-Hill, New York, 2002).

[24] H. Martel and P. R. Shapiro, A convenient set of comoving cosmological variables and their application, Mon. Not. R. Astron. Soc. 297, 467 (1998).

[25] P. Dumond and G. Chabrier, Impact of the shape of the prestellar density fluctuations on the core mass function, Mon. Not. R. Astron. Soc. 536, 2449 (2024).

[26] C. Federrath and R. S. Klessen, On the star formation efficiency of turbulent magnetized clouds, Astrophys. J. 763, 51 (2013).

[27] The cosine and sine integral function are defined as follows: $Ci(z) = - \int_z^\infty \frac{\cos(t)}{t} dt$ and $Si(z) = \int_0^z \frac{\cos(t)}{t} dt$.

[28] B. Dubrulle, Intermittency in fully developed turbulence: Log-Poisson statistics and generalized scale covariance, Phys. Rev. Lett. 73, 959 (1994).

[29] B. Castaing, The temperature of turbulent flows, J. Phys. II 6, 105 (1996).

[30] P. Mocz and B. Burkhart, A Markov model for non-lognormal density distributions in compressive isothermal turbulence, Astrophys. J. Lett. 884, L35 (2019).

[31] F. Z. Molina, S. C. O. Glover, C. Federrath, and R. S. Klessen, The density variance-Mach number relation in supersonic

ence - I. Isothermal, magnetized gas, Mon. Not. R. Astron. Soc. 423, 2680 (2012).

[32] L. Konstandin, W. Schmidt, P. Girichidis, T. Peters, R. Shetty, and R. S. Klessen, Mach number study of supersonic turbulence: The properties of the density field, Mon. Not. R. Astron. Soc. 460, 4483 (2016).

[33] N. Brucy, P. Hennebelle, T. Colman, R. S. Klessen, and C. L. Yhuelic, Inefficient star formation in high Mach number environments. II. Numerical simulations and comparison with analytical models, Astron. Astrophys. 690, A44 (2024).

[34] C. Federrath, R. S. Klessen, L. Iapichino, and J. R. Beattie, The sonic scale of interstellar turbulence, Nat. Astron. 5, 365 (2021).

[35] P. Hennebelle, N. Brucy, and T. Colman, Inefficient star formation in high Mach number environments I. The turbulent support analytical model, Astron. Astrophys. 690, A43 (2024).

[36] E. Vázquez-Semadeni, A. Palau, J. Ballesteros-Paredes, G. C. Gómez, and M. Zamora-Avilés, Global hierarchical collapse in molecular clouds. Towards a comprehensive scenario, Mon. Not. R. Astron. Soc. 490, 3061 (2019).

[37] P. Hennebelle and M. Y. Grudić, The physical origin of the stellar initial mass function, Annu. Rev. Astron. Astrophys. 62, 63 (2024).

[38] E. Jaupart and G. Chabrier, Evolution of the density pdf in star-forming clouds: the role of gravity, Astrophys. J. 903, L2 (2020).

[39] N. Schneider, V. Ossenkopf-Okada, S. Clarke, R. S. Klessen, S. Kabanovic, T. Veltchev, S. Bontemps, S. Dib, T. Csengeri, C. Federrath, J. Di Francesco, F. Motte, Ph. André, D. Arzoumanian, J. R. Beattie, L. Bonne, P. Didelon, D. Elia, V. Könyves, A. Kritsuk et al., Understanding star formation in molecular clouds: IV. Column density PDFs from quiescent to massive molecular clouds, Astron. Astrophys. 666, A165 (2022).

[40] E. Quemener and M. Corvellec, Sidus—the solution for extreme deduplication of an operating system, Linux J. 2013, 3 (2013).