Generalized Transport Equation for the Autocovariance Function of the Density Field and Mass Invariant in Star-forming Clouds

The dynamics of star-forming molecular clouds (MCs) are determined by the statistical properties of their density fluctuations under turbulence and gravity. A fundamental quantity in such a study is the autocovariance function (ACF) of density-field fluctuations, which allows the determination of the characteristic correlation length l_c (ρ) of density structures within the cloud. In this Letter, we study the ACF and the correlation length of density fluctuations in MCs and we show that the latter can be identified as the average size of the most correlated structures within the cloud. Generalizing the transport equation derived by Chandrasekhar (1951a) for static, homogeneous isotropic turbulence to a non-isotropic, time-evolving turbulent flow, we show that, whereas the correlation length decreases with time as gravity proceeds in the cloud, the mass contained within these structures of size l_c (ρ) is an invariant, like invariants found, e.g., in incompressible turbulence (Batchelor 1953). This striking result implies that the average mass contained in the most correlated structures in star-forming gravo-turbulent MCs, which will be ultimately distributed within prestellar cores, is imprinted within the initial conditions of the cloud and is constant during its evolution.

2.1. Evolution of a Molecular Cloud

The dynamics of the cloud are described as in Jaupart & Chabrier (2020; hereafter JC20). The only useful equation for the present study is mass conservation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot \left(\rho {\boldsymbol{v}}\right)=0,\end{eqnarray} \tag{ 1 }$$

where ρ denotes the gas mass density and v the velocity field.

We are interested in clouds that will eventually condense locally to form stars, hence we separate the evolution of the background from that of local density deviations. The velocity field v is thus split into a mean velocity V and a (turbulent) velocity u (Ledoux & Walraven 1958). Introducing the logarithmic excess of density, $s=\mathrm{log}(\rho /\overline{\rho })$, we get by definition:

