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

HAL
open science

# Evolution of the Density PDF in Star-forming Clouds:
The Role of Gravity

Etienne Jaupart, Gilles Chabrier

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

► To cite this version:

Etienne Jaupart, Gilles Chabrier. Evolution of the Density PDF in Star-forming Clouds: The Role of Gravity. The Astrophysical Journal Letters, 2020, 903, 10.3847/2041-8213/abbda8. insu-03711455

HAL Id: insu-03711455

https://insu.hal.science/insu-03711455v1

Submitted on 1 Aug 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

# Evolution of the Density PDF in Star-forming Clouds: The Role of Gravity

Etienne Jaupart¹ and Gilles Chabrier¹,²

¹Ecole normale supérieure de Lyon, CRAL, Université de Lyon, UMR CNRS 5574, F-69364 Lyon Cedex 07, France; etienne.jaupart@ens-lyon.fr, chabrier@ens-lyon.fr

²School of Physics, University of Exeter, Exeter EX4 4QL, UK

Received 2020 June 3; revised 2020 September 15; accepted 2020 October 1; published 2020 October 23

## Abstract

We derive an analytical theory of the PDF of density fluctuations in supersonic turbulence in the presence of gravity in star-forming clouds. The theory is based on a rigorous derivation of a combination of the Navier–Stokes continuity equations for the fluid motions and the Poisson equation for the gravity. It extends previous approaches, first by including gravity and second by considering the PDF as a dynamical system, not a stationary one. We derive the transport equations of the density PDF, characterize its evolution, and determine the density threshold above which gravity strongly affects and eventually dominates the dynamics of turbulence. We demonstrate the occurrence of two power-law tails in the PDF, with two characteristic exponents, corresponding to two different stages in the balance between turbulence and gravity. Another important result of this study is to provide a procedure to relate the observed column density PDFs to the corresponding volume density PDFs. This allows us to infer, from the observation of column densities, various physical parameters characterizing molecular clouds, notably the virial parameter. Furthermore, the theory offers the possibility to date the clouds in units of $t_{\text{coll}}$, the time since a statistically significant fraction of the cloud started to collapse. The theoretical results and diagnostics reproduce very well numerical simulations and observations of star-forming clouds. The theory provides a sound theoretical foundation and quantitative diagnostics to analyze observations or numerical simulations of star-forming regions and to characterize the evolution of the density PDF in various regions of molecular clouds.

Unified Astronomy Thesaurus concepts: Molecular clouds (1072); Star formation (1569); Hydrodynamics (1963)

## 1. Introduction

It has been established by many studies that the volume-weighted probability density function (PDF) of supersonic isothermal turbulence displays a nearly log-normal shape for solenoidically driven turbulence, at least for Mach numbers $M \lesssim 30$ (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Kritsuk et al. 2007; Federrath et al. 2008, 2010; Pan et al. 2019a), even in the presence of magnetic fields (Lemaster & Stone 2008; Collins et al. 2012). In dense star-forming regions, however, the line-of-sight extinction and inferred column density PDFs have been observed to develop a power-law tail at high densities, for extinctions $A_V \gtrsim 2-5$ (e.g., Kainulainen et al. 2006, 2009; Schneider et al. 2012, 2013 and references therein), a feature identified as the signature of gravity. Indeed, a similar feature of the PDF is found in numerical simulations of turbulence that include self-gravity (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).

Whereas these two opposite signatures of the PDF, log-normal versus power law, seem to be clearly identified, when and how precisely gravity starts affecting the dynamics of turbulence and thus the properties and the evolution of the PDF remains to be fully understood. Understanding the statistical properties of supersonic turbulence is at the heart of analytical theories of the star formation process, so understanding the physics governing the shape and the evolution of the density PDF of supersonic turbulence is of prime importance. A few attempts have been made to explain the development of power-law tails in the density PDF (Girichidis et al. 2014; Donkov & Stefanov 2018; Guszejnov et al. 2018). These approaches, however, focus only on the gravitationally unstable parts of a cloud, using self-similar gravitational collapse models and/or

geometrical arguments. While these models derive asymptotic exponents of power-law tails, they lack a complete description of the density fluctuations, in both the gravitationally stable and unstable parts of the cloud, and treat the PDF as a static system, even though Girichidis et al. (2014) follow its time evolution numerically. A first attempt to derive a robust theoretical framework of density fluctuations in compressible turbulence has been addressed by Pan et al. (2018, 2019a) based on the formalism developed by Pope (1981, 1985) and Pope & Ching (1993) for the PDF of any quantity, expressed as the conditional expectations of its time derivatives. Within the framework of this formalism, Pan et al. (2019a) used a probabilistic approach to turbulence to derive a theoretical formulation of the PDF of density fluctuations in a steady state from first principles. In this paper, we follow a similar approach and derive an analytical theory to describe the dynamics of turbulence in dense regions of molecular clouds (MCs) and its interplay with gravity. The approach generalizes the aforementioned ones in two ways. First, we include the impact of gravity on the cloud’s dynamics. Second, we consider the density PDF not as a stationary system but as one evolving with time, implying that the conditional expectation of the flow velocity divergence is time-dependent and nonzero. The theory explains the evolution of the PDF and determines the density thresholds above which gravity strongly affects and eventually dominates the dynamics of the turbulence. We provide a procedure to relate the observed column density PDFs to the underlying volume density PDFs, allowing us to infer various physical parameters characterizing molecular clouds from observations. The theory and its diagnostics are confronted with numerical simulations of gravoturbulent collapsing clouds and with various available observations.
## 2. Mathematical Framework

### 2.1. Description of a Molecular Cloud

We consider an isolated, turbulent, self-gravitating molecular cloud. Neglecting for now the magnetic field, the cloud's evolution is given by the standard Navier-Stokes and Poisson equations:

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

$$ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{1}{\rho} \nabla P + \mathcal{G} + \frac{1}{\rho} \nabla \cdot \underline{\underline{\sigma}}_{\nu}, \qquad (2) $$

$$ \nabla \cdot \mathcal{G} = -4\pi G \rho, \qquad (3) $$

where $\rho$ and $P$ denote respectively the density and pressure of the gas in the cloud, $\mathbf{v}$ the velocity field, $\mathcal{G}$ the gravity field, and $\underline{\underline{\sigma}}_{\nu}$ the viscous stress tensor. We close the system of equations by using a barotropic equation of state $P(\rho) \propto \rho^{\gamma}$ for the gas.

We separate the evolution of the background from that of local density deviations. We thus split the velocity field $\mathbf{v}$ between the mean velocity $\mathbf{V}$ and the (turbulent) velocity $\mathbf{u}$ (Ledoux & Walraven 1958) and we introduce the logarithmic excess of local density $s$,

$$ \mathbf{V} = \frac{1}{\bar{\rho}} \bar{\rho} \mathbf{v}, \qquad (4) $$

$$ \mathbf{u} = \mathbf{v} - \mathbf{V}, \qquad (5) $$

$$ \rho = \bar{\rho}(\mathbf{x}, t) e^s, \qquad (6) $$

where $\bar{\Phi}(\mathbf{x}, t)$ denotes the mathematical expectation, also called statistical average or mean, of any random field $\Phi$ (e.g., Pope 1985; Frisch 1995). We note that $\bar{\mathbf{u}} \neq 0$ a priori but $\bar{\rho}\bar{\mathbf{u}} \equiv 0$ by definition (Equations (4) and (5)). 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, i.e., for $\bar{\rho}$ and $\mathbf{V}$. Averaging Equations (1) and (2) yields an evolution equation for the mean flow written in conservative form,

$$ \frac{\partial(\bar{\rho} \mathbf{V})}{\partial t} + \nabla \cdot (\bar{\rho} \mathbf{V} \otimes \mathbf{V}) = - \nabla \bar{P} + \bar{\rho} \underline{\underline{\mathcal{G}}} + \nabla \cdot (\underline{\underline{\sigma}}_{\nu} - \bar{\rho} \mathbf{u} \otimes \mathbf{u}), \qquad (7) $$

with all quantities replaced by their mean values, except for the appearance of the turbulent Reynolds stress tensor, $-\bar{\rho}\mathbf{u} \otimes \mathbf{u}$. The trace of this tensor corresponds to the turbulent pressure while its traceless part is related to the turbulent viscous tensor. We consider molecular viscous effects to be negligible in molecular clouds and thus neglect the viscous tensor $\underline{\underline{\sigma}}_{\nu}$ in the general equations.

Having obtained the averaged evolution of the cloud, we can obtain the evolution of the density deviations by subtraction. This yields a transport equation for $s$, written in a Lagrangian form:

$$ \frac{Ds}{Dt} = - \nabla \cdot \mathbf{u} - (\mathbf{u} \cdot \nabla) \ln(\bar{\rho}), \qquad (8) $$