$$\begin{eqnarray}&&{\boldsymbol{V}}\equiv \displaystyle \frac{1}{\overline{\rho }}\overline{\rho {\boldsymbol{v}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\boldsymbol{u}}\equiv {\boldsymbol{v}}-{\boldsymbol{V}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\rho \equiv \overline{\rho }({\boldsymbol{x}},t)\,{e}^{s},\end{eqnarray} \tag{ 4 }$$

where $\overline{{\rm{\Phi }}}({\boldsymbol{x}},t)\equiv {\mathbb{E}}\left({\rm{\Phi }}\right)({\boldsymbol{x}},t)$ is the mathematical expectation, also called statistical average or mean, of random field Φ (Pope 1985; Frisch 1995). We note that $\overline{{\boldsymbol{u}}}\ne 0$ a priori but $\overline{\rho {\boldsymbol{u}}}=0$. This ensures that on average there is no transfer of mass due to turbulence and the equation of continuity (1) remains valid for the mean field,

$$\begin{eqnarray}&&\displaystyle \frac{\partial \overline{\rho }}{\partial t}+{\boldsymbol{\nabla }}\cdot \left(\overline{\rho }\ {\boldsymbol{V}}\right)=0.\end{eqnarray} \tag{ 5 }$$

Subtracting the equations for the average variables from the original equations, we obtain the evolution of the density deviations.

2.2. Model for the Statistics of a Turbulent Cloud

2.2.1. Statistically Homogeneous Clouds

In studies of star formation, be it observations of a cloud or numerical simulations, one usually has access to only a small number of samples (only one in most cases). Thus, one has to make a basic assumption, sometimes called "fair-sample hypothesis," that the observed sample is large enough for volumetric (or time) averages over this single sample to provide accurate statistical estimates. For this procedure to be valid, the random field must be ergodic and thus statistically homogeneous (Papoulis & Pillai 1965). Note that statistical homogeneity does not imply spatial homogeneity. Ergodicity, one of the fundamental hypotheses of statistical physics, insures that the average value of a statistical quantity (density fluctuations in the present context) is equal to the mean of a large number (in space or time) of measured quantities (e.g., Penrose 1979). It is commonly made, for instance, in studies of turbulent flows, with or without self-gravity (Chandrasekhar 1951a, 1951b; Batchelor 1953; Pope 1985; Frisch 1995; Pan et al. 2018, 2019a, 2019b; Jaupart & Chabrier 2020) or in cosmology to study the dynamical evolution of structures in the universe (Peebles 1973; Heinesen 2020). This assumption does not constrain fluctuations around the average to be small. Statistical homogeneity implies that, for any stochastic field Φ, $\overline{{\rm{\Phi }}({\boldsymbol{x}},t)}=\overline{{\rm{\Phi }}}(t)$. In particular, $\overline{\rho }({\boldsymbol{x}},t)=\overline{\rho }(t)$ in our context.

With these assumptions, the dynamics of the cloud density and logarithm of density fluctuations are governed by the following equations:

$$\begin{eqnarray}&&-\displaystyle \frac{d\mathrm{ln}(\overline{\rho })}{{dt}}=-\displaystyle \frac{1}{\overline{\rho }}\displaystyle \frac{d\overline{\rho }}{{dt}}={\boldsymbol{\nabla }}\cdot {\boldsymbol{V}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{Ds}}{{Dt}}=-{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 7 }$$

where $\tfrac{d}{{dt}}$ denotes the derivative of a variable that is only a function of time t and $\tfrac{D}{{Dt}}=\tfrac{\partial }{\partial t}+\left({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}\right)$ is the Lagrangian derivative. Equation (6) shows that the statistical homogeneity hypothesis for ρ implies, to be consistent, that the rhs of Equation (6) must be a function of time t only. It thus constrains the flow to belong to a certain class of flows. In order to fulfill this constraint, it suffices that:

$$\begin{eqnarray}&&\boxed{{\boldsymbol{V}}({\boldsymbol{x}},t)=\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{x}}+{{\boldsymbol{c}}}_{V}(t),}\end{eqnarray} \tag{ 8 }$$

where $\underline{\underline{{L}_{V}}}(t)$ is a 3 × 3 matrix and c _V (t) is a spatially constant vector. Enforcing V = 0 yields exactly the equations usually used to prescribe the evolution of a periodic simulation box in an astrophysical context (Federrath & Klessen 2012; Pan et al. 2019b). However, this is not equivalent to applying periodic boundary conditions (see, e.g., Robertson & Goldreich 2012 for an example of a periodic box and V ≠ 0).

2.2.2. Accepted Class of Flows

In our homogeneous model, the bulk flow V is restricted to a certain class of flows. This class, however, contains many kinds of flows relevant to the present study, such as linearized shears, notably galactic shears, homogeneous rotations, and in particular solid rotations, and global homogeneous contractions or expansions, which need not be isotropic. We note that this construction is similar to that used in Newtonian cosmology, where usually V = H(t) x is the Hubble flow and H(t) is Hubble's expansion rate (see also Buchert & Ehlers 1997; Vigneron 2021 for the class of permitted flow in cosmology).

Therefore, these models can properly describe the evolution of the density-field statistics in star-forming clouds.

2.3. Ergodicity and the ACF of the Homogeneous Density Field

In ergodic theory, which specifies under which conditions the ergodic hypothesis is valid and provides an assessment of errors in the estimation of averages, the ACF C_ρ of the statistically homogeneous density field is of prime importance (see, e.g., E. Jaupart & G. Chabrier 2021, in preparation, hereafter JC21). It is defined as

$$\begin{eqnarray}&&{C}_{\rho }{({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} )\equiv {\mathbb{E}}\left(\rho ({\boldsymbol{x}})\rho ({\boldsymbol{x}}^{\prime} )\right)-\overline{\rho }(t)}^{2},\end{eqnarray} \tag{ 9 }$$

and reaches a maximum at ${\boldsymbol{\xi }}={\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} ={\bf{0}}$ : ${C}_{\rho }({\boldsymbol{\xi }})\leqslant {C}_{\rho }({\bf{0}})=\mathrm{Var}\left(\rho \right)$ (see, e.g., Papoulis & Pillai 1965 for a demonstration), where $\mathrm{Var}\left(\rho \right)$ is the variance of ρ. Length scales for which correlations are statistically significant are encoded in the ACF. This statistical object allows thus the extraction of characteristic length scales of physical processes.

2.3.1. Slutsky's Theorem and the Correlation Length

As mentioned above, one assumes statistical homogeneity and builds the following ergodic estimator for the expectation of ρ:

$$\begin{eqnarray}&&{\hat{\rho }}_{L}=\displaystyle \frac{1}{{L}^{3}}{\int }_{{\rm{\Omega }}}\rho ({\boldsymbol{x}})\,d{\boldsymbol{x}},\end{eqnarray} \tag{ 10 }$$

where ${\rm{\Omega }}={\left[-\tfrac{L}{2},\tfrac{L}{2}\right]}^{3}$ is a control volume of linear size L and volume L³, which is sought to be as large as possible. The ergodic estimator ${\hat{\rho }}_{L}$ has variance:

$$\begin{eqnarray}&&\mathrm{Var}({\hat{\rho }}_{L})=\displaystyle \frac{1}{{L}^{3}}{\int }_{2{\rm{\Omega }}}{C}_{\rho }({\boldsymbol{\xi }})\,\prod _{k=1}^{3}\left(1-\displaystyle \frac{| {\xi }_{k}| }{L}\right)\,d{\boldsymbol{\xi }},\end{eqnarray} \tag{ 11 }$$

where the integration volume 2Ω = [ − L, + L]³ stems from the change of variables $({\boldsymbol{x}},{\boldsymbol{x}}^{\prime} )\to ({\boldsymbol{\xi }}={\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} ,{\boldsymbol{y}}={\boldsymbol{x}}+{\boldsymbol{x}}^{\prime} )$. This leads to Slutsky's theorem (Papoulis & Pillai 1965): the stochastic field ρ is mean ergodic in the mean square (MS) sense, if and only if

$$\begin{eqnarray}&&\boxed{\displaystyle \frac{1}{{L}^{3}}{\int }_{2{\rm{\Omega }}}{C}_{\rho }({\boldsymbol{\xi }})d{\boldsymbol{\xi }}\mathop{\to }\limits_{L\to \infty }0.}\end{eqnarray} \tag{ 12 }$$

From this, one derives two sufficient (physical) conditions for ρ to be mean ergodic. Either

$$\begin{eqnarray}&&{\int }_{{{\mathbb{R}}}^{3}}{C}_{\rho }({\boldsymbol{\xi }})d{\boldsymbol{\xi }}\lt \infty ,\end{eqnarray} \tag{ 13 }$$

or

$$\begin{eqnarray}&&{C}_{\rho }({\boldsymbol{\xi }})\mathop{\longrightarrow }\limits_{| {\boldsymbol{\xi }}| \to \infty }0,\end{eqnarray} \tag{ 14 }$$

which means that values of the density field at two points separated by a lag ξ are uncorrelated at infinitely large distance. The first condition leads to the definition of the correlation length l_c (ρ) of the density field ρ (see, e.g., Papoulis & Pillai 1965):

$$\begin{eqnarray}&&\boxed{{({l}_{c}(\rho ))}^{3}=\displaystyle \frac{1}{{2}^{3}\,{C}_{\rho }({\bf{0}})}{\int }_{{{\mathbb{R}}}^{3}}{C}_{\rho }({\boldsymbol{\xi }})d{\boldsymbol{\xi }}=\displaystyle \frac{1}{{2}^{3}}{\int }_{{{\mathbb{R}}}^{3}}{\tilde{C}}_{\rho }({\boldsymbol{\xi }})d{\boldsymbol{\xi }},}\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{\tilde{C}}_{\rho }({\boldsymbol{\xi }})=\displaystyle \frac{{C}_{\rho }({\boldsymbol{\xi }})}{\mathrm{Var}\left(\rho \right)}\end{eqnarray} \tag{ 16 }$$

is the correlation coefficient at lag ξ that generates a measure of how correlated two values of the density field are. Then, using the two physical assumptions Equations (13) and (14), one obtains for l_c (ρ) ≪ L, and from Equation (11):

$$\begin{eqnarray}&&\mathrm{Var}({\hat{\rho }}_{L})\simeq \mathrm{Var}(\rho ){\left(\displaystyle \frac{2\,{l}_{c}(\rho )}{L}\right)}^{3}=\mathrm{Var}(\rho ){\left(\displaystyle \frac{{l}_{c}(\rho )}{R}\right)}^{3},\end{eqnarray} \tag{ 17 }$$

where R = L/2. Comparing Equation (17) with the variance $\mathrm{Var}\left({\hat{\rho }}_{{\boldsymbol{x}},N}\right)$ of the estimator of the expectation $\overline{\rho }$ obtained from a frequency interpretation where the experiment is repeated over N independent trials ω_i ,

$$\begin{eqnarray}&&{\hat{\rho }}_{{\boldsymbol{x}},N}=\displaystyle \frac{1}{N}\sum _{i=1}^{N}\rho ({\boldsymbol{x}},{\omega }_{i}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\mathrm{Var}({\hat{\rho }}_{{\boldsymbol{x}},N})=\displaystyle \frac{\mathrm{Var}(\rho )}{N},\end{eqnarray} \tag{ 19 }$$

we see that one can interpret the ratio ${(R/{l}_{c}(\rho ))}^{3}$ as an effective number of "independent" samples.

2.4. The Correlation Length as an Average Size of the Most Correlated Structures

2.4.1. Integral Scale

In homogeneous and isotropic turbulence, one introduces a quantity similar to the correlation length, called the integral scale l_i (not to be confused with the injection scale), defined as (Batchelor 1953)

$$\begin{eqnarray}&&{l}_{{\rm{i}}}(\rho )=\displaystyle \frac{1}{{C}_{\rho }(0)}{\int }_{0}^{\infty }{C}_{\rho }(\xi )d\xi ={\int }_{0}^{\infty }{\tilde{C}}_{\rho }(\xi )d\xi .\end{eqnarray} \tag{ 20 }$$

In the usual phenomenology of turbulence, this integral scale is used as a measure of the lags for which the velocities are significantly correlated and thus gives a measure of the accuracy of volumetric averages as estimates of actual statistical averages (Frisch 1995). The correlation length, the very quantity that enters Slutsky's theorem (Equation (12)), is given, in this isotropic context, by

$$\begin{eqnarray}&&{l}_{c}{\left(\rho \right)}^{3}=\displaystyle \frac{\pi }{2}{\int }_{0}^{\infty }{\tilde{C}}_{\rho }(\xi ){\xi }^{2}d\xi .\end{eqnarray} \tag{ 21 }$$

One finds that l_c ≃ l_i in many cases. Indeed, for an exponential ACF with ${C}_{\rho }({\boldsymbol{\xi }})=\mathrm{Var}(\rho ){e}^{-| {\boldsymbol{\xi }}| /{l}_{{\rm{i}}}}$, we have l_c = π^1/3 l_i, whereas for a Gaussian ACF (${C}_{\rho }({\boldsymbol{\xi }})=\mathrm{Var}(\rho ){e}^{-| {\boldsymbol{\xi }}{| }^{2}/\lambda }$), l_c = l_i . Moreover, for an ACF of the form ${C}_{\rho }{(\xi )=\mathrm{Var}(\rho )(1-(\xi /{l}_{0})}^{p})$ for r < l₀ and decaying rapidly outward, one gets l_c = (1.9 − 0.8) l_i for p ∈ [0.2, ∞ [ (the typical value in turbulence is p = 2/3 for the velocity field).

The integral scale can thus serve as a proxy for the correlation length, but this latter is the only quantity defined in absence of isotropy, as well as the one entering Slutsky's theorem (Equation (12)).

2.4.2. Average Size of the Most Correlated Structures

If the ACF of the ergodic field ρ is isotropic, the above equation for the integral scale l_i(ρ) can be used to define a weight function W_l (ξ) that measures the correlation of structures of size ξ = ∣ ξ ∣

$$\begin{eqnarray}&&{W}_{l}(r)=\displaystyle \frac{1}{{l}_{{\rm{i}}}(\rho )}\displaystyle \frac{{C}_{\rho }(r)}{\mathrm{Var}\left(\rho \right)}=\displaystyle \frac{{\tilde{C}}_{\rho }(r)}{{l}_{{\rm{i}}}(\rho )}.\end{eqnarray} \tag{ 22 }$$

Note that this weight function does not need to be positive and, in general, can have negative values, but its integral over all possible sizes ξ is 1 by construction. If the ACF of ρ is positive, however, W_l (r) can be further identified as the PDF of the size r of correlated structures. We can then build the weighted average of the size of correlated structures, $\left\langle {l}_{w}\right\rangle$, as :

$$\begin{eqnarray}&&\left\langle {l}_{w}\right\rangle ={\int }_{0}^{\infty }\xi {W}_{l}(\xi )d\xi =\displaystyle \frac{1}{{l}_{{\rm{i}}}(\rho )}{\int }_{0}^{\infty }{\tilde{C}}_{\rho }(\xi )\,\xi \,d\xi .\end{eqnarray} \tag{ 23 }$$

Then, as was the case for the integral scale, l_i(ρ), in many situations

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{\tilde{C}}_{\rho }{(\xi )\,\xi \,d\xi \simeq {l}_{c}(\rho )}^{2},\end{eqnarray} \tag{ 24 }$$

which yields:

$$\begin{eqnarray}&&\boxed{\left\langle {l}_{w}\right\rangle \simeq \displaystyle \frac{{l}_{c}{(\rho )}^{2}}{{l}_{{\rm{i}}}(\rho )}\simeq {l}_{c}(\rho ).}\end{eqnarray} \tag{ 25 }$$

Thus, l_c (ρ) measures the average size of correlated structures, weighted by the correlation coefficients ${\tilde{C}}_{\rho }(\xi )$. We then call this average size the average size of the most correlated structures, in order to indicate that it is a weighted average.

This construction, which relies on the assumption of isotropy, serves to illustrate the physical meaning of l_c (ρ). In the absence of such an assumption, l_c (ρ) is the only quantity that can be defined, but can still be interpreted as a measure of the average size of the most correlated structures. This is in agreement with the picture obtained from Equations (17) and (19), where the ratio ${(R/{l}_{c})}^{3}$ is interpreted as an effective number of "independent" samples in the volume V = (2R)³.

Chandrasekhar (1951a) derived a transport equation for the autocovariance function C_ρ in a statistically homogeneous isotropic and globally static medium with fixed background density $\overline{\rho }(t)={\rho }_{0}$. We generalize his result to our class of statistically homogeneous flows that are not necessarily isotropic and with non-trivial evolution, i.e., for which $\overline{\rho }(t)$ is a function of time and $\overline{{\boldsymbol{v}}}\ne {\bf{0}}$.

3.1. Transport Equation

The derivation of the transport equation follows the lines of Chandrasekhar (1951a) but accounting now for the non-trivial background flow; it is given in Appendix A. Expressing everything in terms of the logarithmic density s (see Equation (4)), we find:

$$\begin{eqnarray}\boxed{\begin{array}{rcl}0 & = & \displaystyle \frac{\partial }{\partial t}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)+{\left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)}^{i}\displaystyle \frac{\partial }{\partial {\xi }_{i}}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)\\ & & +\displaystyle \frac{\partial }{\partial {\xi }_{i}}{\left({R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}\right)}_{{\boldsymbol{\xi }}}+\displaystyle \frac{\partial }{\partial {\xi }_{i}}{\left({R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}\right)}_{-{\boldsymbol{\xi }}},\end{array}}\end{eqnarray} \tag{ 26 }$$

where ${R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}$ is the cross correlation function of the two fields e^s and e^s uⁱ , which depends only on the lag ξ under the assumption of statistical homogeneity. In fact, from the definition of u , $\overline{{e}^{s}{u}_{i}}=0$, so that ${R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}$ is also the cross covariance function of e^s and e^s uⁱ . If one assumes statistical isotropy,

$$\begin{eqnarray}&&{R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}({\boldsymbol{\xi }})={L}_{{e}^{s},{e}^{s}u}(| {\boldsymbol{\xi }}| )\,{\xi }^{i}.\end{eqnarray} \tag{ 27 }$$

Then, the last two terms on the right-hand side of Equation (26) can be combined to give $2{\partial }_{{\xi }_{i}}{\left({R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}\right)}_{{\boldsymbol{\xi }}}$ and we recover the result of Chandrasekhar (1951a) (his Equation (13)). Equation (26) thus generalizes the transport equation for the ACF of ρ derived by Chandrasekhar (1951a) for a non-isotropic, time-evolving flow:

$$\begin{eqnarray*}&&\displaystyle \frac{\partial }{\partial t}\left({C}_{\rho }({\boldsymbol{\xi }})\right)=2{\partial }_{{\xi }_{i}}\left({L}_{\rho ,\rho u}{\xi }^{i}\right),\end{eqnarray*}$$

with the addition of the advection term for relative velocity ${\rm{\Delta }}\overline{{v}^{i}}=\overline{{v}^{i}}({\boldsymbol{x}},t)-\overline{{v}^{i}}({\boldsymbol{x}}^{\prime} ,t)={\left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)}^{i}$, because distortion can only be generated by the relative motion (Kolmogorov 1941; Frisch 1995), and without assuming statistical isotropy at all scales.

As before, we use in the following the two common physical assumptions that enforce ergodicity. The covariance and cross covariance functions C_ρ (or ${C}_{{e}^{s}}$ ) and ${R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}$ are both assumed to decay rapidly to 0 as ∣ξ∣ → ∞ and to be integrable.

3.2. Correlation Length and Conserved Quantity

An important quantity characterizing the statistics of the stochastic field ρ (or e^s ) is the correlation length l_c (ρ) (or l_c (e^s )), defined earlier:

$$\begin{eqnarray}\begin{array}{rcl}{l}_{c}{\left(\rho \right)}^{3} & = & \displaystyle \frac{1}{8\mathrm{Var}\left(\rho \right)}{\displaystyle \int }_{{{\mathbb{R}}}^{3}}{C}_{\rho }({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\\ & = & \displaystyle \frac{1}{8\mathrm{Var}\left({e}^{s}\right)}{\displaystyle \int }_{{{\mathbb{R}}}^{3}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}={l}_{c}{\left({e}^{s}\right)}^{3}.\end{array}\end{eqnarray} \tag{ 28 }$$

Integrating Equation (26) over all possible lags ξ yields the conservation equation:

$$\begin{eqnarray}&&\boxed{{\left(\mathrm{Var}\left({e}^{s}\right){l}_{c}{\left({e}^{s}\right)}^{3}\right)}_{t}\overline{\rho }(t)=\mathrm{const}.}\end{eqnarray} \tag{ 29 }$$

or, in terms of the density field ρ:

$$\begin{eqnarray}&&\displaystyle \frac{{\left(\mathrm{Var}\left(\rho \right){l}_{c}{\left(\rho \right)}^{3}\right)}_{t}}{\overline{\rho }(t)}=\mathrm{const},\end{eqnarray} \tag{ 30 }$$

where quantities of the form ${\left(X\right)}_{t}$ mean that the value of quantity X is taken at time t. These two equations are modified versions of the conservation equation derived by Chandrasekhar (1951a) (his Equation (17)):

$$\begin{eqnarray*}&&{\int }_{0}^{\infty }{C}_{\rho }(r,t){r}^{2}{dr}=\mathrm{const}.\end{eqnarray*}$$

They account for evolution of the average (background) density field and depend explicitly on the correlation length. The detailed derivation of these equations is given in Appendix B.

We note that the conserved quantity in Equation (29) has the dimension of a mass; we will come back to this point later.

To test Equation (29) in astrophysical conditions, we use the numerical simulations presented in Federrath & Klessen (2012, 2013) and used in JC20. These simulations model the isothermal gravo-turbulent evolution of clouds in periodic boxes of size L with different resolutions ${N}_{\mathrm{res}}$, average density ρ₀, where turbulence is driven at fixed rms Mach numbers ${ \mathcal M }$ with solenoidal or compressive forcing or a mixture of both. They belong to the class of statistically homogeneous flows presented in Section 2.2.1 where V = 0. In each simulation, gravity is added after a gravitationless turbulence state has developed. As soon as gravity is switched on, the variance of the density field increases due to the condensation of structures. As shown above, the increase of the variance of ρ is expected to be accompanied by a decrease of the correlation length l_c (ρ). To measure this decrease we use the relation derived in JC21:

$$\begin{eqnarray}&&\mathrm{Var}\left(\displaystyle \frac{{\rm{\Sigma }}}{{\mathbb{E}}\left({\rm{\Sigma }}\right)}\right)\simeq \mathrm{Var}\left(\displaystyle \frac{\rho }{{\mathbb{E}}\left(\rho \right)}\right)\displaystyle \frac{{l}_{c}(\rho )}{R},\end{eqnarray} \tag{ 31 }$$

that relates the variance of the column-density field Σ to the variance of ρ, l_c (ρ) and the half size of the simulation box R = L/2. The derivation of Equation (31) is given in Appendix C. Equation (31) thus yields the estimate ${\hat{l}}_{c}/R$ of the ratio of the correlation length l_c (ρ) to the half size of the simulation box R = L/2:

$$\begin{eqnarray}&&\displaystyle \frac{{\hat{l}}_{c}}{R}\equiv \displaystyle \frac{\mathrm{Var}\left(\tfrac{{\rm{\Sigma }}}{{\mathbb{E}}\left({\rm{\Sigma }}\right)}\right)}{\mathrm{Var}\left(\tfrac{\rho }{{\mathbb{E}}\left(\rho \right)}\right)}\simeq \displaystyle \frac{{l}_{c}(\rho )}{R}.\end{eqnarray} \tag{ 32 }$$

As shown previously, had Equation (32) been an exact equality, we would expect ${\hat{l}}_{c}/R\propto \mathrm{Var}{\left(\rho \right)}^{-1/3}$ (for fixed $\overline{\rho }$). Equation (32), however, is only a proxy to derive an estimate of l_c (ρ) within a factor of order unity which depends on the shape of the ACF of the density field (see Section 2.4 and Appendix C). Furthermore, the ACF is initially that of inertial turbulence and evolves toward an ACF whose shape at short lags is determined by gravity-induced dynamics. We expect the ACF to change with time between these two regimes. Once the dynamics in high-density regions (short scales) start to be dominated by gravity (the regime we are interested in), we expect the ACF at short lags, while evolving with time, to preserve its functional form. In this regime, i.e., for $\mathrm{Var}{\left(\rho \right)}_{t}\gg \mathrm{Var}{\left(\rho \right)}_{t=0}$, we expect ${\hat{l}}_{c}/R\propto \mathrm{Var}{\left(\rho \right)}^{-1/3}$ (as mentioned above). However, as the simulations can only resolve structures larger than ${\rm{\Delta }}{x}_{\min }=L/{N}_{\mathrm{res}}$, resolution issues can prevent the occurrence of this behavior in the simulations. Instead, we expect values of ${\hat{l}}_{c}/R$ to level off at some point in the simulations.

Figure 1 displays estimated values of ${\hat{l}}_{c}/R$ as a function of the ratio $\mathrm{Var}{\left(\rho \right)}_{t}/\mathrm{Var}{\left(\rho \right)}_{t=0}$ (which increases with time) from hydrodynamic simulations for various Mach number ${ \mathcal M }$ and resolution ${N}_{\mathrm{res}}$. As expected, the ${\hat{l}}_{c}/R$ ratio decreases as the variance $\mathrm{Var}\left(\rho \right)$ increases. At high variance values (late times), the correlation length is observed to level off at a value that depends on the resolution ${N}_{\mathrm{res}}$, corresponding to ${\hat{l}}_{c}\simeq 2{\rm{\Delta }}{x}_{\min }={\lambda }_{J}({\rho }_{\max })/2$ where

$$\begin{eqnarray}&&{\lambda }_{J}({\rho }_{\max })=4{\rm{\Delta }}x,\end{eqnarray} \tag{ 33 }$$

with ${\rm{\Delta }}x=L/{N}_{\mathrm{res}}$ the grid resolution, is the Jeans length at density

$$\begin{eqnarray}&&{\rho }_{\max }=\displaystyle \frac{\pi {c}_{{\rm{s}}}^{2}}{16G{\rm{\Delta }}{x}_{\min }^{2}},\end{eqnarray} \tag{ 34 }$$

with c_s the sound speed, above which cloud collapsing features are not resolved (Truelove et al. 1997). For simulations at ${ \mathcal M }=50$ at the highest resolution ${N}_{\mathrm{res}}=1024$, we observe that the scaling ${\hat{l}}_{c}/R\propto \mathrm{Var}{\left(\rho \right)}^{-1/3}$ holds over a decade for $\mathrm{Var}{\left(\rho \right)}_{t}/\mathrm{Var}{\left(\rho \right)}_{t=0}\geqslant 5$. In the other simulations, this scaling law is inhibited by the leveling of ${\hat{l}}_{c}$ (save perhaps for ${ \mathcal M }=10$ where it holds for half a decade). It would thus be interesting to carry out all simulations with the same highest resolution (${N}_{\mathrm{res}}=1024$, for example).

Zoom In Zoom Out Reset image size

The initial values of ${\hat{l}}_{c}/R$ yield ${\hat{l}}_{c}/L=$ ${0.056}_{-0.013}^{+0.01}$, ${\hat{l}}_{c}/L=$ ${0.037}_{-0.006}^{+0.009}$, ${\hat{l}}_{c}/L=$ ${0.025}_{-0.003}^{+0.005}$, and ${\hat{l}}_{c}/L={0.013}_{-0.0018}^{+0.0015}$ for ${ \mathcal M }=3,$ 5, 10, 50, respectively. For simulations with ${ \mathcal M }\in \{3,5,10\}$, one finds that, within a factor of order unity, ${\hat{l}}_{c}\simeq L/{{ \mathcal M }}^{2}={\lambda }_{{\rm{s}}}$, where λ_s is the sonic length, which is found to be close to the average width of filamentary structures in isothermal turbulence (Federrath 2016). This is not surprising because l_c (ρ) describes the average size of the most correlated substructures. For the ${ \mathcal M }=50$ simulations, however, ${\hat{l}}_{c}$ is about 30 times larger than λ_s = L/2500. λ_s is not resolved in these simulations (${N}_{\mathrm{res}}=512$ or 1024), which explains the large discrepancy between ${\hat{l}}_{c}$ and its expected value λ_s.

The above results show that Equation (32) allows a good approximation of the actual ratio l_c (ρ)/R. They also emphasize the fact that correlated substructures are only resolved down to the smallest Jeans length that can be achieved in the simulations. They do not imply that structures larger than l_c (ρ) are not correlated. Such large correlated structures can exist (e.g., large filaments) but they are less correlated than the structures smaller than l_c (ρ) (i.e., they are associated with a lower correlation coefficient ${C}_{\rho }/\mathrm{Var}\left(\rho \right)$). Importantly enough, the simulations for the highest resolutions confirm that the quantity ${l}_{c}{\left(\rho \right)}^{3}\mathrm{Var}\left({e}^{s}\right)$ is indeed conserved, as expected from our theoretical analysis.

Observations and numerical simulations of MCs in dense star-forming regions have reported a significant increase of the density variance compared with the one obtained from gravitationless isothermal turbulence (e.g., Kainulainen et al. 2006, 2009; Schneider et al. 2012, 2013 and references therein for observations and, e.g., Kritsuk et al. 2010; Ballesteros-Paredes et al. 2011; Cho & Kim 2011; Collins et al. 2012; Federrath & Klessen 2013; Lee et al. 2015; Burkhart et al. 2016 for simulations). Such an increase is believed to be the signature of gravity.

5.1. Evolution of the Correlation Length

In JC20, we showed that this increase of variance due to gravity occurs on a short (local) timescale compared with the the typical timescale for variation of the cloud's global mean density $\overline{\rho }$.

This increase of the variance results in a decrease of the product $\overline{\rho }(t){l}_{c}{\left({e}^{s}\right)}^{3}$ in order to meet the constraint of the conservation equation (Equation (29)):

$$\begin{eqnarray}&&\mathrm{Var}\left({e}^{s}\right){\left(2{l}_{c}({e}^{s})\right)}^{3}\,\overline{\rho }(t)=\mathrm{const}.\end{eqnarray} \tag{ 35 }$$

Given the difference of timescales, we can assume, that, during this phase of variance increase, the (background) average density $\overline{\rho }=\mu {m}_{{\rm{H}}}\bar{n}$ (where μ and m_H = 1.66 × 10⁻²⁴ g denote the mean molecular weight and atomic mass unit, respectively) is almost constant and the conservation equation essentially holds

$$\begin{eqnarray}&&\boxed{\displaystyle \frac{{\left({l}_{c}({e}^{s})\right)}_{t}^{3}}{{\left({l}_{c}({e}^{s})\right)}_{t={t}_{0}}^{3}}\simeq \displaystyle \frac{{\left(\mathrm{Var}\left({e}^{s}\right)\right)}_{t={t}_{0}}}{{\left(\mathrm{Var}\left({e}^{s}\right)\right)}_{t}}\ll 1.}\end{eqnarray} \tag{ 36 }$$

It is worth stressing that whereas, by construction, $\bar{n}$ is exactly constant in mass-conserving simulations, it is not necessarily the case in real star-forming clouds, as it depends on the bulk flow (see Equation (6) and, e.g., Robertson & Goldreich 2012). Thus, the growing impact of gravity on the turbulent flow is accompanied by a drastic decrease of the correlation length of the density field l_c (e^s ) = l_c (ρ).

Physically speaking, Equation (36) implies that, during the cloud's evolution, the distribution of matter evolves from being concentrated in weakly correlated structures of average size ${\left({l}_{c}({e}^{s})\right)}_{t={t}_{0}}$ to being concentrated in smaller, denser more and more correlated regions of average size ${\left({l}_{c}({e}^{s})\right)}_{t}\ll {\left({l}_{c}({e}^{s})\right)}_{t={t}_{0}}$.

This picture is consistent with scenarios of star formation where the mass concentration in the cloud evolves from large filamentary structures to smaller, denser ones, and eventually to small prestellar cores (André 2017; André et al. 2019). Within the terminology of the present study, this is described as follows: dense and short-scale tightly correlated substructures (i.e., stellar cores) appear in larger, less correlated ones (i.e., filaments). The former ones correspond to objects of average size l_c (ρ)(t) whereas the latter correspond to objects of average (radial) size l_c (ρ)(t₀), which corresponds to the "initial" correlation length in early collapsing structures. Indeed, t₀ corresponds to the time at which some dense and significant regions within the cloud start to collapse and to deviate from the global evolution (contraction or expansion) of the cloud.

It is important to emphasize that the present theoretical framework, which is based on the hypothesis of statistical homogeneity, does not rely on any assumption regarding the condition or the magnitude of density deviation required for collapse. Furthermore, this framework is able to describe simultaneously a hierarchy of structures spanning a vast range of sizes and densities within the cloud during its evolution.

5.2. The Average Mass of Prestellar Cores

The quantity $\overline{\rho }(t)\mathrm{Var}\left({e}^{s}\right){\left(2{l}_{c}(\rho )\right)}^{3}$ has the dimension of a mass (Section 3.2) and corresponds to the average mass contained in the most correlated structures, M_corr:

$$\begin{eqnarray}&&{M}_{\mathrm{corr}}\propto \overline{\rho }\mathrm{Var}\left({e}^{s}\right){\left(2{l}_{c}(\rho )\right)}^{3},\end{eqnarray} \tag{ 37 }$$

with a proportionality coefficient of the order unity that depends on the geometry and where the 2³ term stems from the definition of l_c (ρ) since this latter corresponds to the half size of correlated structures. Initially, M_corr is located within the correlated structures embedded inside large filaments of average width l_c (ρ)(t₀). As collapse proceeds, this (conserved) amount of mass gets distributed in shorter scale, more correlated substructures of average size l_c (ρ)(t) < l_c (ρ)(t₀). Eventually, these structures will become prestellar cores. Thus, M_corr ultimately represents the average mass that is available to form (prestellar) cores. For a Chabrier-like core mass function (Chabrier 2003, 2005), this average mass is close to the characteristic mass. We calculate below an estimate of its value under typical Milky-Way-like conditions.

Observations and theoretical models of star formation indicate that initially, i.e., before the onset of star formation, the variance characteristic of the PDF of density fluctuations ressembles that of isothermal, fully developed turbulence (Padoan & Nordlund 2002; Mac Low & Klessen 2004; Hennebelle & Chabrier 2008; Hopkins 2012; Schneider et al.2013; De Oliveira et al. 2014; Vázquez-Semadeni et al. 2019), i.e.,

$$\begin{eqnarray}&&\mathrm{Var}\left({e}^{s}\right)({t}_{0})={\left(b{ \mathcal M }\right)}^{2}\end{eqnarray} \tag{ 38 }$$

(Federrath et al. 2008; Molina et al. 2012; Beattie et al. 2021). It remains to determine the correlation length l_c (ρ)(t₀). While, in case of pure gravitationless turbulence, this latter should be about the sonic length λ_s, it is not necessarily the case if gravity initially plays a non-negligible role. A detailed determination of the correlation length will be presented in a forthcoming paper.

Meanwhile, it is safe to take the observed average radial size of correlated filamentary structures, ∼0.1 pc, which is of the order of the sonic length, as an estimate of the "initial" correlation length l_c (ρ)(t₀) (see, e.g., Arzoumanian et al. 2011; Hennebelle & Falgarone 2012; Federrath 2016; Hennebelle & Inutsuka 2019 for a more complete discussion). Under typical Milky-Way-like conditions, this yields a characteristic correlated mass

$$\begin{eqnarray}&&\begin{array}{l}{M}_{\mathrm{corr}}={a}_{{\rm{g}}}\overline{\rho }\mathrm{Var}\left({e}^{s}\right){\left(2{l}_{c}(\rho )\right)}^{3}=1.56{M}_{\odot }\left(\displaystyle \frac{{a}_{{\rm{g}}}}{1}\right)\left(\displaystyle \frac{\mu }{2.0}\right)\\ \ \times \,\left(\displaystyle \frac{\overline{n}}{{10}^{3}\ {\mathrm{cc}}^{-1}}\right){\left(\displaystyle \frac{b\,{ \mathcal M }}{0.4\times 5}\right)}^{2}{\left(\displaystyle \frac{{l}_{c}(\rho )({t}_{0})}{0.1\ \mathrm{pc}}\right)}^{3}.\end{array}\end{eqnarray} \tag{ 39 }$$

It must be kept in mind that the last term of Equation (39), namely the correlation length, entails a dependence upon the cloud's initial main properties, average density, and Mach number (${l}_{c}{(\rho )({t}_{0})\equiv {l}_{c}(\rho )}_{0}[\bar{n},{ \mathcal M }]$).

It must be emphasized that M_corr is set up at the initial stages of gravitational contraction, before gravity starts significantly affecting the properties of turbulence (see JC20). It corresponds to the mass contained in the most correlated regions embedded in the initial filamentary structures generated by turbulence. It is not necessarily the average mass of all filaments in the cloud. Similarly, $\overline{n}$ is the initial average density of the cloud, representative of scales large enough that the ergodic estimates are accurate. It is not the average density in a (small-scale) collapsed subregion.

An important property of M_corr is that it is not expected to vary significantly among clouds which, initially, at large scale, meet the typical observed Larson conditions (i.e., $\overline{n}\sim {L}^{-{\eta }_{d}},{ \mathcal M }\sim {L}^{\eta }$, with η_d ∼ 0.7–1.0, η ∼ 0.4–0.5). Indeed, under such conditions, the quantity $\overline{n}{{ \mathcal M }}^{2}$, thus M_corr, remains approximately constant. This remarkable behavior has been advocated in a different approach, which involves a collapse criterion (namely the virial condition), to explain the apparent universality of the peak of the core mass function for a wide range of stellar cluster conditions (Hennebelle & Chabrier 2008).

As seen from Equation (39), the theory predicts that, under MW conditions, the average mass available to form prestellar cores, which is ultimately located in the most correlated structures of size l_c (ρ), is of the order of ∼ 1M_⊙, in agreement with observations (André et al. 2019).

The theory presented in this Letter, based solely on mass conservation in a statistically homogeneous medium (not necessarily isotropic or spatially homogeneous) with non-trivial evolution, first provides a description of the ACF and of the evolution of the correlation length l_c (ρ) of the density field in star-forming clouds. We show that this correlation length can be identified as the average size of the most correlated structures (see Section 2.4). Then, the theory provides a generalization of the transport equation derived by Chandrasekhar (1951a) for the ACF (Equation (26)) of density fluctuations in a turbulent medium. It demonstrates the occurrence of an invariant in the cloud's evolution, which is the average mass contained in the most correlated structures (Equation (29)). For any initial field of density fluctuations, this mass is conserved no matter what dominates the global dynamics (e.g., turbulence or gravity). Comparison with high-resolution numerical simulations (Federrath & Klessen 2012, 2013) confirms the theoretical relation (Section 4). This gives an original and robust description of the physical process occurring in star-forming clouds. As collapse progresses within (regions of) the cloud, the variance of the density field increases, therefore the correlation length l_c (ρ) decreases (Section 5.1), so collapse affects more and more correlated, shorter and shorter scales, yielding the formation of increasingly smaller and clumpier structures. Within this framework, dense and short-scale correlated substructures (cores) of average size l_c (ρ)(t) form in larger correlated structures (filaments) of average size l_c (ρ)(t₀) ∼ 0.1 pc. It is worth stressing that the theory, which is based on statistical homogeneity, does not constrain fluctuations around the average to be small and is able to simultaneously describe a hierarchy of structures spanning a large range of size and densities in various environments. The theory shows that, under Milky-Way-like typical conditions, the invariant average mass contained in the most correlated structures, which will eventually feed (prestellar) cores is of the order of ∼1 M_⊙, providing an appealing explanation for the universality of the peak of the IMF in MW environments.

Starting from the mass-conservation equation (Equation (1)) and multiplying it by $\rho ^{\prime} \equiv \rho ({\boldsymbol{x}}^{\prime} )$, one obtains:

$$\begin{eqnarray}&&\rho ^{\prime} \displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\rho ^{\prime} \rho \,{v}_{i}\right)=0.\end{eqnarray} \tag{ A1 }$$

Interchanging the primed and unprimed quantities in the above equation yields

$$\begin{eqnarray}&&\rho \displaystyle \frac{\partial \rho ^{\prime} }{\partial t}+\displaystyle \frac{\partial }{\partial x{{\prime} }_{i}}\left(\rho \,\rho ^{\prime} \,v{{\prime} }_{i}\right)=0.\end{eqnarray} \tag{ A2 }$$

Adding the two equations and taking the average, one obtains (Chandrasekhar 1951a):

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{R}_{\rho }({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} ,t)+\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\overline{\rho ^{\prime} \rho \,{v}_{i}}\right)+\displaystyle \frac{\partial }{\partial x{{\prime} }_{i}}\left(\overline{\rho \,\rho ^{\prime} \,v{{\prime} }_{i}}\right)=0.\end{eqnarray} \tag{ A3 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\rho }({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} ,t) & = & \overline{\rho ({\boldsymbol{x}})\rho ({\boldsymbol{x}}^{\prime} )},\\ & = & {C}_{\rho }({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} ,t)+\overline{\rho }{\left(t\right)}^{2},\end{array}\end{eqnarray} \tag{ A4 }$$