where $\frac{D}{Dt} = \frac{\partial}{\partial t} + (\mathbf{v} \cdot \nabla)$ is the Lagrangian derivative.

### 2.2. Transport Equations for the Probability Distribution Function of Logarithmic Density Fluctuations

Assuming that turbulent fields are statistically homogeneous, one can derive two transport equations for the probability distribution function of logarithmic density fluctuations (s-PDF) $f$ (Pope 1981, 1985; Pan et al. 2018, 2019a, 2019b):

$$ \frac{\partial}{\partial t} f(s, t) = \left\{ 1 + \frac{\partial}{\partial s} \right\} \left[ \langle (\nabla \cdot \mathbf{u}) | s \rangle f \right], \qquad (9) $$

$$ \frac{\partial}{\partial t} [\langle (\nabla \cdot \mathbf{u}) | s \rangle f] = \left\{ 1 + \frac{\partial}{\partial s} \right\} [\langle (\nabla \cdot \mathbf{u})^2 | s \rangle f] + f \left\langle \frac{D \nabla \cdot \mathbf{u}}{Dt} | s \right\rangle, \qquad (10) $$

where terms of the form $\langle \Phi | s \rangle \equiv \langle \Phi | s(\mathbf{x}, t) = s \rangle$ denote the conditional expectations of the random field $\Phi$ knowing that $s(\mathbf{x}, t) = s$, and can be computed as the average of the field $\Phi$ in all regions where $s(\mathbf{x}, t) \in [s, s + ds]$.

### 2.3. Stationary Solutions

Equations (9) and (10) give insights into the interplay between dynamical quantities and the steady-state value of the density s-PDF, $f$. Pan et al. (2018, 2019a, 2019b) have shown and tested numerically that

1. $f$ is stationary if and only if $\langle \nabla \cdot \mathbf{u} | s \rangle = 0, \forall s$,

2. In a steady state, $f$ can be formally computed as