is the correlation function.

Decomposing V into the mean velocity V and turbulent component u ( v = V + u ), we obtain:

$$\begin{eqnarray}\begin{array}{rcl}0 & = & \displaystyle \frac{\partial }{\partial t}{C}_{\rho }({\boldsymbol{\xi }},t)-2\,{C}_{\rho }({\boldsymbol{\xi }})\displaystyle \frac{1}{\overline{\rho }}\displaystyle \frac{\partial }{\partial t}\overline{\rho }+({V}_{i}-V{{\prime} }_{i})\displaystyle \frac{\partial }{\partial {\xi }_{i}}{C}_{\rho }({\boldsymbol{\xi }})\\ & & +\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\overline{\rho ^{\prime} \rho \,{u}_{i}}\right)+\displaystyle \frac{\partial }{\partial x{{\prime} }_{i}}\left(\overline{\rho \,\rho ^{\prime} \,u{{\prime} }_{i}}\right),\end{array}\end{eqnarray} \tag{ A5 }$$

where $\xi ={\boldsymbol{x}}-{\boldsymbol{x}}^{\prime}$ and where we have used Equation (6). Then, dividing both sides by $\overline{\rho }{\left(t\right)}^{2}$ and using Equation (8), we obtain:

$$\begin{eqnarray}\begin{array}{rcl}0 & = & \displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{{C}_{\rho }({\boldsymbol{\xi }})}{{\overline{\rho }}^{2}}\right)+{\left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)}^{i}\displaystyle \frac{\partial }{\partial {\xi }_{i}}\left(\displaystyle \frac{{C}_{\rho }({\boldsymbol{\xi }})}{{\overline{\rho }}^{2}}\right)\\ & & +\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\overline{\displaystyle \frac{\rho ^{\prime} \rho }{{\overline{\rho }}^{2}}\,{u}_{i}}\right)+\displaystyle \frac{\partial }{\partial x{{\prime} }_{i}}\left(\overline{\displaystyle \frac{\rho \,\rho ^{\prime} }{{\overline{\rho }}^{2}}\,u{{\prime} }_{i}}\right).\end{array}\end{eqnarray} \tag{ A6 }$$

Expressing everything in terms of the logarithmic density s (see Equation (4)), we find:

$$\begin{eqnarray}\boxed{\begin{array}{rcl}0 & = & \displaystyle \frac{\partial }{\partial t}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)+{\left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)}^{i}\displaystyle \frac{\partial }{\partial {\xi }_{i}}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)\\ & & +\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\overline{{e}^{s^{\prime} }\,{e}^{s}\,{u}_{i}}\right)+\displaystyle \frac{\partial }{\partial x{{\prime} }_{i}}\left(\overline{{e}^{s}\,{e}^{s^{\prime} }\,u{{\prime} }_{i}}\right)\end{array}},\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & \displaystyle \frac{\partial }{\partial t}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)+{\left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)}^{i}\displaystyle \frac{\partial }{\partial {\xi }_{i}}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)\\ & & +\displaystyle \frac{\partial }{\partial {\xi }_{i}}{\left({R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}\right)}_{{\boldsymbol{\xi }}}+\displaystyle \frac{\partial }{\partial {\xi }_{i}}{\left({R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}\right)}_{-{\boldsymbol{\xi }}},\end{array}\end{eqnarray} \tag{ A8 }$$

where ${R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}$ is the cross correlation function of the two fields e^s and e^s uⁱ , which depends only on the lag ${\boldsymbol{\xi }}={\boldsymbol{x}}-{\boldsymbol{x}}^{\prime}$ under the assumption of statistical homogeneity.