$$ f(s) = \frac{Ce^{-s}}{\langle (\nabla \cdot \mathbf{u})^2 | s \rangle} \exp \left( - \int_0^s \frac{\left\langle \frac{D \nabla \cdot \mathbf{u}}{Dt} \right\rangle | s' \rangle}{\langle (\nabla \cdot \mathbf{u})^2 | s' \rangle} ds' \right), \qquad (11) $$

enabling us to discuss the impact of dynamical effects on the density PDF $f(s)$.

### 2.4. Effects of Gravity on the Density PDF

The effect of gravity on the density PDF without assuming a steady state can be inferred by recasting Equation (10) as an equation for $\ln f$:

$$ \begin{aligned} \langle \nabla \cdot \mathbf{u} | s \rangle \frac{\partial}{\partial t} \ln f - \langle (\nabla \cdot \mathbf{u})^2 | s \rangle \frac{\partial}{\partial s} \ln f + \frac{\partial}{\partial t} \langle \nabla \cdot \mathbf{u} | s \rangle \\ = \left\langle \frac{D \nabla \cdot \mathbf{u}}{Dt} | s \right\rangle + \left\{ 1 + \frac{\partial}{\partial s} \right\} \langle (\nabla \cdot \mathbf{u})^2 | s \rangle, \end{aligned} \qquad (12) $$

where the terms on the right-hand side (rhs) are then treated as source terms. We note that, due to Equation (9), the term $\partial t \langle \nabla \cdot \mathbf{u} | s \rangle$ on the left-hand side of Equation (12) is in fact seen as an operator acting on $f$; in a similar way the pressure gradient is seen as a nonlocal operator acting on the velocity field in standard studies of incompressible hydrodynamics with periodic boundary conditions (see, e.g., Frisch 1995). We then split $\langle \frac{D \nabla \cdot \mathbf{u}}{Dt} | s \rangle$ as

$$ \left\langle \frac{D \nabla \cdot \mathbf{u}}{Dt} | s \right\rangle = S_{\text{turb}}(s, t) + S_{\text{grav}}(s, t) + S_{\text{th}}(s, t) \qquad (13) $$

with

$$ S_{\text{grav}}(s, t) \equiv -4\pi G \bar{\rho}(e^s - 1), \qquad (14) $$
$$S_{\text{th}}(s, t) \equiv - \left\langle \nabla \cdot \left( \frac{1}{\rho} \nabla P \right) |s \right\rangle, \quad (15)$$

$$S_{\text{turb}}(s, t) \equiv - \langle \nabla \mathbf{u}: \nabla \mathbf{u} | s \rangle - 2 \langle \nabla \mathbf{V}: \nabla \mathbf{u} | s \rangle, \quad (16)$$

where $\nabla x: \nabla y = (\partial_i x_j) (\partial_j y_i)$ using Einstein's summation convention. Equations (14)–(16) are obtained by taking the divergence of Equation (2) and subtracting its average, knowing that the turbulent fields $\rho$ and $\mathbf{u}$ are statistically homogeneous. Their explicit derivation is given in Appendix A. Then, using Equation (12), the statistics of the flow within the cloud will be dominated by gravity (i.e., will differ from the statistics of pure gravitationless turbulence), whenever

$$|S_{\text{grav}}(s)| \gtrsim \max(|S_{\text{turb}}|, |S_{\text{th}}|, |\{1 + \partial_s\} \langle (\nabla \cdot \mathbf{u})^2 |s\rangle|). \quad (17)$$

Note that if the dynamics is dominated by gravity, we expect $\langle (\nabla \cdot \mathbf{u}) |s\rangle$ to be amplified in collapsing regions such that $|S_{\text{grav}}(s)| \sim \langle (\nabla \cdot \mathbf{u})^2 |s\rangle$ (see Section 3.3).

Physically, Equation (17) expresses the fact that gravity dominates whenever one of the two following conditions is fulfilled. Either (1) it overcomes thermal (pressure) or turbulent contributions to the dynamics of the cloud ($|S_{\text{grav}}(s)| \gtrsim \max(|S_{\text{turb}}|, |S_{\text{th}}|$); or (2) either convergent flows are produced by gravitational collapse ($|S_{\text{grav}}(s)| \sim \langle (\nabla \cdot \mathbf{u})^2 |s\rangle$) or divergent flows are forced to collapse, regardless of their their initial expansion ($S_{\text{grav}}(s) > \langle (\nabla \cdot \mathbf{u})^2 |s\rangle$).

As the aim of our study is to know when gravity will yield significant departures from pure (gravitationless) turbulence, we can evaluate the terms on the rhs of Equation (17) as for standard steady-state turbulence without gravity (which we denote with the subscript $\mathcal{G}$):

$$|S_{\text{grav}}(s)| \gtrsim \max(|S_{\text{turb}}|, |S_{\text{th}}|, |\{1 + \partial_s\} \langle (\nabla \cdot \mathbf{u})^2 |s\rangle|)_{\mathcal{G}}. \quad (18)$$

Pan et al. (2019a) performed such an analysis and found that $\langle (\nabla \cdot \mathbf{u})^2 |s\rangle_{\mathcal{G}} \sim (\nabla \cdot \mathbf{u})^2_{\mathcal{G}}$, while the other terms have no straightforward functional forms. To further simplify Equation (18), we start from Equation (12) for turbulence without gravity,

$$-\left(\langle (\nabla \cdot \mathbf{u})^2 |s\rangle \frac{\partial}{\partial s} \ln f\right)_{\mathcal{G}} = (S_{\text{turb}})_{\mathcal{G}} + (S_{\text{th}})_{\mathcal{G}} + \langle (1 + \partial_s) \langle (\nabla \cdot \mathbf{u})^2 |s\rangle \rangle_{\mathcal{G}}, \quad (19)$$

then use of the triangle inequality,

$$\left| \langle (\nabla \cdot \mathbf{u})^2 |s\rangle \frac{\partial}{\partial s} \ln f \right|_{\mathcal{G}} \le |S_{\text{turb}}|_{\mathcal{G}} + |S_{\text{th}}|_{\mathcal{G}} + \langle (1 + \partial_s) \langle (\nabla \cdot \mathbf{u})^2 |s\rangle \rangle_{\mathcal{G}}, \quad (20)$$

yields for the condition given by Equation (18)

$$|S_{\text{grav}}(s)| \gtrsim \left| \langle (\nabla \cdot \mathbf{u})^2 |s\rangle \frac{\partial}{\partial s} \ln f \right|_{\mathcal{G}}. \quad (21)$$

Making the standard approximation that $f_{\mathcal{G}}$ is a log-normal form of variance $\sigma_s$ yields the simplified condition

$$|S_{\text{grav}}(s)| \gtrsim \overline{(\nabla \cdot \mathbf{u})^2}_{\mathcal{G}} \times \left| \frac{s + \frac{1}{2} \sigma_s^2}{\sigma_s^2} \right|, \quad (22)$$

where $\sigma_s$ is given in terms of the rms Mach number $M$ and forcing parameter $b$ as (e.g., Federrath et al. 2008)

$$\sigma_s^2 = \ln(1 + (bM)^2). \quad (23)$$

Then, approximating $(\nabla \cdot \mathbf{u})^2_{\mathcal{G}}$ within an order-of-magnitude estimate as

$$\overline{(\nabla \cdot \mathbf{u})^2}_{\mathcal{G}} = \frac{1}{\bar{\rho}^2} \frac{(\Delta \rho)^2}{\tau_{\text{turb}}^2}, \quad (24)$$

where $\overline{(\Delta \rho)^2} = \overline{(\rho - \bar{\rho})^2} \simeq (bM\bar{\rho})^2$ and $\tau_{\text{turb}}$ is a typical turbulent timescale, of the order of the crossing time $\tau_c = L_c/(2\sigma_v)$, with $\sigma_v$ the 3D velocity dispersion and $L_c$ the diameter of the cloud, Equation (22) reduces to

$$\begin{aligned} |e^s - 1| &\gtrsim (bM)^2 \times \left(\frac{\tau_{G,0}}{\tau_{\text{turb}}}\right)^2 \times \left| \frac{s + \frac{1}{2}\sigma_s^2}{\sigma_s^2} \right| \\ &\gtrsim (bM)^2 \times \alpha_{\text{vir}}(t) \times \left| \frac{s + \frac{1}{2}\sigma_s^2}{\sigma_s^2} \right|, \end{aligned} \quad (25)$$

where $\tau_{G,0} = 1/\sqrt{4\pi G\bar{\rho}}$ and $\alpha_{\text{vir}}(t) = 5\sigma_v^2/(\pi G L_c^2 \bar{\rho}(t))$ is the virial parameter, equal to $= 2E_{\text{kin}}/E_{\text{grav}}$ for a homogeneous spherical cloud.

This equation introduces a new characteristic timescale, $\tau_{G,0} \equiv 1/\sqrt{4\pi G\bar{\rho}}$. This timescale characterizes the impact of gravity upon turbulence in the PDF evolution of the cloud, as formalized by Equations (13) and (25). It is roughly half the mean freefall time of the cloud, $\tau_{\text{ff},0} \equiv \sqrt{\frac{3\pi}{32G\bar{\rho}}}$. Equation (25) then allows a determination, within a factor of a few, of the density above which gravity is expected to change the statistics of turbulence significantly.

Furthermore, following Pan et al. (2019a) and using Equation (11), we see that because $S_{\text{grav}}(s) < 0$ when $s > 0$ and $S_{\text{grav}}(s) > 0$ when $s < 0$, respectively, gravity tends to broaden the PDF at both small and large densities, resulting in a larger variance compared with the case with no gravity. This can be understood by considering that gravity acts as an extra compressive forcing. This is equivalent to increasing $(bM)$ in compressible turbulent simulations (Equation (23)).

Therefore, according to the present analysis, we expect to have typically two regions (in terms of density) with different contributions governing the statistics of turbulence:

1. a first region, corresponding to $s < s_G$, where $s_G$ is given by Equation (25), where the statistics is similar to that of gravitationless turbulence but with a (more or less) increased variance due to gravity. The s-PDF in this region is log-normal-like,

2. a second region, corresponding to densities $s > s_G$, where gravity has a dominant impact on the statistics of turbulence, and the PDF will depart from (Gaussian) log-normal statistics.
The threshold density, $s_G$, between the two regions evolves with time on the same timescale $\bar{\tau}$ as the global, average properties of the cloud ($\bar{\rho}(t)$, $\alpha_{vir}(t)$, $(\nabla \cdot \mathbf{u})^2(t)$),<fnref n="3" /> according to Equation (25). However, as will be shown in Section 3, at densities $s > s_G$, the PDF will start departing from a log-normal form and develop a power law on shorter timescales, of the order of a typical local freefall time, $\tau_{ff}(s) < \bar{\tau}$.

### 3. Evolution of the Density PDF in Star-forming Clouds

Observations of column density PDFs in MCs show that regions where star formation has not yet occurred exhibit log-normal PDFs whereas regions with numerous prestellar cores exhibit power-law tails at high column densities (Kainulainen et al. 2009; Schneider et al. 2013). Similarly, numerical simulations of star formation in turbulent clouds show that density PDFs develop power-law tails as the simulations evolve (Klessen 2000; Federrath & Klessen 2013). This suggests that the density PDF in star-forming clouds is not stationary but evolves with time, implying $\langle \nabla \cdot \mathbf{u} | s \rangle \neq 0$ (see Section 2.3).

#### 3.1. Mathematical Derivation: Equivalence of the Velocity Divergence and s-PDF Power-law Tail Exponents

Finding solutions of Equation (9) for any function $\langle \nabla \cdot \mathbf{u} | s, t \rangle$ is not straightforward. Assuming, as a simplification, separability of the time and density variables $\langle \nabla \cdot \mathbf{u} | s, t \rangle = h(t) \times g(s)$ yields, from the method of characteristics, the solution

$$f(s, t) = \Phi \left( \int hdt' + \int \frac{1}{g} ds' \right) \frac{e^{-s}}{g(s)}, \quad (26)$$

with $\Phi$ any differentiable function.

We prove now that a nonstationary s-PDF develops a power-law tail of exponent $\alpha_s = a + 1$, with $a > 0$, if and only if the conditional expectation of the velocity divergence $\langle \nabla \cdot \mathbf{u} | s \rangle$ scales at large $s \gg 1$ as $\langle \nabla \cdot \mathbf{u} | s, t \rangle \approx h(t) \times e^{as}$. Indeed, if for $s \ge s_c$, for some $s_c$, one has

$$\langle \nabla \cdot \mathbf{u} | s, t \rangle = h(t) \times e^{as}, \quad (27)$$

with $a > 0$, then

$$f(s, t) = \Phi \left( a \int hdt' - e^{-as} \right) e^{-(1+a)s}, \quad (28)$$

for $s > s_c$, from Equation (27). Hence, as $a > 0$, the PDF is expected to develop a power-law tail with an exponent $-(1 + a)$ at a given time $t$ at $s > s_c$ sufficiently large that $\Phi(a \int hdt' - e^{-as}) \approx \Phi(a \int hdt')$. The proof of the reciprocal of this result (Equations (27) and (28)) is given in Appendix B. Thus, observed power-law tails, $f(s, t) \propto e^{-\alpha_s s}$, with exponents $\alpha_s = 3/2$ and 2 correspond to underlying expectations $\langle \nabla \cdot \mathbf{u} | s, t \rangle \propto e^{s/2}$ and $\langle \nabla \cdot \mathbf{u} | s, t \rangle \propto e^s$, respectively.

#### 3.2. Physical Interpretation: Transitory Regime and Short Time Evolution

At any time in a cloud, we can compute the threshold value $s_G$ above which gravity starts altering the statistics of fully developed turbulence significantly. For diffuse, hot, and/or

turbulent clouds ($\alpha_{vir} \gg 1$), however, this value can be so large that the probability $P(s > s_G)$ of finding regions $s > s_G$ becomes very small. In such cases, one can completely neglect the effect of gravity. To be more quantitative, let us assume that gravity can be neglected if $P(s > s_G) \le 10^{-9}$. Assuming a log-normal PDF, this yields $s_G \ge 6\sigma_s - 0.5\sigma_s^2$ (where $\sigma_s$ is the variance in Equation (23)). In hot and turbulent clouds, where $T \sim 8000$ K, ($bM \sim 1$ (e.g., Draco, Miville-Deschênes et al. 2017), this yields $\alpha_{vir} \gtrsim 5.5$ from Equation (25). As the cloud cools down and contracts, $\alpha_{vir}$ decreases, resulting in a value of $s_G$ small enough to observe significant departures from a log-normal PDF. We can thus define a time $t_0$ in the lifetime of the cloud as the time at which the volume fraction of the cloud corresponding to (dense) regions with $s > s_G$, where the gas PDF starts departing from the statistics of pure turbulence, $\langle \nabla \cdot \mathbf{u} | s \rangle \approx 0$, under the growing influence of gravity, becomes noticeable, i.e., statistically significant. This fixes the “zero of time” in the lifetimes of star-forming clouds, whatever the (undefinable) time since which they have been formed. The time $t_0$ thus corresponds to the time at which some dense regions start to collapse and depart from the global evolution (contraction or expansion) of the cloud, which is described by the time variation of $\bar{\rho}(t)$. This time $t_0$ then enables us to determine a physically motivated value to fix the indefinite integral in Equation (28), as the one being equal to 0 at $t_0$.

For regions with $s > s_G$, we expect from Equations (12) and (13) at short times $t = t_0 + t_{coll}$ after $t_0$, i.e., in the linear regime, to have $\langle \nabla \cdot \mathbf{u} | s \rangle \approx -4\pi G \bar{\rho} e^s (t - t_0) = -4\pi G \bar{\rho} e^{s_{coll}}$ (i.e., $a = 1$ in Equation (28)), yielding for the PDF $f$

$$f(s, t) \approx \tilde{\Phi} \left( \tau_{G,0}^{-2} \frac{t_{coll}^2}{2} + e^{-s} \right) e^{-2s}. \quad (29)$$

Therefore, for densities $s > s_G$, we expect to see the onset of a first power-law tail in the s-PDF, $f(s, t) \propto e^{-\alpha_s s}$, with a steep exponent $\alpha_s \approx 2$ in a typical timescale $\tau_G(s) = \tau_{G,0} e^{-s/2}$. As seen from Equation (29), the onset of this first power-law tail occurs, for a given time $t$, at a density $s_t = \rho_t(t_{coll})/\bar{\rho} \approx (\tau_{G,0}/t_{coll})^2 \approx 0.25(\tau_{ff,0}/t_{coll})^2$, as found in numerical calculations (Girichidis et al. 2014). At later time (a few $\tau_{ff}(s)$, see Section 3.3) for a given density, or at higher densities for a given time, a second power law develops with $\alpha_s = 3/2$, a signature of regions in freefall collapse, as seen in Section 3.3.

#### 3.3. Asymptotic Case: Evolution in Regions of “Freefall” Collapse

The densest regions in star-forming clouds are expected to collapse under their own gravity on a timescale of the order of the local freefall time $\tau_{ff}(\rho) \propto (G\rho)^{-1/2}$. For these regions we thus expect a scaling

$$-\langle \tau_{ff}^{-1}(\rho) | s \rangle \propto \langle \nabla \cdot \mathbf{u} | s \rangle = -c \sqrt{4\pi G \bar{\rho}} e^{s/2}, \quad (30)$$

where $c$ is a constant of proportionality of order unity. This yields, from Equation (28),

$$f(s, t) = \Phi \left( \frac{c}{2} \sqrt{4\pi G} \int_{t_0}^{t} \sqrt{\bar{\rho}(t')} dt' + e^{-s/2} \right) e^{-\frac{3}{2}s}, \quad (31)$$

where $t = t_0 + t_{coll}$. Then, if the time after which a dense region of the cloud started to collapse, $t_{coll}$, is short compared to the characteristic time of variation of $\bar{\rho}$, $t_{coll} \ll \bar{\tau}$, meaning that

<footnote n="3">³ The timescale $\bar{\tau}$ of variation of $\bar{\rho}$ is not necessary equal to $\tau_{ff,0}$. If there is enough turbulent support for example, it can be larger. It depends on what drives the global evolution of the cloud.</footnote>
the global properties of the cloud did not have time to evolve significantly, we can write

$$f(s, t) \approx \Phi\left(\frac{c}{2}\sqrt{4\pi G \bar{\rho}(t)} t_{\text{coll}} + e^{-s/2}\right)e^{-\frac{3}{2}s}. \quad (32)$$

Therefore, the PDF develops a power-law tail with a specific exponent $-3/2$ for $s \ge s_G$ within a typical time $t(s) \equiv 2c^{-1}\tau_{G,0} e^{-s/2} \simeq c^{-1}\tau_{ff,0} e^{-s/2} \simeq c^{-1}\tau_{ff}(s)$.

This analysis thus shows that the onset of power-law tails, $f(s) \propto e^{-\alpha_s s}$, in the PDF reflects the growing impact of gravity on the turbulent flow, with a first power-law exponent $\alpha_s \lesssim 2$, which reaches the asymptotic value $\alpha_s = 3/2$ in freefall collapsing regions.

### 4. From Volume to Column Densities

Observations of dense MCs trace the density integrated along the line of sight, and thus reveal the PDF of the column density $\Sigma$ or its logarithmic deviations $\eta = \ln(\Sigma/\bar{\Sigma})$. Many efforts have been made to link the observed $\eta$-PDF to properties of the underlying s-PDF (Vazquez-Semadeni & Garcia 2001; Brunt et al. 2010; Burkhart & Lazarian 2012; Federrath & Klessen 2013). In the present study, we will use the relation of Burkhart & Lazarian (2012) to illustrate our findings. Furthermore, we will adopt the relation given by Federrath & Klessen (2013) to link the exponents $\alpha_\eta$ and $\alpha_s$ of the $\eta$-PDF and s-PDF, respectively:

$$\alpha_{\eta} = -\frac{2}{1 - \frac{3}{\alpha_s}}. \quad (33)$$

Hence, for regions in freefall collapse we expect a power-law tail in the $\eta$-PDF with an asymptotic exponent $\alpha_\eta = 2$ ($\alpha_s = 3/2$), with a transition domain with $\alpha_\eta \ge 4$ ($\alpha_s \ge 2$) (see Section 3).

In order to make comparison between our theory and numerical simulations or observations for the $\eta$-PDF, we have derived a way to relate the volume density at which the s-PDF, $f(s)$, develops power laws to the column density at which the $\eta$-PDF, $p(\eta)$, develops a similar behavior. We use $s_{crit}$ and $\eta_{crit}$ to denote the critical value corresponding to the beginning of a power-law tail in the two respective PDFs. Assuming ergodicity and statistical isotropy, we obtain (see Appendix C for details)

$$\int_{\eta_{\text{crit}}}^{\infty} p(\eta) d\eta \approx \left( \int_{s_{\text{crit}}}^{\infty} f(s) ds \right)^{2/3}. \quad (34)$$

In the case where there are two power-law tails, starting at $s_1$ and $s_2$, Equation (34) remains a good approximation as long as $s_2 - s_1 \gtrsim 1$, so that the upper bound in the integrals is not very important. This procedure is tested against numerical simulations in Section 5 and confronted with observations in Section 6. More details on how one obtains Equation (34) are given in Appendix C.

### 5. Comparison with Numerical Simulations

#### 5.1. Numerical Setup

To understand how gravity affects the s- and $\eta$-PDFs in star-forming clouds and compare with our theoretical formulation, we use the numerical simulations of isothermal self-gravitating turbulence on 3D periodic grids presented in Federrath & Klessen (2012, 2013), kindly provided by the authors. These

simulations model isothermal self-gravitating magneto-hydrodynamic turbulence on 3D periodic grids with resolution $N_{res}^3 = 128^3-1024^3$. Here, we will only consider simulations with no magnetic field. In the simulations, turbulence is driven with solenoidal or compressive forcing or with a mixture of the two. Sink particles are used (see Federrath & Klessen 2012 or our Appendix D for details).

After Equation (25), we expect gravity to make a dominant contribution at densities $s > s_G$, which yields here

$$|e^{s_G} - 1| \equiv (bM)^2 \times \alpha_{\text{vir},0} \times \left| \frac{s_G + \frac{1}{2}\sigma_s^2}{\sigma_s^2} \right|, \quad (35)$$

where $\alpha_{vir,0} = 5\sigma_v^2/(6GL_b^2\rho_0)$ is the virial parameter suited for a box of size $L_b$ and 3D velocity dispersion $\sigma_v$, as in Federrath & Klessen (2012, 2013).$^4$ On the other hand, the maximum density $\rho_{max}$ above which the simulations do not properly resolve the collapse and describe the statistics of the cloud reads (Truelove et al. 1997)

$$\rho_{\text{max}} = \rho_0 e^{s_{\text{max}}} = \frac{\pi c_s^2}{16G \Delta x_{\text{min}}^2}, \quad (36)$$

with $\Delta x_{min}$ the size of the most resolved cell. This condition can be rewritten as

$$s_{\text{max}} = \ln(\alpha_{\text{vir},0}) + 2 \ln\left(\frac{N_{\text{res}}}{M}\right) + \ln\left(\frac{6\pi}{80}\right). \quad (37)$$

For $s > s_{max}$, the lack of resolution will yield the development of shallow power laws corresponding to the spurious fragmentation of these regions (see Figure 1 below and Federrath & Klessen 2013).

#### 5.2. Evolution of the PDFs

Figure 1 compares our analytical calculations of the s- and $\eta$-PDFs with the solenoidal simulations for $M \simeq 3$, $N_{res} = 512$ and $M = 5$, $N_{res} = 256$ at initial time $t = 0$ (orange), and at star formation efficiencies SFE = 0% (red) and SFE = 20% (dark red). A major result of Section 2.4 (Equations (25) and (35)) is the determination of a density threshold, $s_G$ (resp. $\eta_G$), above which the s-PDF (resp. $\eta$-PDF) is expected to develop power-law tails. Similarly, spurious shallow power laws will develop above $s_{max}$ (resp. $\eta_{max}$). In both cases $\eta_G$ and $\eta_{max}$ can be obtained from the determination of the corresponding values on the s-PDF with Equation (34). We note the excellent agreement between the theoretical and numerical curves over the whole range of densities. Notably, the theoretical determinations of $s_G$ from Equation (35), the threshold of the gravity-impacted domain, agree very well with the onset of a power law in the simulations. It is worth stressing that for the $M \sim 3$ simulation, $s_G \sim 0.1$ and $\alpha_{vir} \ll 1$, and thus we do not expect the power-law tail with exponent $\alpha_s = 3/2$ to develop up to $s_G$ in a time $\tilde{t} \simeq 1$ (since this requires typically a few $\tau_{ff,0}$). However, departures from a log-normal behavior are indeed seen to start at about $s \sim 0.1$.

Figure 2 compares $\eta_{mes} \equiv \eta_{G_{mes}}$, $\eta_{max_{mes}}$, directly measured on various simulations, to $\eta_{th} \equiv \eta_{G_{th}}$, $\eta_{max_{th}}$ derived from Equation (34) and the values of $s_G$ and $s_{max}$. As seen in the figure, the agreement between the theoretical value $\eta_{th}$ and the measured one $\eta_{mes}$ is remarkable.

<footnote>$^4$ Note that this differs from $\alpha_{vir,0}$ defined in Section 2.4 by a factor $\pi/6 \simeq 1/2$, if the cloud size $L_c$ is taken to be the box size $L_b$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Evolution of the s-PDF (top row) and η-PDF (bottom row) for solenoidal simulations with $M \simeq 3$, $N_{\text{res}} = 512$ (left) and $M = 5$, $N_{\text{res}} = 256$ (right) at $t = 0$ (orange circles), SFE = 0% (red stars), and SFE = 20% (dark red squares). Horizontal error bars represent bin spacing and vertical error bars indicate the uncertainty in the η-PDFs corresponding to three different projection directions of the simulation box. Log-normal fits of the low-density parts of PDFs at $t = 0$ and SFE = 20% are shown by dotted lines. The vertical dotted red lines correspond to the value of $s_G$ calculated from Equation (35) with values of ($bM$) and $\sigma_s$ calculated at time $\tilde{t} = t/\tau_{\text{ff,0}}$. For $s > s_G$, the s-PDFs and η-PDFs first develop power-law tails with exponents $\alpha_s = 2$, $\alpha_{\eta} = 4$ (black dashed lines) and then $\alpha_s = 3/2$, $\alpha_{\eta} = 2$ (green dotted-dashed lines) at higher density. The vertical dotted gray lines at $s > 6$ correspond to $s_{\text{max}}$ from Equation (37).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** $\eta_{th}$, calculated from Equation (34) and the s-PDF, vs. $\eta_{mes}$, directly measured on the η-PDFs, for runs with $M = 3, 5, 10$ (light to dark red). Runs with $M \sim 3$ and $M = 10$ have two values of $\eta_{max}$, one for each resolution $N_{\text{res}} = 256$ and $N_{\text{res}} = 512$. We note the excellent agreement between $\eta_{mes}$ and $\eta_{th}$.</figcaption></figure>

## 6. Comparison with Observations

In this section, we confront our theory with observations of column density PDFs in various MCs. We use a simple model with one or two power-law tails, characterized by one or two transition densities, $s_1$ and $s_2$, between log-normal and power-

laws, as described in Appendix E. From the determination of the variance $\sigma_{s,\eta}$ in the log-normal parts of the PDF, we get an estimate of the product ($bM$) (Equation (23)), while from the determinations of $s_1$ and $s_2$ we get an estimate of $\alpha_{\text{vir}} = 5\sigma_v^2/(\pi GL_c^2 \bar{\rho})$, and of the time since the first regions started to collapse, in units of mean freefall time $\tilde{t}_{\text{coll}} = t_{\text{coll}}/\tau_{\text{ff,0}}$ (Equations (25), (32), (29)). The values are given in Table 1.

We apply our method to two different clouds: Orion B (Schneider et al. 2013; Orkisz et al. 2017) and Polaris (André et al. 2010; Miville-Deschênes et al. 2010; Schneider et al. 2013). The data were kindly provided by Nicola Schneider. For Orion B, the average column density is $\bar{N}(H_2) = 2.06 \times 10^{21} \text{ cm}^{-2}$ and the cloud’s total mass and area above an extinction $A_v \ge 1$ are $M_{c,A_v \ge 1} = 29.69 \times 10^3 M_{\odot}$ and $A_{c,A_v \ge 1} = 651 \text{ pc}^2$. For Polaris this yields ($\bar{N}(H_2)$, $M_{c,A_v \ge 1}$, $A_{c,A_v \ge 1}$) = (1.73 $\times 10^{21} \text{ cm}^{-2}$, 1.21 $\times 10^3 M_{\odot}$, 3.9 pc$^2$).

The first one, Orion B, contains numerous prestellar cores. Its η-PDF displays a log-normal part at low column densities and a power-law tail at high densities with exponent $\alpha_{\eta} \simeq 2$, corresponding to an underlying s-PDF with exponent $\alpha_s = 3/2$, a signature of collapsed regions, as seen in Figure 3 (left). The power-law tail develops for $s > s_1 = 1.73_{-0.23}^{+0.25}$ pc$^2$. We can thus estimate that in this cloud (statistically significant, see Section 3) collapse of the densest regions has occurred since $\tilde{t}_{\text{coll}} \approx (2-5) \times e^{-s_1/2} \gtrsim 1$ (see Section 3.3). Note here that $\tau_{\text{ff,0}}$ corresponds to the region under study in the cloud, not to the global cloud itself. Estimation of ($bM$) from the determination of $\sigma_s$ combined with estimation of $s_1$ yield $\alpha_{\text{vir}} \sim 1$. The estimated
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Left: observed η-PDF of the cloud Orion B. Dashed-dotted black line: log-normal fit of the low-η part of the PDF. Dashed red and green lines: power law with exponent α<sub>η</sub> = 2, corresponding to an underlying s-PDF with a power-law exponent α<sub>s</sub> = 3/2, a signature of collapsed regions. Right: observed η-PDF of the cloud Polaris. Dashed-dotted black line: log-normal fit of the low-η part of the PDF. Dashed red and green lines: power laws respectively with exponents α<sub>η</sub> = 4 and 2, corresponding to an underlying s-PDF transiting from log-normal to power laws with exponents α<sub>s</sub> = 2 and 3/2, respectively.</figcaption></figure>

<table><caption>**Table 1**

Properties of the Clouds</caption><thead><tr><th>Name<br>(1)</th><th>Func. form<br>(2)</th><th>σ<sub>s</sub><br>(3)</th><th>(bM)<br>(4)</th><th>μ<br>(5)</th><th>α<sub>1</sub><br>(6)</th><th>α<sub>2</sub><br>(7)</th><th>s<sub>1</sub><br>(8)</th><th>s<sub>2</sub><br>(9)</th><th>α<sub>vir</sub><br>(10)</th><th>t&#770;<sub>coll</sub><br>(11)</th></tr></thead><tbody><tr><td>Orion B</td><td>Ln+1Pl</td><td>1.2<sup>+0.09</sup><sub>-0.1</sub></td><td>1.8<sup>+0.26</sup><sub>-0.26</sub></td><td>-0.92<sup>+0.13</sup><sub>-0.11</sub></td><td>2</td><td>...</td><td>1.73<sup>+0.25</sup><sub>-0.23</sub></td><td>...</td><td>0.88<sup>+0.26</sup><sub>-0.23</sub></td><td>≥1</td></tr><tr><td>Polaris</td><td>Ln+2Pl</td><td>0.98<sup>+0.07</sup><sub>-0.08</sub></td><td>1.27<sup>+0.16</sup><sub>-0.15</sub></td><td>-0.55<sup>+0.08</sup><sub>-0.07</sub></td><td>4</td><td>2</td><td>1.68<sup>+0.38</sup><sub>-0.34</sub></td><td>6.3<sup>+0.1</sup><sub>-0.15</sub></td><td>≤1.2</td><td>0.2 ± 0.1</td></tr></tbody></table>

**Note.** Columns: (1) cloud's name; (2) functional form: Ln+1Pl, Ln+2Pl (log-normal and 1 or 2 power laws); (3) standard deviation of the log-normal part σ<sub>s</sub>; (4) (bM) associated with σ<sub>s</sub>; (5) most probable s-value μ; (6) exponent of the first power law α<sub>1</sub>; (7) exponent of the second power law α<sub>2</sub>; (8) transition between the log-normal part and the first power law s<sub>1</sub>; (9) transition between the two power laws s<sub>2</sub>; (10) virial parameter α<sub>vir</sub> associated with s<sub>1</sub>; (11) time since the first region started to collapse in units of mean freefall time, t&#770;<sub>coll</sub>.

(bM) is compatible with mean Mach numbers M = 5.4<sup>+0.8</sup><sub>-0.8</sub> for b = 1/3 and M = 3.6<sup>+0.5</sup><sub>-0.5</sub> for b = 1/2, in agreement with Orkisz et al. (2017).

The second cloud, Polaris, where detectable star formation does not seem to have occurred yet, exhibits an extended power-law tail with a steep exponent, α<sub>η</sub> ≃ 4, corresponding to an s-PDF power-law tail of exponent α<sub>s</sub> = 2 for s > s<sub>1</sub>, before reaching the asymptotic values α<sub>η</sub> ≃ 2, i.e., α<sub>s</sub> = 3/2 at high density, s > s<sub>2</sub>, as seen in Figure 3 (right). Carrying out the same analysis as for Orion B, we get (s<sub>1</sub>, s<sub>2</sub>) = (1.68<sup>+0.38</sup><sub>-0.34</sub>, 6.3<sup>+0.1</sup><sub>-0.15</sub>). The value of s<sub>1</sub> yields here t&#770;<sub>coll</sub> ≈ 0.5 e<sup>-s<sub>1</sub>/2</sup> = 0.22<sup>+0.03</sup><sub>-0.04</sub> for this cloud. The determination of the density s<sub>2</sub>, which corresponds to collapsing regions, yields t&#770;<sub>coll</sub> ≈ (2–5) × e<sup>-s<sub>2</sub>/2</sup> ≈ 0.09–0.21, consistent with the above estimate of t&#770;<sub>coll</sub>, which we finally estimate as t&#770;<sub>coll</sub> = 0.2 ± 0.1. The theory thus suggests that gravity has started dominating dense regions, corresponding to the onset of the first power law at s = s<sub>1</sub>, only recently, i.e., for a short time t&#770;<sub>coll</sub>. According to these determinations, the quiescent Polaris region is quite young and has not even reached half its mean freefall time yet. Eventually, we expect it to start forming detectable prestellar cores on a timescale of the order of its mean freefall time, most likely in the “Saxophone” region, which entails most of the power-law part of its PDF (Schneider et al. 2013). Taking s<sub>1</sub> = s<sub>G</sub> yields an upper limit α<sub>vir</sub> ≲ 1.2. The estimated (bM) for Polaris yields mean Mach numbers M = 3.8<sup>+0.4</sup><sub>-0.4</sub> and M = 2.5<sup>+0.3</sup><sub>-0.3</sub> for

b = 1/3 and b = 1/2, respectively, consistent with the estimation of Schneider et al. (2013).

## 7. Conclusion

In this Letter, we have derived an analytical theory of the PDF of density fluctuations in supersonic turbulence in the presence of a gravity field in star-forming molecular clouds. The theory is based on a derivation of a combination of the coupled Navier–Stokes equations for the fluid motions and the Poisson equation for the gravity. The theory extends previous approaches (Pope 1981, 1985; Pan et al. 2018, 2019a, 2019b), first by including gravity and second by considering the PDF as a dynamical system, not a stationary one. We derive rigorously the transport equations of the PDF, characterize its evolution, and determine the density threshold above which gravity strongly affects and eventually dominates the dynamics of the turbulence. The theoretical results and diagnostics reproduce very well numerical simulations of gravoturbulent collapsing clouds (Section 5) and various available observations (Section 6).

A major result of the theory is the characterization of two density regions in the PDF (see Section 2.4). A low density region where gravity does not affect the dynamics of turbulence significantly, so the PDF is that of pure gravitationless turbulence, which resembles a log-normal form for isothermal, dominantly solenoidal turbulence. Then, above a density
threshold, $s_G$, given by Equations (25) and (35), gravity starts affecting the turbulence significantly, essentially by increasing the velocity dispersion (thus the variance). Above this threshold, $s > s_G$, power-law tails develop over time in the s-PDF, $f(s, t) \propto e^{-\alpha_s s}$, i.e., $p(\eta, t) \propto e^{-\alpha_\eta \eta}$ for the $\eta$-PDF of the surface density, as a direct consequence of the rising impact of gravity upon turbulence (see Section 3). Within a typical timescale $\sim \tau_G(s) = \tau_{G,0} e^{-s/2}$, with $\tau_{G,0} = 1/\sqrt{4\pi G\bar{\rho}}$, this yields the onset of a first power-law tail with $\alpha_s \ge 2$, i.e., $\alpha_\eta \ge 4$. Later on, after a few $\tau_{ff}(s)$ for a given density $s$, and/or at higher density, i.e., smaller scales, a second power law develops, with $\alpha_s = 3/2$, i.e., $\alpha_\eta = 2$. This is the signature of regions in freefall collapse.

Another important result of this study is to provide a procedure to relate the observed thresholds in column density, corresponding to the onset of the two power-law tails in the $\eta$-PDF, to the corresponding ones in volume density in the s-PDF (see Section 4 and Appendix A). Combined with the results of Sections 2.4 and 3, this allows us to infer, from the observation of column densities, various physical parameters characterizing molecular clouds (or regions of them), notably the virial parameter $\alpha_{vir}$. Moreover, the theory offers the possibility to date the clouds in units of $\tilde{t}_{coll}$, i.e., the time since a statistically significant fraction of dense regions of the cloud started to collapse, normalized to the cloud's mean freefall time. This explains why clouds exhibiting $\eta$-PDF with steep power laws ($\alpha_\eta \ge 3$) or extended log-normal parts are quiescent, since they have a short “age” $\tilde{t}_{coll}$. This applies to Polaris (André et al. 2010; Miville-Deschênes et al. 2010; Schneider et al. 2013) (Section 6) but could also explain the quiescence of Chamaeleon III (De Oliveira et al. 2014).

The theory derived in this study allows the determination of the aforementioned volume density and column density thresholds, $s_G$, $\eta_G$, and the characteristic timescales $\tau_G(s)$, $\tilde{t}_{coll}$ (Equations (25), (32), (29), (34), (35)). This yields quantitative, predictive diagnostics, from either simulations or observations, to determine precisely the relative impact of gravity upon turbulence within star-forming clouds/clumps and their evolutionary status. The theory thus provides a precise scale and clock to numericists and observers exploring star formation in MCs. It provides a sound theoretical foundation and quantitative diagnostics to analyze observations or numerical simulations of star-forming regions and to characterize the evolution of the density PDF in various regions of MCs. This theoretical framework provides a new vision on how gravitational collapse initiates and evolves within turbulent dense star-forming regions.

The authors are grateful to Christoph Federrath for sending us the PDFs of his simulations and to Nicola Schneider for the observational column densities presented in this article. We thank the anonymous referee for his/her insightful remarks that helped to improve the manuscript. We are also thankful to Benoit Commerçon, Jérémy Fensch, Guillaume Laibe, and Quentin Vigneron for helpful conversations.

## Appendix A
### Derivation of the Expression of the Source Terms, Equations (14)–(16)

To obtain the source terms, Equations (14)–(16), that appear in Equation (12), we start by taking the divergence of

Equation (2):

$$ 
\begin{aligned}
\partial_t(\nabla \cdot \mathbf{v}) + \nabla \cdot (\{\mathbf{v} \cdot \nabla\}\mathbf{u}) &+ \nabla \cdot (\{\mathbf{u} \cdot \nabla\}\mathbf{V}) + \nabla \cdot (\{\mathbf{V}\nabla\}\mathbf{V}) \\
&= -4\pi G \rho - \nabla \cdot \left(\frac{1}{\rho}\nabla P\right).
\end{aligned}
\quad (A1) $$

We then take the average of Equation (A2) to obtain

$$ \partial_t(\nabla \cdot \mathbf{V}) + \nabla \cdot (\{\mathbf{V} \cdot \nabla\}\mathbf{V}) = -4\pi G \bar{\rho}, \quad (A2) $$

where $\partial_t(\nabla \cdot \bar{\mathbf{u}}) = \nabla \cdot (\{\nabla \cdot \nabla\}\bar{\mathbf{u}}) = \nabla \cdot (\{\bar{\mathbf{u}} \cdot \nabla\}\mathbf{V}) = \nabla \cdot \left(\frac{1}{\rho}\nabla P\right) = 0$, because the turbulent fields $\rho$ and $\mathbf{u}$ are statistically homogeneous and because of the barotropic equation of state $P = P(\rho)$. Then, by subtraction, we obtain

$$ 
\begin{aligned}
\partial_t(\nabla \cdot \mathbf{u}) + \nabla \cdot (\{\mathbf{v} \cdot \nabla\}\mathbf{u}) + \nabla \cdot (\{\mathbf{u} \cdot \nabla\}\mathbf{V}) \\
&= -4\pi G \bar{\rho}(e^s - 1) - \nabla \cdot \left(\frac{1}{\rho}\nabla P\right).
\end{aligned}
\quad (A3) $$

We then note that $\nabla \cdot (\{\mathbf{u} \cdot \nabla\}\mathbf{V}) = \mathbf{u} \cdot \nabla(\nabla \cdot \mathbf{V}) + (\partial_i u_j)(\partial_j V_i) = (\partial_i u_j)(\partial_j V_i) = \nabla V: \nabla \mathbf{u}$, because $\rho$ is homogeneous, and by expanding $\nabla \cdot (\{\nabla \cdot \nabla\}\mathbf{u})$ we finally obtain

$$ 
\begin{aligned}
\frac{D \nabla \cdot \mathbf{u}}{Dt} = & - \nabla \mathbf{u}: \nabla \mathbf{u} - 2 \nabla V: \nabla \mathbf{u} - 4\pi G \bar{\rho}(e^s - 1) \\
& - \nabla \cdot \left(\frac{1}{\rho} \nabla P\right),
\end{aligned}
\quad (A4) $$

giving the expression of the source terms, Equations (14)–(16).

## Appendix B
### Reciprocal of Equations (27) and (28)

In Section 3 we have shown that a conditional expectation $\langle \nabla \cdot \mathbf{u}|s, t \rangle = h(t) \times e^{as}$, with $a > 0$, would produce an s-PDF with a power-law tail with exponent $\alpha_s = a + 1$, i.e. $f(s) \propto e^{-(a+1)s}$ (Equations (27)and (28)). We show here the reciprocal.

Let us assume that the s-PDF, $f$, is nonstatic and has a power-law tail with exponent $\alpha_s = a + 1$ with $a > 0$. More precisely, let us assume that

$$ f(s, t) = A(s, t) e^{-(a+1)s}, \quad (B1) $$

with a function $A(s, t)$ such as $A(s, t) \approx B(t)$ for $s \ge s_c$, for some $s_c$, where $B(t)$ is a $C^1$ function of the time variable only, with a bounded derivative. Rewriting Equation (9) as

$$ \left\{ \frac{\partial}{\partial s} + \frac{\partial \ln A}{\partial s} - a \right\} \langle \nabla \cdot \mathbf{u} | s, t \rangle = \frac{\partial \ln A}{\partial t}, \quad (B2) $$

one obtains

$$ 
\begin{aligned}
\langle \nabla \cdot \mathbf{u} | s, t \rangle = C(t) \times A(s, t)^{-1} \times e^{a s} + A(s, t)^{-1} e^{as} \\
\times \int_{s_i}^{s} e^{-as'} \frac{\partial A}{\partial t}(s', t) \, ds',
\end{aligned}
\quad (B3) $$

with $C(t)$ a function of the time variable only and $s_i$ some fixed density. As $f$ is not stationary, $\langle \nabla \cdot \mathbf{u} | s, t \rangle$ is not zero everywhere but, at any time $t$, there exists $s_0(t)$ such that $\langle \nabla \cdot \mathbf{u} | s_0(t), t \rangle = 0$ to ensure $\nabla \cdot \bar{\mathbf{u}} = 0$. Then we can fix the
function $C(t)$ to write without any loss of generality

$$ \langle \nabla \cdot \mathbf{u} | s, t \rangle = A(s, t)^{-1} e^{as} \int_{s_0(t)}^{s} e^{-as'} \frac{\partial A}{\partial t}(s', t) ds'. \quad (B4) $$

Then, because $a > 0$ and $\partial_t B(t)$ is bounded, the integral on the rhs of Equation (B4) is bounded and converges rapidly toward $\int_{s_0(t)}^{+\infty} e^{-as'} \partial_t A(s', t) ds' = I(t)$. The asymptotic behavior of $\langle \nabla \cdot \mathbf{u} | s, t \rangle$ for large $s \ge s_c$ is thus

$$ \langle \nabla \cdot \mathbf{u} | s, t \rangle \approx I(t) \times B(t)^{-1} \times e^{as} = h(t) \times e^{as}. \quad (B5) $$

This shows that to a nonstationary s-PDF with a power-law tail of exponent $\alpha_s = a + 1$ there corresponds a conditional expectation $\langle \nabla \cdot \mathbf{u} | s, t \rangle \approx h(t) \times e^{as}$ for $s \gg 1$.

## Appendix C

### Transitions to Power-law Tails

In Section 4 we derived a way to relate the volume density at which the s-PDF, $f(s)$, develops power laws to the column density at which the $\eta$-PDF, $p(\eta)$, develops a similar behavior. We call $s_{crit}$ the critical value corresponding to the beginning of a power-law tail in the s-PDF (Equation (34)).

Assuming ergodicity, we relate the volume fraction of regions with densities exceeding $s_{crit}$ to the probability of finding a density exceeding $s_{crit}$:

$$ \frac{V(s \ge s_{\text{crit}})}{V(\text{cloud})} = \int_{s_{\text{crit}}}^{\infty} f(s) ds. \quad (C1) $$

We now want to evaluate the area of this volume projected onto the plane perpendicular to the line of sight $S(s \ge s_{crit})$. Assuming statistical isotropy, we get

$$ \frac{S(s \ge s_{\text{crit}})}{S(\text{cloud})} \approx \left( \frac{V(s \ge s_{\text{crit}})}{V(\text{cloud})} \right)^{2/3}. \quad (C2) $$

We then identify regions in the observed area of the cloud contributing to the power law in the $\eta$-PDF with regions included in the projected area $S(s \ge s_{crit})$. This yields for the critical surface density $\eta_{crit}$ at which the $\eta$-PDF transits to a power law

$$ \int_{\eta_{\text{crit}}}^{\infty} p(\eta) d\eta = \frac{S(s \ge s_{\text{crit}})}{S(\text{cloud})} \approx \left( \int_{s_{\text{crit}}}^{\infty} f(s) ds \right)^{2/3}, \quad (C3) $$

which is Equation (34).

## Appendix D

### Numerical Models

In each simulation, gravity is switched on and sink particles are allowed to form after a state of fully developed turbulence has been reached, which determines the initial conditions at $t = t_0 = 0$ in the simulation. The associated transport equations for these simulations are

$$ V = 0, \quad (D1) $$

$$ \bar{\rho} = \rho_0, \quad (D2) $$

$$ s = \ln(\rho/\rho_0), \quad (D3) $$

$$ \frac{D \nabla \cdot \mathbf{u}}{Dt} = - \nabla \mathbf{u}: \nabla \mathbf{u} - 4\pi G \rho_0 (e^s - 1) \Theta(t) \\ - c_s^2 \nabla^2 s + \nabla \cdot \mathbf{F}_{\text{stir}}, \quad (D4) $$

$$ \frac{\partial}{\partial t} [\langle (\nabla \cdot \mathbf{u})^n | s \rangle f] = \left\{ 1 + \frac{\partial}{\partial s} \right\} [\langle (\nabla \cdot \mathbf{u})^{n+1} | s \rangle f] \\ + f \left\langle \frac{D(\nabla \cdot \mathbf{u})^n}{Dt} | s \right\rangle, \quad (D5) $$

where $\rho_0$ is constant, $c_s = 0.2 \text{ km s}^{-1}$ is the sound speed, $\nabla \cdot \mathbf{F}_{\text{stir}}$ is the divergence of the turbulent forcing, which is 0 for a solenoidal driving, and $\Theta(t)$ is the Heaviside step function ensuring that gravity is plugged in at $t = 0$. In all models the Mach number $M$ slightly increases with time because of collapsing regions. For most models, this only amounts to a few per cent, except for the $M \simeq 3$ simulations, which start at $M \simeq 2$ and end up at $M \simeq 3-4$ because the virial parameter $\alpha_{vir,0} = 5\sigma_v^2/(6GL_b^2\rho_0)$ is very small (see Table 1 in Federrath & Klessen 2012 or 2013). We note that the aforementioned definition of $\alpha_{vir,0}$ taken from Federrath & Klessen (2012, 2013) differs from the one we have introduced in Section 2.4 by a factor $\pi/6 \simeq 1/2$, if the cloud size $L_c$ is taken to be the box size $L_b$. As there is no unique way of translating the dimension of a cubic box into that of a spherical cloud, and in order to simplify the comparison between the simulations and our calculations, we keep their notation and definition.

Finally, to be consistent with the authors we describe the time evolution of the simulations by means of the reduced time $\tilde{t} = t/\tau_{\text{ff},0}$, which is the time in units of mean freefall time $\tau_{\text{ff},0} \equiv \sqrt{\frac{3\pi}{32G\rho_0}}$, and by means of the star formation efficiency (SFE), which is set at 0% at the formation of the first sink particle. The authors only extracted the PDFs up to SFE = 20%, which we will thus refer to as the “long time” of the runs.

## Appendix E

### Model with One or Two Power-law Tails

In this appendix, we develop a simple model that allows us to infer the global s-PDF of molecular clouds from the observations of $\eta$-PDFs. We assume that the PDFs are simply continuous and have only one power law at high densities and a log-normal cutoff at low densities:

$$ \begin{aligned} f(s) &= A_1 e^{-\frac{(s-\mu)^2}{2\sigma_s^2}}, \quad s \le s_{\text{crit}} \\ &= A_2 e^{-\alpha_s(s-s_{\text{crit}})}, \quad s \ge s_{\text{crit}}. \end{aligned} \quad (E1) $$

Enforcing the continuity and normalization of $f$ as well as the necessary condition $\bar{e}^s = 1$ (from our definition of $s$ in Equation (6)), we obtain

$$ A_1 = A_2 e^{\frac{(s_{\text{crit}} - \mu)^2}{2\sigma_s^2}} \quad (E2) $$

$$ 1 = \frac{1}{2} A_1 \sqrt{2\pi\sigma_s^2} \left[ 1 + \operatorname{erf}\left( \frac{s_{\text{crit}} - \mu}{\sqrt{2}\sigma_s} \right) \right] + \frac{A_2}{\alpha_s} \quad (E3) $$
$$1 = A_1 \sqrt{\frac{\pi}{2} \sigma_s^2} e^{\mu + \frac{\sigma_s^2}{2}} \left[ 1 + \operatorname{erf} \left( \frac{s_{\text{crit}} - \mu - \sigma_s^2}{\sqrt{2} \sigma_s} \right) \right] + \frac{A_2 e^{s_{\text{crit}}}}{\alpha_s - 1} \quad (E4)$$

We now assume that the variance $\sigma_s$ in the log-normal part and the exponent $\alpha_s$ of the power-law tail are inferred from the observations of the $\eta$-PDF following Section 4. More precisely, to obtain the variance $\sigma_s$, we use the formula of Burkhart & Lazarian (2012), $\sigma_{\eta}^2 = A_{\eta s} \times \sigma_s^2$, where $A_{\eta s}$ may depend on the forcing parameter $b$. For simulations of compressible turbulence without gravity and with solenoidal driving ($b = 1/3$) they found from their best fit $A_{\eta s} \simeq 0.11$, while observations of molecular clouds yield a value $A_{\eta s} \simeq 0.12-0.16$ for a forcing parameter $b = 0.5$, corresponding to a mixture of solenoidal and compressive driving. From Section 4, $A_2/\alpha_s$ and thus $A_2$ are obtained from the observations. We are now left with a system of three equations for three unknown quantities, namely $s_{\text{crit}}$, $\mu$, and $A_1$. We note that, in this model, the parameter $\mu$, which determines the peak of the log-normal part, is shifted to lower densities to ensure $\bar{e}^3 = 1$. Injecting Equation (E2) into (E3) yields a closed equation for the variable $x = \frac{s_{\text{crit}} - \mu}{\sqrt{2}\sigma_s}$:

$$1 = A_2 e^{x^2 \sqrt{2\pi\sigma_s^2}} \Phi(x) + \frac{A_2}{\alpha_s}, \quad (E5)$$

with $\Phi(x) = \frac{1}{2}[1 + \operatorname{erf}(x)]$ the cumulative distribution function for the normal distribution. Equation (E4) is then used to obtain $\mu$ and then $s_{\text{crit}}$.

In case where the $\eta$-PDF exhibits two power-law tails with exponents $\alpha_{\eta} = 4$ and $\alpha_{\eta} = 2$ we simply assume the following functional form:

$$\begin{aligned} f(s) &= A_1 e^{-\frac{(s-\mu)^2}{2\sigma_s^2}}, \quad s \le s_1 \\ &= A_2 e^{-\alpha_1(s-s_1)}, \quad s_1 \le s \le s_2 \\ &= A_2 e^{-\alpha_1(s_2-s_1)} e^{-\alpha_2(s-s_2)}, \quad s_2 \le s, \end{aligned} \quad (E6)$$

with $\alpha_1 = 2$ and $\alpha_2 = 3/2$, and change the procedure as follows. First, we build the $s$-PDF as if there were only one power-law with exponent $\alpha_{\eta} = 4$ in the $\eta$-PDF with the aforementioned procedure to obtain $A_1$, $A_2$, $\mu$, and $s_1$. We then use Equation (34) to obtain $s_2$:

$$\frac{A_2 e^{-\alpha_1(s_2-s_1)}}{\alpha_2} = \left( \int_{\eta_2}^{\infty} p(\eta) d\eta \right)^{3/2}, \quad (E7)$$

where $\eta_2$ is the column density at the beginning of the second power law with exponent $\alpha_{\eta} = 2$. This modified procedure, while simple to implement, is sufficiently accurate because $s_2$ is

large and thus regions with $s > s_2$ only represent a particularly small fraction of the total volume ($\lesssim 10^{-5}$).

We confront this procedure with observations in Section 6. Errors arising from the determination of $\eta_{\text{crit}}$ and $\sigma_s$ from the observations yield an error $\Delta s_{\text{crit}} = \pm 0.3$ on $s_{\text{crit}}$, which is reasonable.

## ORCID iDs

Etienne Jaupart https://orcid.org/0000-0001-6853-295X

Gilles Chabrier https://orcid.org/0000-0002-8342-9149

## References

* André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102

* Ballesteros-Paredes, J., Vázquez-Semadeni, E., Gazol, A., et al. 2011, MNRAS, 416, 1436

* Brunt, C. M., Federrath, C., & Price, D. J. 2010, MNRAS, 403, 1507

* Burkhart, B., & Lazarian, A. 2012, ApJL, 755, L19

* Burkhart, B., Stalpes, K., & Collins, D. C. 2016, ApJL, 834, L1

* Cho, W., & Kim, J. 2011, MNRAS: Letters, 410, L8

* Collins, D. C., Kritsuk, A. G., Padoan, P., et al. 2012, ApJ, 750, 13

* De Oliveira, C. A., Schneider, N., Merín, B., et al. 2014, A&A, 568, A98

* Donkov, S., & Stefanov, I. Z. 2018, MNRAS, 474, 5588

* Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

* Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51

* Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJL, 688, L79

* Federrath, C., Roman-Duval, J., Klessen, R., Schmidt, W., & Mac Low, M.-M. 2010, A&A, 512, A81

* Frisch, U. 1995, Turbulence: The Legacy of A.N. Kolmogorov (Cambridge: Cambridge Univ. Press)

* Girichidis, P., Konstandin, L., Whitworth, A. P., & Klessen, R. S. 2014, ApJ, 781, 91

* Guszejnov, D., Hopkins, P. F., & Grudić, M. Y. 2018, MNRAS, 477, 5139

* Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, A&A, 508, L35

* Kainulainen, J., Lehtinen, K., & Harju, J. 2006, A&A, 447, 597

* Klessen, R. S. 2000, ApJ, 535, 869

* Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416

* Kritsuk, A. G., Norman, M. L., & Wagner, R. 2010, ApJL, 727, L20

* Ledoux, P., & Walraven, T. 1958, Astrophysics II: Stellar Structure/ Astrophysik II: Sternaufbau (Berlin: Springer), 353

* Lee, E. J., Chang, P., & Murray, N. 2015, ApJ, 800, 49

* Lemaster, M. N., & Stone, J. M. 2008, ApJL, 682, L97

* Miville-Deschênes, M.-A., Martin, P., Abergel, A., et al. 2010, A&A, 518, L104

* Miville-Deschênes, M.-A., Salomé, Q., Martin, P., et al. 2017, A&A, 599, A109

* Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99

* Pan, L., Padoan, P., & Nordlund, Å. 2018, ApJL, 866, L17

* Pan, L., Padoan, P., & Nordlund, Å. 2019a, ApJ, 881, 155

* Pan, L., Padoan, P., & Nordlund, Å. 2019b, ApJ, 876, 90

* Passot, T., & Vázquez-Semadeni, E. 1998, PhRvE, 58, 4501

* Pope, S. 1981, PhFl, 24, 588

* Pope, S., & Ching, E. S. 1993, PhFlA, 5, 1529

* Pope, S. B. 1985, PrECS, 11, 119

* Schneider, N., André, P., Könyves, V., et al. 2013, ApJL, 766, L17

* Schneider, N., Csengeri, T., Hennemann, M., et al. 2012, A&A, 540, L11

* Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJL, 489, L179

* Vázquez-Semadeni, E. 1994, ApJ, 423, 681

* Vázquez-Semadeni, E., & Garcia, N. 2001, ApJ, 557, 727