To obtain the conserved quantity, Equation (29), one starts by noting that

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d}{{dt}}\left({\iiint }_{{{\mathbb{R}}}^{3}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\right) & = & {\iiint }_{{{\mathbb{R}}}^{3}}{\partial }_{t}{C}_{{e}^{s}}({\boldsymbol{\xi }},t)\,d{\boldsymbol{\xi }}\\ & = & -{\iiint }_{{{\mathbb{R}}}^{3}}{\left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)}^{i}\displaystyle \frac{\partial }{\partial {\xi }_{i}}\left({C}_{{e}^{s}}({\boldsymbol{\xi }})\right)\,d{\boldsymbol{\xi }}\\ & & -2{\iint }_{{\prime\prime} \partial {{\mathbb{R}}}^{3}{\prime\prime} }{R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}({\boldsymbol{\xi }})\,{{dS}}^{i},\end{array}\end{eqnarray} \tag{ B1 }$$

where the surface integral (the second term on the right-hand side of the equation) vanishes due to the assumption on ${R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}$. The first term on the right-hand side can be rewritten such that:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d}{{dt}}\left({\iiint }_{{{\mathbb{R}}}^{3}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\right) & = & -{\iiint }_{{{\mathbb{R}}}^{3}}{\boldsymbol{\nabla }}\cdot \left({C}_{{e}^{s}}({\boldsymbol{\xi }})\,\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)\,d{\boldsymbol{\xi }}\\ & & +{\iiint }_{{{\mathbb{R}}}^{3}}({\boldsymbol{\nabla }}\cdot {\boldsymbol{V}}){C}_{{e}^{s}}({\boldsymbol{\xi }})d{\boldsymbol{\xi }}.\end{array}\end{eqnarray} \tag{ B2 }$$

The first term on the right-hand side can be turned into a surface integral, which also vanishes due to the assumption on ${C}_{{e}^{s}}$. We are thus left with:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left({\iiint }_{{{\mathbb{R}}}^{3}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\right)=-\displaystyle \frac{d\mathrm{ln}(\overline{\rho })}{{dt}}\times \left(8\mathrm{Var}\left({e}^{s}\right){l}_{c}{\left({e}^{s}\right)}^{3}\right),\end{eqnarray} \tag{ B3 }$$

which yields:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left(\mathrm{Var}\left({e}^{s}\right){l}_{c}{\left({e}^{s}\right)}^{3}\right)=-\left(\mathrm{Var}\left({e}^{s}\right){l}_{c}{\left({e}^{s}\right)}^{3}\right)\times \displaystyle \frac{d\mathrm{ln}(\overline{\rho })}{{dt}},\end{eqnarray} \tag{ B4 }$$

and

$$\begin{eqnarray}&&\boxed{{\left(\mathrm{Var}\left({e}^{s}\right){l}_{c}{\left({e}^{s}\right)}^{3}\right)}_{t}\overline{\rho }(t)=\mathrm{const}.}\end{eqnarray} \tag{ B5 }$$

In principle, the integral in Equation (28) must be carried out over all possible lags ξ and hence over the whole space ${{\mathbb{R}}}^{3}$, which may seem conceptually problematic as we want to deal with a cloud of finite size. Regarding the bulk flow, however, we rely on the same line of reasoning as in statistical mechanics: if the actual subspace of permitted lags is large enough, it can be assimilated into the whole space ${{\mathbb{R}}}^{3}$. The argument is the following. If Ω, the subspace of permitted lags, is such that its volume ∣Ω∣ is $\gg {l}_{c}{\left({e}^{s}\right)}^{3}$, i.e., contains a large number of correlation volumes, and if C_ρ (or ${C}_{{e}^{s}}$ ) tends to 0 as ∣ξ∣ → ∞ and is integrable, the integral over Ω can be seen as an integral over ${{\mathbb{R}}}^{3}$.

To understand the meaning of the conserved quantity in Equation (29) (or Equation (B5)) and the approximation made, we now consider a finite subspace of permitted lags. Let Ω_t be the "average" volume of space describing the cloud under study, evolving with the average velocity field $\overline{{\boldsymbol{v}}}={\boldsymbol{V}}({\boldsymbol{x}},t)+\overline{{\boldsymbol{u}}}(t)=\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{x}}+{{\boldsymbol{c}}}_{V}(t)+\overline{{\boldsymbol{u}}}(t)$. Ω_t is hence a mass-conserving domain and, like $\overline{\rho }(t)$, is allowed to evolve with time. If Ω_t possesses point symmetry, then the subspace of permitted lags is simply Ω_t,ξ = 2 Ω_t . This subspace is evolving with the relative velocity field ${\rm{\Delta }}\overline{{\boldsymbol{v}}}=\overline{{\boldsymbol{v}}}$ $({\boldsymbol{x}},t)-\overline{{\boldsymbol{v}}}$ $({\boldsymbol{x}}^{\prime} ,t)=\underline{\underline{{L}_{V}}}$ (t) · ξ , because distorsion can only be generated by the relative motion (Kolmogorov 1941; Frisch 1995). Due to the Reynolds' transport theorem, one has:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\right)\\ \quad =-\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}{| }^{2}}\displaystyle \frac{d| {{\rm{\Omega }}}_{t}| }{{dt}}{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}+\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\times \left({\boldsymbol{\nabla }}\cdot \,{\rm{\Delta }}\overline{{\boldsymbol{v}}}\right)\,d{\boldsymbol{\xi }}\\ \quad +\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}\left({\partial }_{t}+(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }})\cdot {\boldsymbol{\nabla }}\right){C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }},\end{array}\end{eqnarray} \tag{ B6 }$$

where:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \,{\rm{\Delta }}\overline{{\boldsymbol{v}}}={\boldsymbol{\nabla }}\cdot \left(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }}\right)={\boldsymbol{\nabla }}\cdot \,{\boldsymbol{V}}=-\displaystyle \frac{d\mathrm{ln}(\overline{\rho })}{{dt}},\end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }\displaystyle \frac{d| {{\rm{\Omega }}}_{t}| }{{dt}}=-\displaystyle \frac{d\mathrm{ln}(\overline{\rho })}{{dt}}={\boldsymbol{\nabla }}\cdot \,{\boldsymbol{V}}.\end{eqnarray} \tag{ B8 }$$

This leads to:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\right)\\ \quad =\,\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}\left({\partial }_{t}+(\underline{\underline{{L}_{V}}}(t)\cdot {\boldsymbol{\xi }})\cdot {\boldsymbol{\nabla }}\right){C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\\ \quad =\,-2\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iint }_{\partial {{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}{R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}({\boldsymbol{\xi }})\,{{dS}}^{i}.\end{array}\end{eqnarray} \tag{ B9 }$$

Assuming now that the contribution from the surface integral at the rhs of Equation (B9) is negligible, i.e., assuming that ${R}_{{e}^{s},{e}^{s}{\boldsymbol{u}}}^{i}$ decays rapidly to 0 at large lags ξ and that Ω_t (and hence Ω_{t, ξ} ) is large enough (for example such that $| {{\rm{\Omega }}}_{t}| \gg {l}_{c}{\left({e}^{s}\right)}^{3}$), we are left with:

$$\begin{eqnarray}&&\displaystyle \frac{1}{| {{\rm{\Omega }}}_{t}| }{\iiint }_{{{\rm{\Omega }}}_{t,{\boldsymbol{\xi }}}}{C}_{{e}^{s}}({\boldsymbol{\xi }})\,d{\boldsymbol{\xi }}\simeq 8\,\mathrm{Var}\left({e}^{s}\right)\displaystyle \frac{{l}_{c}{\left({e}^{s}\right)}^{3}}{| {{\rm{\Omega }}}_{t}| }=\mathrm{const}.\end{eqnarray} \tag{ B10 }$$

Using the fact that $\overline{\rho }(t)| {{\rm{\Omega }}}_{t}| =M({{\rm{\Omega }}}_{t})=\mathrm{const}$, we obtain:

$$\begin{eqnarray}&&\boxed{\mathrm{Var}\left({e}^{s}\right)\,{l}_{c}{\left({e}^{s}\right)}^{3}\,\overline{\rho }(t)=\mathrm{const},}\end{eqnarray} \tag{ B11 }$$

which is Equation (B5) (or Equation (29)). These calculations are valid for any (mass-conserving) sub-domain Ω_t that is large enough for the surface integral on the rhs of Equation (B9) to be negligible. Equation (B11) therefore implies that the fundamental quantity $\mathrm{Var}\left({e}^{s}\right)\,{l}_{c}{\left({e}^{s}\right)}^{3}\,\overline{\rho }(t)$ is conserved.

We give the derivation of Equation (C7). For a cubic simulation domain of size L, projecting the density field along one of the three principal directions of the cube leads to a statistically homogeneous column-density field such that:

$$\begin{eqnarray}&&{\mathbb{E}}\left({\rm{\Sigma }}(x,y)\right)={\mathbb{E}}\left(\rho \right)\times L.\end{eqnarray} \tag{ C1 }$$

In a cubic box, the ACF of Σ is

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\rm{\Sigma }}}({\boldsymbol{r}}) & = & {\mathbb{E}}\left(\left({\rm{\Sigma }}({\boldsymbol{u}}+{\boldsymbol{r}})-{\mathbb{E}}\left(\rho \right)\,L\right)\left({\rm{\Sigma }}({\boldsymbol{u}})-{\mathbb{E}}\left(\rho \right)\,L\right)\right)\\ & = & {\displaystyle \int }_{{[-L/2,L/2]}^{2}}{C}_{\rho }({\boldsymbol{r}},z-z^{\prime} ){dz}\,{dz}^{\prime} \\ & = & {\displaystyle \int }_{[-L,L]}{C}_{\rho }({\boldsymbol{r}},u){du}{\displaystyle \int }_{-L+| u| }^{L-| u| }\displaystyle \frac{{dv}}{2}\\ & = & L{\displaystyle \int }_{[-L,L]}{C}_{\rho }({\boldsymbol{r}},u)\left(1-\displaystyle \frac{| u| }{L}\right){du},\end{array}\end{eqnarray} \tag{ C2 }$$

while the variance is

$$\begin{eqnarray}&&\mathrm{Var}({\rm{\Sigma }})={C}_{{\rm{\Sigma }}}({\bf{0}})=L{\displaystyle \int }_{[-L,L]}{C}_{\rho }({\bf{0}},u)\left(1-\displaystyle \frac{| u| }{L}\right){du}.\end{eqnarray} \tag{ C3 }$$

Thus, assuming that the density field is statistically isotropic at small scales (i.e., the ACF is isotropic at short lags), one obtains:

$$\begin{eqnarray}&&\mathrm{Var}({\rm{\Sigma }})\simeq L{\int }_{[-L,L]}{C}_{\rho }(| u| )\left(1-\displaystyle \frac{| u| }{L}\right){du}.\end{eqnarray} \tag{ C4 }$$

Provided the correlation length of the density field is much smaller than the size of the box L (i. e., l_c (ρ) ≪ L), one can approximate the integral on the rhs of Equation (C4) by the following expression:

$$\begin{eqnarray}&&{\int }_{[-L,L]}{C}_{\rho }(| u| )\left(1-\displaystyle \frac{| u| }{L}\right){du}\simeq 2{l}_{{\rm{i}}}(\rho )\mathrm{Var}(\rho )\simeq 2{l}_{c}(\rho )\mathrm{Var}(\rho ),\end{eqnarray} \tag{ C5 }$$

where l_i(ρ) is the integral scale of the density field. Thus,

$$\begin{eqnarray}&&\mathrm{Var}({\rm{\Sigma }})\simeq 2\,L\,{l}_{c}(\rho )\mathrm{Var}(\rho ).\end{eqnarray} \tag{ C6 }$$

This yields

$$\begin{eqnarray}&&\boxed{\mathrm{Var}\left(\displaystyle \frac{{\rm{\Sigma }}}{{\mathbb{E}}\left({\rm{\Sigma }}\right)}\right)\simeq \mathrm{Var}\left(\displaystyle \frac{\rho }{{\mathbb{E}}\left(\rho \right)}\right)\displaystyle \frac{2\,{l}_{c}(\rho )}{L}=\mathrm{Var}\left(\displaystyle \frac{\rho }{{\mathbb{E}}\left(\rho \right)}\right)\displaystyle \frac{{l}_{c}(\rho )}{R},}\end{eqnarray} \tag{ C7 }$$

where R = L/2. This is an important result because it provides a measure of l_c (ρ)/R independently of the ACF.