#  The density structure of supersonic self-gravitating turbulence

Shivan Khullar,1,2 Christoph Federrath,3,4 Mark R. Krumholz, 3,4 Christopher
D. Matzner1
1David A. Dunlap Dept. of Astronomy & Astrophysics, University of Toronto, 50
St. George St., Toronto, ON M5S 3H4, Canada
2Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George St., Toronto, ON M5S 3H8, Canada
3Research School of Astronomy and Astrophysics, Australian National
University, Cotter Rd., Weston ACT 2611, Australia
4ARC Centre of Excellence for Astronomy in Three Dimensions (ASTRO-3D),
Canberra ACT 2601, Australia E-mail<EMAIL_ADDRESS>(SK)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We conduct numerical experiments to determine the density probability
distribution function (PDF) produced in supersonic, isothermal, self-
gravitating turbulence of the sort that is ubiquitous in star-forming
molecular clouds. Our experiments cover a wide range of turbulent Mach number
and virial parameter, allowing us for the first time to determine how the PDF
responds as these parameters vary, and we introduce a new diagnostic, the
dimensionless star formation efficiency versus density ($\epsilon_{\rm
ff}(s)$) curve, which provides a sensitive diagnostic of the PDF shape and
dynamics. We show that the PDF follows a universal functional form consisting
of a log-normal at low density with two distinct power law tails at higher
density; the first of these represents the onset of self-gravitation, and the
second reflects the onset of rotational support. Once the star formation
efficiency reaches a few percent, the PDF becomes statistically steady, with
no evidence for secular time-evolution at star formation efficiencies from
about five to 20 percent. We show that both the Mach number and the virial
parameter influence the characteristic densities at which the log-normal gives
way to the first power-law, and the first to the second, and we extend (for
the former) and develop (for the latter) simple theoretical models for the
relationship between these density thresholds and the global properties of the
turbulent medium.

###### keywords:

hydrodynamics – turbulence – stars:formation – dust, extinction – ISM: clouds

††pubyear: 2021††pagerange:  The density structure of supersonic self-
gravitating turbulence –B

##  1 Introduction

The interplay of turbulence and gravity lies at the heart of many
astrophysical problems. Over the last few decades, there has been tremendous
progress toward understanding astrophysical turbulence in a statistical sense
(e.g., Vazquez-Semadeni 1994; Padoan et al. 1997; Nordlund & Padoan 1999;
Elmegreen 2002b; Kritsuk et al. 2007). The density probability distribution
function (PDF) has shown particular promise as a tool to study turbulent gas
flows in the presence of gravity (e.g. Ballesteros-Paredes et al. 1998;
Elmegreen 2002a; Krumholz & McKee 2005; Dib et al. 2007; Federrath et al.
2008a; Federrath & Banerjee 2015; Nolan et al. 2015; Murray & Chang 2015;
Federrath et al. 2016a; Burkhart et al. 2016; Pan et al. 2018, 2019; Jaupart &
Chabrier 2020). The density PDF serves as an important input to analytical
theories of turbulence-regulated star formation, which link estimates of the
star formation rate or the initial mass function of stars directly to the PDF
(see e.g., Padoan & Nordlund 2002, 2004; Krumholz & McKee 2005; Padoan et al.
2007; Hennebelle & Chabrier 2008, 2009, 2011; Padoan & Nordlund 2011;
Federrath & Klessen 2012; Hopkins 2013a). It is also closely related to the
directly-observable column density distribution function, the N-PDF (Brunt et
al. 2010a, b; Ginsburg et al. 2013; Kainulainen et al. 2014). For these
reasons, it is critical to understand the physical processes shaping the
density PDF.

A number of studies have shown that, for supersonic isothermal turbulence, the
density PDF assumes a log-normal (LN) distribution, at least for low Mach
numbers (Vazquez-Semadeni, 1994; Padoan et al., 1997; Ostriker et al., 2001;
Kritsuk et al., 2007; Lemaster & Stone, 2008; Federrath et al., 2008a, 2010a;
Konstandin et al., 2012; Molina et al., 2012; Nolan et al., 2015). This would
be expected from an application of the central limit theorem: the build up of
density ($\rho$) fluctuations (such as shocks) should be random and
multiplicative in nature, or equivalently, for the variable
$s=\ln(\rho/\rho_{0})$, this is a random additive process which will result in
a Gaussian distribution (LN in $\rho$: Pope & Ching 1993; Passot & Vázquez-
Semadeni 1998). However, some numerical simulations have shown deviations from
the LN distribution as well (Federrath et al. 2008b; Schmidt et al. 2009;
Federrath et al. 2010a; Price & Federrath 2010; Konstandin et al. 2012;
Federrath & Klessen 2013; Pan et al. 2019). Hopkins (2013b) argues that these
appear because the density fluctuations are not uncorrelated and thus the
central limit theorem is not fully applicable. Despite these caveats, the LN
distribution remains a good approximation for turbulent supersonic media, and
is consistent with observations of molecular clouds (Berkhuijsen & Fletcher
2008; Hill et al. 2008; Burkhart et al. 2010; Burkhart & Lazarian 2012; Maier
et al. 2016; Federrath et al. 2016b; Maier et al. 2017; Chen et al. 2018;
however, for a contrary view, see Lombardi et al. 2015; Alves et al. 2017 who
suggest that the distribution follows a power law modified at low densities by
the finite size of a given map).

Real molecular clouds are self-gravitating in addition to being magnetized and
turbulent, and virtually all numerical studies that include gravity find that
it modifies the density PDF (e.g. Ballesteros-Paredes et al. 2011; Kritsuk et
al. 2011; Collins et al. 2012; Federrath & Klessen 2013; Girichidis et al.
2014; Burkhart et al. 2015). Self-gravity causes the turbulent (LN)
distribution to develop a power-law (PL) tail at high densities. An analogous
tail is observed in the N-PDF of several observed molecular clouds
(e.g.,Schneider et al. 2013, 2015c, 2015b; Pokhrel et al. 2016; Dib et al.
2020). This PL tail forms when gravitational forces become important, such as
in the process of collapse. Kritsuk et al. (2011) invoke models based on the
collapse of singular isothermal spheres (Penston 1969a, b; Larson 1969; Shu
1977; Hunter 1977; Whitworth & Summers 1985) to explain the exponent of the PL
tail. Collins et al. (2012) find that the inclusion of self-gravity leads to
the division of the density PDF into two distinct phases, a low-density
turbulent cloud and high-density, self-gravitating cores. Girichidis et al.
(2014) develop an analytical description for how this density PDF evolves
during the free-fall collapse of a uniform density sphere. These authors find
that the PL tail of the (volume-weighted) density PDF asymptotes to a slope
$-1.54$. They argue that physical processes such as driven converging flows
can accelerate collapse and flatten the slope of the PL tail. Jaupart &
Chabrier (2020) build on the analytical results of Pan et al. (2018, 2019) and
develop a theory for the shape and evolution of the density PDF from the
Navier-Stokes and Poisson equations. They determine a density threshold above
which gravity would dominate the dynamics of turbulence, which depends on the
cloud parameters such as the Mach number and virial parameter. In their
theory, PL tails develop naturally over time-scales of a few mean free-fall
times.

High-resolution simulations also suggest the existence of a second PL tail at
higher densities, which is distinguished from the first tail by having a
different slope (Kritsuk et al., 2011; Collins et al., 2011, 2012; Murray et
al., 2017; Veltchev et al., 2019; Marinkova et al., 2020). However, the origin
of this feature remains disputed. Kritsuk et al. (2011) and Murray et al.
(2017) suggest that it is due to pile up of mass caused by the formation of
centrifugally-supported discs, but Kritsuk et al. also speculate that it could
be a numerical artefact created when simulations violate the Truelove et al.
(1997) resolution condition at the highest densities. Jaupart & Chabrier
(2020) argue that a second PL tail, with a characteristic slope of 1.5, is a
signature of regions in free-fall collapse. Donkov et al. (2020) model clouds
as self-gravitating spheres, and treat turbulent support as a modified
equation of state (following Donkov & Stefanov 2018, 2019). In their model, a
second PL tail emerges as a result of the transition between predominantly
turbulent support at larger radii and predominantly thermal support at smaller
radii. There is also some observational evidence for second PL tails in
N-PDFs, mostly from analysis of Herschel-derived dust column density maps
(Tremblin et al., 2014; Schneider et al., 2015a; Pokhrel et al., 2016; Dib et
al., 2020). It is unclear, however, if these features are due to the same
processes that drive the second PL features in simulated volume density PDFs.

Current uncertainties about the functional form of the density PDF, and its
physical origin, provide the primary motivation for this paper. To test these,
we perform a suite of simulations of isothermal gravo-turbulent fluids in
which we systematically vary the dimensionless parameters – the Mach number,
$\mathcal{M}$, and the virial parameter, $\alpha_{\rm vir,0}$ – that describe
the gas. We characterise the variation of the PDF with these parameters and
over time within a single simulation, and present analytic models to interpret
these results. In the process, we make extensive use of the intimate links
between density PDFs and the star formation rate per free-fall time,
$\epsilon_{\rm ff}$, which is one of the central quantities that any
successful theory of star formation should be able to predict (e.g. Krumholz &
McKee 2005; Federrath & Klessen 2012, 2013; Salim et al. 2015; Vutisalchavakul
et al. 2016; Heyer et al. 2016; Leroy et al. 2017; Grisdale et al. 2019).
Khullar et al. (2019) show that the variation of this quantity with density,
$\epsilon_{\rm ff}(\rho)$, is a powerful diagnostic that can be used to
identify changes in behaviour at characteristic densities, often referred to
as star formation thresholds, because it characterizes the rate at which
gravity drives the evolution from low to high densities. We extend this
analysis to show that $\epsilon(\rho)$ can also be used to understand the
origin and variation of the density PDF, and that the two functions can be
regarded as two sides of the same coin.

The rest of the paper is organized as follows: in Section 2, we describe our
numerical simulations and our method for extracting PDF information from them.
In Section 3, we discuss the physical processes controlling the density PDF
and look at the variation of the PDF parameters with the simulation
parameters. In Section 4.2, we present models for the dependence of the PDF on
Mach number and virial parameter. We summarize our conclusions in Section 5.

##  2 Simulations and Methods

###  2.1 Simulations

Table 1: Summary of the simulation suite. Column 1: Simulation name. Columns 2-4: Virial Parameter as defined by equation 5, (3D) Velocity dispersion, Mach number $\mathcal{M}$. Columns 5-7: Box size, Effective grid size (resolution), Effective minimum cell size. Columns 8-10: Total mass, Mean density, Sink particle threshold density. Columns 11-12: Autocorrelation Time, Free-fall time at the mean density. Name | $\alpha_{\rm vir,o}$ | $\sigma_{V}$ | $\mathcal{M}$ | $L$ | $N_{\rm res}$ | $\Delta x$ | $M$ | $\rho_{0}$ | $\rho_{\rm sink}$ | $T_{\rm auto}$ | $t_{\rm ff}$
---|---|---|---|---|---|---|---|---|---|---|---
|  | [km/s] |  | [pc] |  | [AU] | [M⊙] | [g cm-3] | [g cm-3] | [Myr] | [Myr]
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12)
V0.5M2.5 | 0.5 | 0.5 | 2.5 | 2 | 20483 | 201.45 | 193.86 | 1.64 $\times$ 10-21 | 8.29 $\times$ 10-17 | 1.957 | 1.644
V1M2.5 | 1 | 0.5 | 2.5 | 2 | 20483 | 201.45 | 96.93 | 0.82 $\times$ 10-21 | 8.29 $\times$ 10-17 | 1.957 | 2.326
V0.5M5 | 0.5 | 1.0 | 5 | 2 | 20483 | 201.45 | 775.44 | 6.56 $\times$ 10-21 | 8.29 $\times$ 10-17 | 0.978 | 0.822
V1M5 | 1 | 1.0 | 5 | 2 | 20483 | 201.45 | 387.72 | 3.28 $\times$ 10-21 | 8.29 $\times$ 10-17 | 0.978 | 1.163
V2M5 | 2 | 1.0 | 5 | 2 | 20483 | 201.45 | 193.86 | 1.64 $\times$ 10-21 | 8.29 $\times$ 10-17 | 0.978 | 1.645
V0.5M10 | 0.5 | 2.0 | 10 | 2 | 20483 | 201.45 | 3101.75 | 2.62 $\times$ 10-20 | 8.29 $\times$ 10-17 | 0.489 | 0.411
V1M10 | 1 | 2.0 | 10 | 2 | 81923 | 50.35 | 1550.87 | 1.31 $\times$ 10-20 | 1.33 $\times$ 10-15 | 0.489 | 0.582
V2M10 | 2 | 2.0 | 10 | 2 | 20483 | 201.45 | 775.44 | 6.56 $\times$ 10-21 | 8.29 $\times$ 10-17 | 0.489 | 0.822
V2M5res256 | 2 | 1.0 | 5 | 2 | 2563 | 1611.61 | 193.86 | 1.64 $\times$ 10-21 | 1.30 $\times$ 10-18 | 0.978 | 1.645
V2M5res512 | 2 | 1.0 | 5 | 2 | 5123 | 805.80 | 193.86 | 1.64 $\times$ 10-21 | 5.18 $\times$ 10-18 | 0.978 | 1.645
V2M5res1024 | 2 | 1.0 | 5 | 2 | 10243 | 402.90 | 193.86 | 1.64 $\times$ 10-21 | 2.08 $\times$ 10-17 | 0.978 | 1.645
V1M10res2048 | 1 | 2.0 | 10 | 2 | 20483 | 201.45 | 1550.87 | 1.31 $\times$ 10-20 | 8.29 $\times$ 10-17 | 0.489 | 0.582
V2M5res4096 | 2 | 1.0 | 5 | 2 | 40963 | 100.73 | 193.86 | 1.64 $\times$ 10-21 | 3.32 $\times$ 10-16 | 0.978 | 1.645

We perform a suite of 13 simulations using the Adaptive Mesh Refinement (AMR,
Berger & Colella 1989) code FLASH (Fryxell et al., 2000; Dubey et al., 2008).
We use a setup similar to that of Federrath & Klessen (2012, hereafter, FK12),
but with higher resolution, and encompassing a range of dimensionless
parameters necessary to study the variation of the density PDF with them. We
solve the following equations of hydrodynamics and self-gravity in three
dimensions:

$\displaystyle\frac{\partial\rho}{\partial
t}+\mathbf{\nabla}\cdot(\rho\mathbf{v})=0,$ (1)
$\displaystyle\rho\frac{\partial\mathbf{v}}{\partial
t}+\rho(\mathbf{v}\cdot\mathbf{\nabla})\mathbf{v}=-\rho(\mathbf{g}+\mathbf{F_{\rm
turb}})-\mathbf{\nabla}(P_{\rm th}),$ (2)
$\displaystyle\mathbf{g}=-\mathbf{\nabla}\phi_{\rm gas}+\mathbf{g_{\rm
sink}},$ (3) $\displaystyle\nabla^{2}\phi_{\rm gas}=4\pi G\rho,$ (4)

where $\rho$ and $\mathbf{v}$ denote the gas density and velocity
respectively, and $\mathbf{g}_{\rm sink}$ is the gravity of the sink particles
(§ 2.1.2). Finally, to close the hydrodynamics equations we use an isothermal
equation of state $P_{\rm th}=c_{s}^{2}\rho$ with uniform sound speed
$c_{s}=({k_{B}T/\mu m_{\rm H}})^{1/2}$. While each simulation can be scaled to
different physical parameters, we nominally set $c_{s}=0.2$ km s-1,
corresponding to gas temperature $T\approx 11$ K and mean molecular weight
$\mu=2.3$. The isothermal approximation is reasonable for molecular cloud
environments (Krumholz, 2015) up to a characteristic critical density $\sim
10^{14}$ cm-3 (Masunaga & Inutsuka, 2000). Although our simulations do exceed
this density at times, we maintain the isothermal approximation in order to
focus on the hydrodynamics of the density distribution, rather than the
thermal physics that determines fragmentation.

Our simulation domain is a periodic box of size $L=2$ pc on each side. We
discretise the box with maximum effective resolutions of $N_{\rm
res}=512^{3}$–$8192^{3}$, with a base grid of 5123 with up to 5 levels of AMR
and a refinement criterion that kicks in when the local Jeans length,
$\lambda_{J}=(\pi c_{s}^{2}/G\rho)^{1/2}$, falls below 16 grid cells. This
refinement criterion ensures that the local Jeans length is well-resolved at
any point in space and time, and that the turbulence can be reasonably well
resolved on the Jeans scale. We use a positive-definite Riemann solver (Waagan
et al., 2011) for solving the system of hydrodynamic equations and a multigrid
Poisson solver (Ricker, 2008) for the self-gravity of the gas with periodic
boundary conditions.

We caution readers that our simulations do not include magnetic fields or
stellar feedback, both of which affect the star formation efficiency and thus
likely the density PDF as well (e.g., Krumholz et al., 2012; Myers et al.,
2014; Federrath, 2015; Schneider et al., 2015a; Lin et al., 2016; Mocz et al.,
2017; Grudić et al., 2018; Cunningham et al., 2018). However, including these
effects would dramatically expand the parameter space, and in particular
including feedback effects would likely make the results dependent on the
details of the subgrid feedback recipe adopted. For this reason, purely
hydrodynamic simulations, although less realistic, offer the advantage of a
clean, more easily interpreted experiment. They are also valuable in that they
form a baseline case against which magnetic and feedback effects can be
identified in future work.

####  2.1.1 Turbulence Driving and Initial Conditions

Our simulations run in two phases: we first drive the turbulence to
statistical steady state with gravity turned off, and then turn on gravity and
allow collapse. The procedure used to drive turbulence is similar to previous
simulations (FK12; Federrath 2015; Mathew & Federrath 2020), and uses the
basic method described in detail in Federrath et al. (2010a). We apply a
stochastic Ornstein-Uhlenbeck process (Eswaran & Pope, 1988) to construct an
acceleration field $\mathbf{F}_{\rm turb}$, which appears as a source term in
equation 2. $\mathbf{F}_{\rm turb}$ contains only large-scale modes, $1\leq
k\leq 3$ (where the wavenumber $k$ is in units of $2\pi/L$), and most of the
power is injected at the $k=2$ mode, corresponding to half the box size. Our
turbulence driving module allows us to control the mixture of modes
$\mathbf{F}_{\rm turb}$ excites, so we can drive with purely solenoidal modes
($\nabla\cdot\mathbf{F}_{\rm turb}=0$), purely compressive modes
($\nabla\times\mathbf{F}_{\rm turb}=0$), or a natural mixture of the two. We
use natural mixtures for all the simulations presented here (Federrath et al.,
2010a).

We initialize the simulation box with gas having uniform density $\rho_{0}$
and zero velocity, then drive turbulence for two turnover (or auto-
correlation) times, $T_{\rm auto}=L/2\sigma_{V}$, to allow the turbulent
cascade to build up. Here $\sigma_{V}$ is the steady-state velocity dispersion
produced by our driving. We switch on self-gravity when $t_{\rm tot}=2T_{\rm
auto}$ and denote this time as $t=0$. We continue to drive turbulence for all
times during the evolution of the simulated clouds.

Our simulations can be characterized by three key dimensionless parameters:
(a) the Mach number of the flow, $\mathcal{M}=\sigma_{V}/c_{s}$, (b) the
virial ratio, $\alpha_{\mathrm{vir}}=2E_{\rm kin}/|E_{\rm pot}|$, given as
(under the assumption of a uniform, spherical gas cloud of radius $L/2$ in
isolation),

$\alpha_{\rm vir,0}=\frac{5\sigma_{V}^{2}L}{6GM_{c}},$ (5)

and (c) the turbulence driving parameter, $b$ (Federrath et al. 2008b, 2010a).
The symbols $E_{\rm kin}$, $E_{\rm pot}$, $\sigma_{V}$ and $M_{c}$ denote the
kinetic energy, potential energy, three dimensional velocity dispersion, and
mass of the gas cloud at $t=0$. Our choice of a natural mix of compressive and
solenoidal modes corresponds to $b=0.4$ in all cases, so the simulations can
be described by $\alpha_{\rm vir,0}$ and $\mathcal{M}$ alone. We run
simulations with $\alpha_{\rm vir,0}=0.5$–$2$ and $\mathcal{M}\approx
2.5$–$10$. We summarize these and other key simulation parameters in Table 1.

####  2.1.2 Sink Particles

To model collapsing gas and stars in our simulations, we use sink particles
(Bate & Burkert 1997; Krumholz et al. 2004; Federrath et al. 2010b). When the
gas exceeds a specified density $\rho_{\rm sink}$ (at the maximum level of
refinement), we create a spherical control volume with radius $r_{\rm sink}$,
and perform additional checks on the collapsing region to determine if the gas
in the control volume, for example, is bound and Jeans unstable. (For details,
see Federrath et al. 2010b.) The sink particle threshold density is given by

$\rho_{\rm sink}=\frac{\pi c_{s}^{2}}{G\lambda_{\rm J}^{2}},$ (6)

where $\lambda_{\rm J}$ is the Jeans length at $\rho=\rho_{\rm sink}$; we give
numerical values of $\rho_{\rm sink}$ for our simulations in Table 1. The
control volume radius $r_{\rm sink}=\lambda_{\rm J}/2$ at $\rho=\rho_{\rm
sink}$ is set to 2.5 grid cell lengths at the maximum AMR level. This ensures
that the Jeans length is resolved by at least 5 grid cells on the highest
level of AMR, consistent with the Truelove et al. (1997) criterion to avoid
artificial fragmentation. For more details about the implementation of sink
particles, we refer the reader to Federrath et al. (2010b, 2014).

###  2.2 Density PDF fitting model

A primary goal of this work is to study the shape of the volumetric density
PDF, and its dependence on $\alpha_{\rm vir,0}$ and $\mathcal{M}$. For this
purpose, it is helpful to describe the PDF in terms of a few simple parameters
that we can extract by fitting to the distribution of cell densities found in
our simulations (either in a single time snapshot or averaged over a number of
snapshots). As we show below, and consistent with the conclusions of previous
works (Section 1), the PDFs in our simulations can be described reasonably
well by a piecewise function consisting of a lognormal (LN) component at low
density and one or two power-law (PL) components at high density. We first
introduce the single power-law form for the PDF since it has been used
traditionally by previous works (see e.g., Myers 2015; Burkhart et al. 2017;
Burkhart 2018; Burkhart & Mocz 2019). In terms of the natural logarithm of the
density, $s=\ln(\rho/\rho_{0})$, where $\rho_{0}$ is the mean density in the
computational volume, the single PL form, to which we shall henceforth refer
as LN+PL, is given by

$p(s)=\begin{cases}\frac{N}{\sqrt{2\pi\sigma_{\rm
s}^{2}}}\exp\left[-\frac{(s-s_{0})^{2}}{2\sigma_{\rm s}^{2}}\right]&s<s_{\rm
g},\\\ Np_{0}e^{-s\alpha_{\rm g}}&s\geq s_{\rm g},\end{cases}$ (7)

where $\alpha_{\rm g}$ is the slope of the power-law tail, $\sigma_{\rm s}$ is
the width of the LN part, $s_{0}$ is the mean of the LN part and $N$ and
$p_{0}$ are normalisation factors. This functional form contains 6 variables,
$N,p_{0},s_{0},\sigma_{\rm s},s_{\rm g},\alpha_{\rm g}$, but only 2 of these
are independent because we impose 4 constraints:

  1. 1.

The integral of the PDF over all $s$ is unity.

  2. 2.

The PDF is continuous at the transition point $s_{\rm g}$.

  3. 3.

The PDF is differentiable at $s_{\rm g}$ i.e., $dp/ds$ is continuous).

  4. 4.

The total mass is conserved, so $\int_{-\infty}^{\infty}e^{s}p(s)=1$.

In fitting this functional form to the PDFs produced by our simulations, we
leave the slope of the PL part ($\alpha_{\rm g}$) and the width of the LN part
($\sigma_{\rm s}$) as free parameters, and solve for the remaining parameters
from the 4 constraints; details of the calculation can be found in Appendix
A.1. It is important to note that the well-known relation $s_{0}=-\sigma_{\rm
s}^{2}/2$ (see, e.g., Federrath et al. 2008b) holds only for a pure LN PDF,
and does not hold for a LN+PL form.

We also introduce a double power-law form of the density PDF, which we show
below (Section 3.1) is a better fit to our simulations. This double power-law
form, to which we shall refer as LN+2PL, is a straightforward extension of the
LN+PL model, but now with two power-law sections,

$p(s)=\begin{cases}\frac{N}{\sqrt{2\pi\sigma_{\rm
s}^{2}}}\exp\left[-\frac{(s-s_{0})^{2}}{2\sigma_{\rm s}^{2}}\right]&s<s_{\rm
g},\\\ Np_{0}\exp\left(-s\alpha_{\rm g}\right)&s_{\rm g}\leq s\leq s_{\rm
d},\\\ Np_{1}\exp\left(-s\alpha_{\rm d}\right)&s_{\rm d}\leq s\leq s_{\rm
sink},\end{cases}$ (8)

where $\alpha_{\rm d}$ is the slope of the second PL tail, $s_{\rm d}$ is the
(log) density at which the first power law changes over to the second, $s_{\rm
sink}$ is a maximum cut-off density and all other symbols have the same
meanings as in the LN+PL case.111Note that, in contrast to the LN+PL form, for
LN+2PL we include an explicit truncation of the PDF at $s=s_{\rm sink}$; this
is necessary because, for many of the fits we will perform below, we obtain
$\alpha_{\rm d}<1$, in which case allowing the PDF to continue all the way to
$s=\infty$ causes the normalisation to diverge. This is not necessary for the
LN+PL models, because we can simply impose the constraint $\alpha_{\rm g}>1$,
and we find that in practice our best fitting values never approach this
limit. However, this difference in how we handle the upper limit matters only
for the normalisation, since in practice when we carry out fits to either the
LN+PL or LN+2PL forms, we only fit the data on cells with $s<s_{\rm sink}$.
This model has nine free parameters, $N$, $\sigma_{\rm s}$, $s_{0}$, $s_{\rm
g}$, $s_{\rm d}$, $p_{0}$, $p_{1}$, $\alpha_{\rm g}$, and $\alpha_{\rm d}$,
out of which five are constrained by requiring the function

  1. 1.

be normalized, $\int_{-\infty}^{s_{\rm sink}}p(s)ds=1$, where $s_{\rm sink}$
is the sink particle creation threshold,

  2. 2.

conserve mass, $\int_{-\infty}^{s_{\rm sink}}e^{s}p(s)ds=1$,

  3. 3.

be continuous at $s_{\rm g}$

  4. 4.

be continuous at $s_{\rm d}$, and

  5. 5.

be differentiable at $s_{\rm g}$.

We choose to fit the four parameters $\sigma_{\rm s}$, $\alpha_{\rm g}$,
$\alpha_{\rm d}$ and $s_{\rm d}$, while determining the others from the
constraints – see Appendix A.2.

The mechanics of fitting either the LN+PL or LN+2PL functions to a
distribution of cell densities from a simulation require some care. One option
would be to bin the simulation data in $s$, but this would produce an
undesirable dependence of the results on the choice of binning. We therefore
choose to fit using the cumulative distribution functions (CDFs)
$P(s)=\int_{-\infty}^{s}p(s^{\prime})\,ds^{\prime}$ rather than the
differential PDFs. In addition, to ensure that the PL tail(s) are fitted well,
we use

${}y(s)=\ln\left[\frac{P(s)}{1-P(s)}\right],$ (9)

rather than just $P(s)$ as our fitting function. When performing the fit, we
omit cells with $s>s_{\rm sink}$, since in these cells the density
distribution has been artificially altered by our sink particle method. Our
code to carry out this fit is available
publicly222https://github.com/shivankhullar/PDF_Fit.

###  2.3 Relation between density PDF and star formation efficiency

The star formation efficiency per free-fall time $\epsilon_{\rm ff}$ is a
useful tool that has been measured and studied extensively on the scales of
entire molecular clouds (Heiderman et al. 2010; Lada et al. 2010; Krumholz
2014; Federrath 2013; Salim et al. 2015; Onus et al. 2018; Pokhrel et al.
2021). Burkhart (2018) demonstrates the importance of the form of the PDF and
in particular its PL slope for the dimensionless rate of star formation
($\epsilon_{\rm ff}$).

For an entire cloud, $\epsilon_{\rm ff}\equiv\mbox{SFR}\times t_{\rm
ff}/M_{\rm gas}$ (Krumholz & McKee 2005; FK12), where SFR is the total star
formation rate inside the cloud, $t_{\rm ff}$ is the cloud’s free-fall time,
and $M_{\rm gas}$ is the total gas mass; intuitively, $\epsilon_{\rm ff}$ is
the ratio of the actual star formation rate to the rate that would be expected
if the gas were in free-fall collapse. Khullar et al. (2019) showed that one
can obtain a useful diagnostic for simulations by extending this definition to
consider only the material above some specified (log) density $s$,

$\epsilon_{\rm ff}(s)=\frac{{\rm SFR}\times t_{\rm ff}(s)}{M_{\rm
gas}(s)}=\frac{\mbox{SFR}\times t_{\rm ff,0}}{M_{\rm
tot}}\frac{\int_{s}^{\infty}e^{-s/2}p(s)ds}{\int_{s}^{\infty}p(s)\,ds\times\int_{s}^{\infty}e^{s}p(s)\,ds},$
(10)

where $t_{\rm ff}(s)$ and $M_{\rm gas}(s)$ are the mean free-fall time and
total mass of gas with density $>s$, and $t_{\rm
ff,0}=\sqrt{3\pi/32G\rho_{0}}$ and $M_{\rm tot}$ are the mean density free-
fall time and total gas mass, respectively.

Figure 1: Bottom: The density PDF described by equation 7, i.e., the
combination of a lognormal and power law (LN+PL case). The different colors
indicate PDFs having the same width of the LN part, $\sigma_{\rm s}=1.8$, but
different PL slopes $\alpha$. Star markers indicate the transition point,
$s_{\rm g}$, from the LN part to the PL part. Top: The corresponding
$\epsilon_{\rm ff}(s)$ curve found using equation 10. While it may be
difficult to distinguish the PL part from the tail of the LN part, the
$\epsilon_{\rm ff}(s)$ curves are very distinct and very sensitive to changes
in the PL slope.

In Figure 1, we show sample density PDFs $p(s)$ and their corresponding
$\epsilon_{\rm ff}(s)$ curves for LN+PL models with different PL slopes. We
normalize the $\epsilon_{\rm ff}$ curves by $\epsilon_{\rm
ff,0}=\lim_{s\to-\infty}\epsilon_{\rm ff}(s)$, which eliminates the constants
SFR, $t_{\rm ff,0}$, and $M_{\rm tot}$. Figure 1 shows that the slope of the
power-law tail of $\epsilon_{\rm ff}(s)$ is set by the density PDF slope
$\alpha_{\rm g}$. For $\alpha_{\rm g}=1.5$, we approach constant
$\epsilon_{\rm ff}(s)$ as $s\to\infty$. The PL slope can be related to a
density profile as $\rho(r)\propto r^{-3/\alpha_{\rm g}}$ (Kritsuk et al.
2011; Federrath & Klessen 2013). For $\alpha_{\rm g}=1.5$, $\rho(r)\propto
r^{-2}$ which is typically seen during the early stages of gravitational
collapse (Larson 1969; Penston 1969a; Shu 1977).

A PDF slope shallower than $1.5$ indicates slower than free-fall collapse, so
there is an accumulation of mass at high densities. For the PL part, the slope
of $\epsilon_{\rm ff}(s)$ is given as

$\frac{d\epsilon_{\rm ff}(s)}{ds}=\frac{(3/2-\alpha_{\rm g})\alpha_{\rm
g}(\alpha_{\rm g}-1)}{(\alpha_{\rm g}+1/2)}e^{(\alpha_{\rm g}+\frac{1}{2})s},$
(11)

which goes to zero when $\alpha_{\rm g}=1.5$.

##  3 Results

Throughout this section, we will analyse all results as a function of star
formation efficiency (SFE), in order to study the effects of cloud evolution.
The SFE is defined as

${\rm SFE}(t)=\frac{M_{*}(t)}{M_{*}(t)+M_{\rm gas}(t)},$ (12)

where $M_{*}(t)$ and $M_{\rm gas}(t)$ are the mass in sink particles and gas,
respectively, at any given time $t$. Note that in our simulations, the
denominator $M_{*}(t)+M_{\rm gas}(t)=M$, the (constant) total cloud mass.

###  3.1 Physical processes controlling the density PDF

Here we analyse the density structure, PDF and $\epsilon_{\rm ff}$ for
simulation v1M5 ($\alpha_{\rm vir,0}=1,\mathcal{M}=5$), to better understand
the basic physical processes that control the shape of the density PDF. In
Section 3.2 we analyse the dependence of the PDF on SFE and on the Mach number
and virial parameter of the clouds.

####  3.1.1 The density PDF $p(s)$ and its relation to the $\epsilon_{\rm
ff}(s)$ curve

As discussed in Section 2.3, the density PDF $p(s)$ can be used to calculate
the corresponding $\epsilon_{\rm ff}(s)$ curve. In Figure 2 we show the
density PDF (bottom panel) and the $\epsilon_{\rm ff}(s)$ curve (top panel) at
different times during the evolution of the cloud in the v1M5 simulation. At
$t=0$ when gravity is switched on, the shape of the PDF is close to a
lognormal (LN) distribution (Passot & Vázquez-Semadeni 1998; Klessen 2000;
Kritsuk et al. 2007; Federrath et al. 2008b, 2010a), typical of supersonic
turbulence. At SFE = $0$, just before the formation of the first sink
particle, the width of the LN part of the PDF has broadened and the PDF has
deviated enough from the LN to form a power-law (PL) tail; Jaupart & Chabrier
(2020) report this phenomenon as well. This tail has a slope close to
$\alpha=1.5$, as is evident from the flat portion in the $\epsilon_{\rm
ff}(s)$ curve between $s\sim 5$ and $s\sim 7$.

Figure 2: Dimensionless star formation efficiency $\epsilon_{\rm ff}(s)$ (top)
and density PDF $p(s)$ (bottom) as a function of logarithmic density $s$ for
the v1M5 simulation ($\alpha_{\rm vir,0}=1,\mathcal{M}=5$) at different times
during the evolution of the cloud. The different lines correspond to times
$t=0$ when gravity is switched on (blue line), SFE = 0%, just before the
formation of the first sink particle (green line), SFE = 5% (orange line), SFE
= 10% (red line) and SFE = 20% (cyan line). A LN+2PL fit (equation 8) to the
PDF at SFE = 5% is shown as the black dash-dotted line. The vertical dash-
dotted lines in the bottom panel indicate the transition points between the LN
part and the first PL, $s_{\rm g}$ (purple), from the first PL to the second
PL, $s_{\rm d}$ (green), and at which sink particles form (black). The shaded
regions correspond to the error bars in the fitted parameters.

At $\mbox{SFE}\gtrsim 5\%$, the PDF is relatively steady and exhibits only
minor variations in the slope of the PL tail(s). In this steady state,
$\epsilon_{\rm ff}(s)$ is fairly flat up to densities close to the sink
particle threshold, where it artificially rises for numerical reasons, as mass
is removed from the simulation grid and added to the sink particles. However,
superimposed on this mostly flat behaviour, there is a distinct bump and dip
structure, with a weak local maximum at $s\sim 5$ and a minimum at $s\sim 7$.
In terms of the density PDF, this corresponds to changes in the PL slope, from
a steeper PL between $s\sim 3-7$ to a shallower one at $s\gtrsim 7$.

We can demonstrate quantitatively that the a LN+2PL model is a better fit to
the simulation PDF than a LN+PL, at least once star formation begins, by
comparing the reduced $\chi^{2}$ goodness-of-fit statistics for the two
functional forms. We carry out LN+PL and LN+2PL fits to the density PDFs
measured for every snapshot of simulation v1M5, and we plot the ratio of the
reduced $\chi^{2}$ values, which we denote

$T\equiv\frac{\chi^{2}_{\rm red,LN+2PL}}{\chi^{2}_{\rm red,LN+PL}}$ (13)

as a function of SFE in the top panel of Figure 3. We find that $T<1$ for
almost the entire simulation, indicating that the LN+2PL model is a better
description of the simulation PDF. There are, however, exceptions: periods in
the simulation (e.g., at $\mbox{SFE}\approx 13-16\%$) when $T$ comes close to
or exceeds unity, indicating that a LN+PL fit is equivalent or preferred. At
these times, we find that the error bars on the second PL slope, $\alpha_{\rm
d}$, become very large, as illustrated in the lower panel of Figure 3.
Nonetheless, because the LN+2PL form is preferred at almost all times, we will
adopt this functional form for the purposes of all our analysis in the
remainder of this paper.

Figure 3:  Top: The quantity $T$, which is the ratio $\chi^{2}_{\rm
red,LN+2PL}/\chi^{2}_{\rm red,LN+PL}$, as a function of SFE for the simulation
v1M5. For large parts of the simulation, the LN+2PL model is a significantly
better fit than the LN+PL model. Bottom: The fitted values of the slope of the
second PL tail, $\alpha_{\rm d}$, with the vertical orange lines indicating
the error bars on the estimated values. When $\chi^{2}_{\nu,{\rm
2PL}}/\chi^{2}_{\nu,{\rm 1PL}}$ is close to 1, the LN+2PL density PDF
parameters have large error bars or are highly variable.

####  3.1.2 Gas structures at the power-law transition points $s_{\rm g}$ and
$s_{\rm d}$

Figure 4: Left: The density PDF for a snapshot from the simulation v1M5
($\alpha_{\rm vir,0}=1,\mathcal{M}\approx 5$) shown in blue, alongside our
LN+2PL model fit (red). The vertical lines and corresponding shaded regions
(errors in the fitted parameters) indicate the transition points from the LN
part to the first PL, $s_{\rm g}$ (purple), first PL to second PL, $s_{\rm d}$
(green), and sink particle creation threshold density, $s_{\rm sink}$ (black).
The change in the color-bar occurs at $s=s_{\rm g}$ to indicate different
regions in the gas flow. Centre: The volume-weighted column density,
$\rho_{\rm mw}$ (equation 14), for the same simulation snapshot. The colorbar
change occurs at $s=s_{\rm g}$. Right: A 0.1 pc box centred around the maximum
density in the simulation box (green square in the middle panel). We only show
gas that is denser the second transition point, $s_{\rm d}$. The gas denser
than $s_{\rm g}$ consists of dense cores and filamentary-like structures,
whereas gas denser than $s_{\rm d}$ is almost always found in rotating discs.
Figure 5: Virial parameter $\alpha_{\rm vir}(s)$ as a function of gas density
$s$ for simulation v1M5 at SFE=6%. The orange curve include the gravitational
potential energy of the gas alone, whereas the blue curve includes the
contributions of the stellar gravitational potential as well. As in Figure 4,
dash-dotted vertical lines and bands show the transition points and error bars
in the best-fitting density PDF, $s_{\rm g}$ (purple) and $s_{\rm d}$ (green),
and the sink particle threshold density $s_{\rm sink}$ (black). The horizontal
dashed line indicates $\alpha_{\rm vir}=1$, which separates unbound gas
($\alpha_{\rm vir}>1)$ from gas that is bound and collapsing ($\alpha_{\rm
vir}<1)$. Gas denser than $s_{\rm g}$ is bound and collapsing.

The LN+2PL model has two breakpoints: $s_{\rm g}$, where the LN part
transitions to the first PL (or, equivalently, the $\epsilon_{\rm ff}(s)$
curve begins to bend down), and $s_{\rm d}$, where the first power law
converts to the second, and the $\epsilon_{\rm ff}$ curve has its second
inflection point. We next attempt to understand the nature of these transition
points. We begin by looking at the column density morphology at different gas
densities in Figure 4, which shows the state of the v1M5 simulation at SFE =
6%. In the left panel we show the density PDF and the corresponding LN+2PL
model fit to the simulation data. The middle panel shows the mass-weighted
mean density of the gas through the simulation volume, given by

$\rho_{\rm mw}=\frac{\int\rho^{2}\,dx}{\int\rho\,dx},$ (14)

where $x$ is the line-of-sight distance. We show the mass-weighted density
projection because it picks out the dense gas structures along the line of
sight better than a volume-weighted projection. The change in the colourbar
from white-blue to black-red-yellow marks the transition point $s_{\rm g}$.

The right panel shows a zoomed-in 0.1 pc box centred around the maximum
density (green box in the middle panel) and only shows gas denser than $s_{\rm
d}$.

From the figure, it is immediately obvious that the higher density threshold,
$s_{\rm d}$ is picking out structures found in rotationally-flattened discs.
The connection between the second PL and rotationally-supported discs had
previously been made by Kritsuk et al. (2011) and Murray et al. (2017), and
our analysis is consistent with theirs. The meaning of $s_{\rm g}$ is somewhat
less obvious from the morphology, but can be clarified by examining the
physical state of the gas.

Specifically, we calculate the virial parameter $\alpha_{\rm vir}$ (dropping
the subscript to distinguish it from the virial ratio of the whole simulation
box, $\alpha_{\rm vir,0}$) as a function of gas (log over-)density, $s$. The
virial parameter is defined as the ratio of two terms in the virial theorem:
twice the kinetic energy, divided by the absolute value of the potential
energy: $\alpha_{\rm vir}=2E_{\rm kin}/|E_{\rm grav}|$. Both $E_{\rm kin}$ and
$E_{\rm grav}$ depend on the structure and distribution of the gas. To find
$\alpha_{\rm vir}(s)$, we extract connected regions above a density $s$, which
we decompose into clumps using yt’s Clump Finder333https://yt-
project.org/doc/analyzing/analysis_modules/clump_finding.html (Turk et al.,
2011). We compute the kinetic energy for each clump in its centre-of-mass
frame, and we compute the gravitational potential both for the gas alone, and
for the combined gravity of the gas and the stars.444For the gravitational
energy, we use the full potential due to all gas and stars within the
simulation, rather than just the self-gravity of the region in the clump. This
has the effect of suppressing $\alpha_{\rm vir}$ somewhat for clumps deep in
an external potential well (such as overdensities within an accretion disc).
However, this is a minor effect, which does not change our qualitative
conclusions.

We plot $\alpha_{\rm vir}(s)$ in Figure 5. The figure shows that $\alpha_{\rm
vir}>1$ for densities comparable to $\rho_{0}$ ($s\approx 0$), regardless of
whether we include the stellar contribution to the potential. As the density
increases, $\alpha_{\rm vir}$ decreases, reaching a minimum $\alpha_{\rm
vir}<1$ at $s\approx s_{\rm g}$, indicating that this gas is bound and
collapsing; at this transition, stellar gravity contributes roughly half the
total potential energy. As we go to even higher densities, $\alpha_{\rm vir}$
rises until it reaches $s\approx s_{\rm d}$. The rise is sharp if we neglect
stellar gravity, but only very shallow if we include it. It is therefore clear
that stellar gravity largely dictates the behaviour of the gas density PDF for
$s\geq s_{\rm g}$. At densities higher than $s_{\rm g}$, $\alpha_{\rm vir}$
increases slightly (while still being $<1$ throughout), primarily due to
increased turbulent and rotational support as we approach the disc density
$s_{\rm d}$. It reaches a maximum at $s\approx s_{\rm d}$, before flattening
out.

We are now in a position to interpret the features of the PDF: $s_{\rm g}$
marks the transition between bound and unbound gas, and at this transition gas
and stellar gravity are roughly equally important. For $s>s_{\rm g}$, stellar
gravity becomes increasingly dominant, such that the gas is not bound to
itself, but is bound to stars. At $s\approx s_{\rm d}$, rotational support
becomes increasingly important, such that gas at $s\gtrsim s_{\rm d}$ is found
largely in rotationally-supported discs.

###  3.2 Variation of the density PDF parameters

Having described our density PDF model and the physical processes controlling
it, we now look at how the parameters describing the PDF vary over time and as
a function of Mach number and (global) virial parameter.

####  3.2.1 Dependence of the PDF on SFE

Figure 6: The evolution of the various density PDF parameters: (a) the width
of the LN part of the PDF, $\sigma_{\rm s}$ (in red), (b) the two PL slopes,
$\alpha_{\rm g}$ and $\alpha_{\rm d}$ (blue and orange respectively), (c) the
two transition points, $s_{\rm g}$ and $s_{\rm d}$ (purple and green
respectively), as a function of Star Formation Efficiency (SFE) in different
panels for different simulation parameters. The heavy central line shows the
best-fitting value, and the vertical lines around them show the $1\sigma$
uncertainty. Various trends are discussed in detail in Section 3.2. The
parameters $\sigma_{\rm s}$ and $s_{\rm g}$ tend to increase with increasing
SFE, whereas the PL slopes $\alpha_{\rm g}$ and $\alpha_{\rm d}$ remain fairly
constant.

We fit our LN+2PL model described in Section 2 to every snapshot in each of
our simulations. In Figure 6, we show how the six important (four fitted, two
derived) PDF parameters vary with SFE for each of the simulations. The plot is
arranged so that each row has a constant $\alpha_{\rm vir,0}$ and varying
$\mathcal{M}$, whereas each column has constant $\mathcal{M}$ and varying
$\alpha_{\rm vir,0}$.

Starting with the width of the LN part of the PDF, $\sigma_{\rm s}$ (red
lines), we note that $\sigma_{\rm s}$ appears to be nearly constant throughout
the star formation process in most of the simulations. The only secular trend
visible is that $\sigma_{\rm s}$ increases slightly with SFE in the
simulations with $\alpha_{\rm vir,0}=0.5$. For these simulations, the
gravitational collapse tends to widen the distribution of the LN part of the
PDF, as the turbulence becomes increasingly gravity-driven and thus
compressive555If we consider the time between $t=0$ and $\mbox{SFE}=0\%$, all
the simulations show an increase in the width of the LN part of the PDF in
this duration.. This is consistent with the predictions of Jaupart & Chabrier
(2020); we discuss comparisons with this model, and extensions to it, further
below.

The PL slopes, $\alpha_{\rm g}$ (blue lines) and $\alpha_{\rm d}$ (orange
lines), remain fairly constant with SFE for all simulations, independent of
$\mathcal{M}$ or $\alpha_{\rm vir,0}$. As already shown in Figure 3, there
exist certain windows of time (or SFE) for which a LN+PL model fits the data
almost as well as our fiducial LN+2PL, and these are reflected in the
intermittently large error bars on the estimates of $\alpha_{\rm g}$ or
$\alpha_{\rm d}$. Nonetheless, it is clear that there are no secular trends in
$\alpha_{\rm g}$ or $\alpha_{\rm d}$ with SFE.

The transition points between the LN and the first PL, $s_{\rm g}$ (purple
lines), and the first PL and the second PL, $s_{\rm d}$ (green lines), exhibit
more interesting behaviour. In the $\mathcal{M}\approx 2.5$ and
$\mathcal{M}\approx 5$ simulations, considering the range of SFEs where
$s_{\rm g}$ is well defined (i.e., ignoring times when $s_{\rm g}$ has large
error bars because the PDF is ambiguous between LN+PL and LN+2PL), we find
that $s_{\rm g}$ tends to increase with SFE. This is especially true for the
$\alpha_{\rm vir,0}=0.5$ cases, since gravity plays a more important role for
them. This trend is not surprising given that $s_{\rm g}=s_{0}+\alpha_{\rm
g}\sigma_{\rm s}^{2}$ (see equation 26), and that $\sigma_{\rm s}$ increases
slightly with SFE, whereas $\alpha_{\rm g}$ and $s_{0}$ do not show a
systematic increase or decrease with SFE. Our estimation of the break point
between the first and second PL, $s_{\rm d}$, reveals that for the
$\mathcal{M}\approx 2.5$ and $\mathcal{M}\approx 5$ cases, there are
relatively large fluctuations, which make it hard to identify any trends
across the range of SFEs probed here. For the $\mathcal{M}\approx 10$
simulations, our resolution is not sufficient to resolve the discs in detail,
especially for the $\alpha_{\rm vir,0}=0.5$ and $\alpha_{\rm vir,0}=2$ cases
(see Appendix B for more details on the resolution criterion). Consequently,
the uncertainties on $s_{\rm g}$ and $s_{\rm d}$ are very large throughout the
simulations, which precludes us from analysing these cases further.

Figure 7: The variation of the PDF parameters on the simulation parameter
$\mathcal{M}$, keeping $\alpha_{\rm vir,0}$ constant. (Top row:) $\sigma_{\rm
s}$ (red) and $s_{0}$ (cyan), (Middle row:) $\alpha_{\rm g}$ (blue) and
$\alpha_{\rm d}$ (orange), (Bottom row:) $s_{\rm g}$ (purple) and $s_{\rm d}$
(green), with the Mach number $\mathcal{M}$, for 3 different virial parameters
($\alpha_{\rm vir,0}=0.5$ (left column), $\alpha_{\rm vir,0}=1$ (centre
column), $\alpha_{\rm vir,0}=2$ (right column). The different shapes of the
points indicate the average over an SFE range; plus shapes represent the
average between SFE=1% and 3%, pentagons between SFE=4% and 6%, squares
between SFE=9% and 11%, crosses between SFE=14% and 16% (or <16% for
simulations stopped before then), and lastly, circles represent the average
over the entire SFE range. The red dashed line and the dash-dotted line in the
top and bottom row respectively represent our fit to equation 17 and equation
18 for each of the 3 different cases $\alpha_{\rm vir,0}=0.5,1,2$ (see Section
4.1). The dashed line in the bottom row shows the global fit to our model for
the disc density (see Section 4.2, equation 22). The shaded regions indicate
the 2$\sigma$ error bars for each of our fits. Various trends are discussed in
Section 3.2.2.  Figure 8: The variation of the PDF parameters on the
simulation parameter $\alpha_{\rm vir,0}$, keeping $\mathcal{M}$ constant.
(Top row:) $\sigma_{\rm s}$ (red) and $s_{0}$ (cyan), (Middle row:)
$\alpha_{\rm g}$ (blue) and $\alpha_{\rm d}$ (orange), (Bottom row:) $s_{\rm
g}$ (purple) and $s_{\rm d}$ (green), with the simulation virial parameter
$\alpha_{\rm vir,0}$, for 3 different Mach numbers ($\mathcal{M}\approx 2.5$
(left column), $\mathcal{M}\approx 5$ (centre column), $\mathcal{M}\approx 10$
(right column). The different shapes of the points indicate the average over
an SFE range; plus shapes represent the average between SFE=1% and 3%,
pentagons between SFE=4% and 6%, squares between SFE=9% and 11%, crosses
between SFE=14% and 16% (or <16% for simulations stopped before then), and
lastly, circles represent the average over the entire SFE range. The red
dashed line and the dash-dotted line in the top and bottom row respectively
represent our global fit to equation 17 and equation 18 for all the
simulations (see Section 4.1). The dashed line in the bottom row shows the
global fit to our model for the disc density (see Section 4.2, equation 22).
The shaded regions indicate the 2$\sigma$ error bars for each of our fits.
Various trends are discussed in Section 3.2.2.

####  3.2.2 Dependence of the PDF on the Mach number and virial parameter

We now look at the dependence of the fitted PDF parameters on $\mathcal{M}$
and $\alpha_{\rm vir,0}$, since ultimately, any efforts to systematically
study the density PDF should enable us to understand how its shape varies as a
function of quantities measurable through observations. For this purpose, for
each simulation we average the density PDFs for all snapshots over three
intervals of SFE: 1% – 3%, 9% – 11%, and 14% – 16% (or to the highest SFE
attained, for simulations that do not reach 16%). We then fit LN+2PL
functional forms to these averaged PDFs, extract the best-fitting parameters,
and plot the results in Figure 7 and Figure 8; these two figures show the same
underlying data, but the former shows variation with $\mathcal{M}$ at fixed
$\alpha_{\rm vir,0}$, while the latter shows variation with $\alpha_{\rm
vir,0}$ at fixed $\mathcal{M}$.

First consider the dependence on $\mathcal{M}$. The top row of Figure 7 shows
that the width and the peak of the LN part of the PDF, $\sigma_{\rm s}$ and
$s_{0}$, clearly depend on $\mathcal{M}$. This is to be expected given that,
in the absence of self-gravity, the former obeys the well-known relationship
(Federrath et al., 2008b, 2010a; Price et al., 2011; Konstandin et al., 2012;
Molina et al., 2012; Federrath & Banerjee, 2015; Nolan et al., 2015;
Kainulainen & Federrath, 2017)

$\sigma_{\rm s}^{2}=\ln(1+b^{2}\mathcal{M}^{2}).$ (15)

While this relationship predicts the correct qualitative trend of $\sigma_{\rm
s}$ with $\mathcal{M}$, quantitatively equation 15 systematically under-
predicts the measured value of $\sigma_{\rm s}$. To take one example, run
V0.5M2.5 (top left panels in Figure 7 and Figure 8) has $\mathcal{M}=2.5$, and
for our natural driving mix, we expect $b\approx 0.4$; plugging in, the
predicted value of $\sigma_{\rm s}\approx 0.8$, whereas the values that we
obtain by fitting the observed PDF are $\approx 1.1-1.2$, roughly 30% larger.
We will return to the source of this undershoot below.

The peak of the PDF, $s_{0}$, also varies with $\mathcal{M}$, which is again
expected since conservation of mass requires $s_{0}=-\sigma_{\rm s}^{2}/2$ if
the PDF is a pure LN. For our LN+2PL model, $s_{0}$ does not have a
straightforward dependence on $\sigma_{\rm s}$ (see Section A.2), but is still
strongly anti-correlated with it, and consequently with $\mathcal{M}$.

The middle row of Figure 7 shows the two PL slopes, $\alpha_{\rm g}$ and
$\alpha_{\rm d}$, as a function of $\mathcal{M}$. From the figure it is
difficult to identify any clear trends with $\mathcal{M}$, especially since
the uncertainties on $\alpha_{\rm g}$ and $\alpha_{\rm d}$ are relatively
large. This is because, while the LN+2PL model is a better fit to the data
most of the time, there are times when a single PL is a better description of
the density PDF, and during these intervals the values of the slopes become
highly uncertain.

In the bottom row of Figure 7 we show the dependence of $s_{\rm g}$ and
$s_{\rm d}$ on $\mathcal{M}$. We see a clear trend for $s_{\rm g}$ to increase
with $\mathcal{M}$, also seen in Burkhart & Mocz (2019), which is expected
given that $s_{\rm g}\propto\sigma_{\rm s}^{2}$ (see Appendix A.2). The
increased width of the LN for higher $\mathcal{M}$ pushes $s_{\rm g}$ to
higher $s$. Since $\alpha_{\rm g}$ does not show systematic trends, the trend
in $s_{\rm g}$ with $\mathcal{M}$ is determined solely by how $s_{0}$ and
$\sigma_{\rm s}$ vary. We return to this variation and the physical reasons
for it in Section 4.1. We also find that $s_{\rm d}$ increases with
$\mathcal{M}$, with the exception of the $\mathcal{M}\approx 10$ case for
which the separation of the two scales $s_{\rm g}$ and $s_{\rm d}$ is limited
by numerical resolution constraints (see Appendix B).

We now consider the dependence of the fitted PDF parameters on the simulation
virial parameter, $\alpha_{\rm vir,0}$, in Figure 8. Once again, the top panel
shows how $\sigma_{\rm s}$ and $s_{0}$ vary with $\alpha_{\rm vir,0}$, the
middle panel shows the variation of $\alpha_{\rm g}$ and $\alpha_{\rm d}$, and
the bottom panel shows the variation of $s_{\rm g}$ and $s_{\rm d}$ with
$\alpha_{\rm vir,0}$. There is a weak anti-correlation of $\sigma_{\rm s}$
with $\alpha_{\rm vir,0}$ such that increasing $\alpha_{\rm vir,0}$ decreases
$\sigma_{\rm s}$. The reason for this is easy to understand: in the
low-$\alpha_{\rm vir,0}$ cases, the dynamics are determined mostly by self-
gravity, and this results in an increase in the width of the LN part of the
density PDF because it makes the turbulence more compressive (as discussed
above in the context of Figure 6 above). The dependence of $s_{0}$ on
$\alpha_{\rm vir,0}$ can be understood similarly: $s_{0}$ increases with
increasing $\alpha_{\rm vir,0}$ since $s_{0}$ and $\sigma_{\rm s}$ are anti-
correlated. The virial parameter $\alpha_{\rm vir,0}$ has no significant
effect on the PL slopes, $\alpha_{\rm g}$ and $\alpha_{\rm d}$. We do not
expect any such dependence, since they are set by how the process of collapse
proceeds in self-gravitating gas, and this is fairly independent from the
varying strength of turbulence relative to self-gravity on larger scales.
However, we caution that our ability to draw any strong conclusions is once
again limited by the large error bars on our fitted PL slopes.

Lastly, the transition points, $s_{\rm g}$ and $s_{\rm d}$, do not show any
strong dependence on $\alpha_{\rm vir,0}$, except a weak increase with
$\alpha_{\rm vir,0}$ for the well-resolved cases (i.e., except for
$\mathcal{M}=10$). Burkhart & Mocz (2019) expect a weak dependence of $s_{\rm
g}$ on $\alpha_{\rm vir,0}$. We will seek to explain these trends in the next
section.

##  4 Physical models for the PDF parameters

Several authors have proposed models to relate cloud properties to density PDF
parameters. For instance, the slopes of the PL tails, $\alpha_{\rm g}$ (or
$\alpha_{\rm d}$) can be related to the density profile of the collapsing gas
cloud as $\rho(r)\propto r^{-3/\alpha_{\rm g}}$ (e.g., Federrath & Klessen
2013; Kritsuk et al. 2011). In this section we shall seek to explain the
trends observed in Section 3.2.2 for the width of the LN part and the values
of the two transition points.

###  4.1 Lognormal width and first powerlaw transition

As noted earlier, in the presence of gravity, the width of the LN part of the
density PDF, $\sigma_{\rm s}$, increases, as does the log density $s_{\rm g}$
at which the LN gives way to powerlaw behaviour. As discussed above, the
former trend is qualitatively consistent with the relationship between
$\sigma_{\rm s}$ and $\mathcal{M}$ that has long been known for non-self-
gravitating turbulence (equation 15). The latter trend is also qualitatively
consistent with the recent model proposed by Jaupart & Chabrier (2020), who
predict a relationship

$|e^{s_{\rm g}}-1|=(b\mathcal{M})^{2}\times\alpha_{\rm
vir,0}\times\left|\frac{s_{\rm g}+\frac{1}{2}\sigma_{\rm s}^{2}}{\sigma_{\rm
s}^{2}}\right|.$ (16)

However, quantitative evaluation of this prediction shows a result similar to
that for $\sigma_{\rm s}$: the predicted value of $s_{\rm g}$ is substantially
smaller than the measured one. Again using run V0.5M2.5 as an example (bottom
left panel of Figure 7 and Figure 8), the value of $s_{\rm g}$ predicted by
equation 16 is 0.46, whereas the measured value is, depending on the SFE,
between 1 and 2.666This discrepancy can be slightly reduced if we fit to LN+PL
functional forms, as assumed by Jaupart & Chabrier (2020), rather than the
LN+2PL ones that we have shown are a better description of the data. Doing so
yields fits for $s_{\rm g}$ that are typically slightly lower than the values
obtained from the LN+2PL fit, and thus in better agreement with the
predictions of equation 16. However, a significant discrepancy still remains
even for LN+PL fits.

Fortunately, both of these discrepancies can be explained by a simple
observation: the value of $b$ that we have estimated for turbulence with a
natural mix of modes is probably not applicable to a system where self-gravity
makes an appreciable contribution to the driving of turbulence, which in turn
renders the flow more compressible. Indeed, a number of authors have already
pointed out that, locally, the driving parameter can vary and change the star
formation rate (see e.g. Menon et al. 2020; Körtgen 2020; Menon et al. 2021).
We therefore consider a modification of equation 15 in which we relate the
width of the LN part of the PDF with the cloud parameters as

$\sigma_{\rm s}^{2}=\ln(1+b_{\rm eff}^{2}\mathcal{M}^{2})$ (17)

where $b_{\rm eff}$ is an effective driving parameter in the presence of
gravity, which is to be determined. Similarly, in the model of Jaupart &
Chabrier (2020), the transition point between the LN and the first PL tail
also depends on $b$ according to equation 16. We therefore suggest an
analogous modification to this relationship,

$|e^{s_{\rm g}}-1|=(b_{\rm eff}\mathcal{M})^{2}\times\alpha_{\rm
vir,0}\times\left|\frac{s_{\rm g}+\frac{1}{2}\sigma_{\rm s}^{2}}{\sigma_{\rm
s}^{2}}\right|,$ (18)

where here $\sigma_{\rm s}$ is understood to be that given by equation 17, to
include our effective driving parameter.

To investigate whether this single-parameter modification is sufficient to
explain the discrepancies between the models and the data, we perform a non-
linear least-squares fit for $b_{\rm eff}$ using the values of $\sigma_{\rm
s}$ and $s_{\rm g}$ predicted by equation 17 and equation 18, respectively,
and those obtained from the LN+2PL fits to our simulations; for the purposes
of this fit, we weight errors in $\sigma_{\rm s}$ and $s_{\rm g}$ equally, so
the best fit is the value of $b_{\rm eff}$ that minimises the sum of the
squared difference between the predicted and observed values of $\sigma_{\rm
s}$ and $s_{\rm g}$. We carry out this fit in two ways; first, we fit
different values of $b_{\rm eff}$ for the three cases $\alpha_{\rm
vir,0}=0.5,1,2$, and second we fit all simulations together; we refer to this
latter case as a global fit. We show the results of the first fitting method
in Figure 7, and the second in Figure 8. (Recall that the underlying data are
the same in the two figures; they are simply organised differently.) We report
the best-fitting values of $b_{\rm eff}$ for both cases in Table 2.

Examining the figures, it is clear that either the global or the case-by-case
fits yield predicted values of $\sigma_{\rm s}$ and $s_{\rm g}$ that are in
reasonably good agreement with the measurements, and capture both the correct
mean values and the dependence on $\mathcal{M}$; the model-predicted
dependence on $\alpha_{\rm vir}$ is also consistent with the simulation
measurements, although in this case our limited dynamic range in $\alpha_{\rm
vir}$ and the noisiness of the simulation results means that no dependence on
$\alpha_{\rm vir}$ at all would probably be almost an equally-good model.
Nonetheless, this experiment demonstrates that our simple procedure of
modifying the value of $b$ seems to repair most of the problems in the models.
Table 2 shows that the values of $b_{\rm eff}$ requires for this purpose are
significantly greater than the $b=0.4$ expected for a natural mix of
compressive and solenoidal modes, and are substantially closer to the $b=1$
expected for purely compressive turbulence. In the case where we allow
different values of $b_{\rm eff}$ for different $\alpha_{\rm vir}$, the
effective driving parameter shows a clear anti-correlation with the virial
parameter of the simulation. As the relative strength of gravity increases
($\alpha_{\rm vir,0}$ decreases), the effective driving parameter increases,
consistent with the idea that, in these simulations, gravity is providing a
substantial contribution to the turbulent velocity, and that it drives
predominantly compressive modes. We acknowledge, however, that this is just
one possible explanation for the discrepancies between the models and the
data, one that we have adopted largely on the basis of simplicity and the
principle of making minimal modifications to existing models. Other
modifications to equation 15 are certainly plausible, and with our limited and
noisy sampling of parameter space, we are not in a position to rule out such
approaches.

Case | $b_{\rm eff}$
---|---
$\alpha_{\rm vir,0}=0.5$ | 0.95$\pm$0.04
$\alpha_{\rm vir,0}=1.0$ | 0.77$\pm$0.02
$\alpha_{\rm vir,0}=2.0$ | 0.70$\pm$0.01
Global fit | 0.82$\pm$0.02
Table 2: Best-fitting effective turbulence driving parameter $b_{\rm eff}$
from equation 17 and equation 18, obtained using the data in Figure 7. In the
top three rows, we show fits obtained by separately fitting the simulations
with $\alpha_{\rm vir}=0.5,1$, and 2, as indicated; these 3 cases correspond
to 3 columns in Figure 7. The bottom row, Global fit, shows the value of
$b_{\rm eff}$ obtained by fitting for all the simulations simultaneously.

###  4.2 Disc transition scale

Having determined the scalings of $\sigma_{\rm s}$ and $s_{\rm g}$ with cloud
properties, we next seek to do the same for the characteristic density of disc
formation, $s_{\rm d}$. For this purpose we adapt an argument from Burkert &
Bodenheimer (2000), in which motions at the onset of collapse are assumed to
sample the background line width-size relation. We estimate the scale of a
collapsing region to be the sonic scale of the turbulence cascade. Approximate
angular momentum conservation then yields a characteristic disc scale, and the
importance of gravitational instability determines the outer disc density.
Similar estimates appear in previous works by Matzner & Levin (2005) and
Kratter & Matzner (2006), and this calculation is a minor extension of the one
presented by Burkhart & Mocz (2019).

We begin with the line-width size relation

$\sigma_{V}(r)=\sigma_{c}\left(\frac{r}{R_{c}}\right)^{\beta},$ (19)

where $R_{c}$ is the driving scale for turbulent motions,
$\sigma_{c}=\sigma_{V}(R_{c})=\mathcal{M}c_{s}$. Across the supersonic
turbulent cascade, $\beta\approx 1/2$, and this suffices for our purposes
despite the fact that the exponent tends toward Kolmogorov value of
($\beta=1/3$) around the sonic transition (Federrath et al., 2021). Here
$r_{s}$ is the sonic scale, at which $\sigma_{V}(r_{s})=c_{s}$; therefore

$r_{s}=\frac{R_{c}}{\mathcal{M}^{1/\beta}}.$ (20)

To collapse, a region on this scale must achieve a gas density $\rho_{s}$
above the critical density required for collapse:

$r_{s}^{2}=\frac{\phi_{1}c_{s}^{2}}{G\rho_{s}}$ (21)

where $\phi_{1}$ is a factor of order unity. This implies a critical sonic
mass given by $m_{s}=\phi_{2}\rho_{s}r_{s}^{3}$, where we account for geometry
and complicated dynamics with the factor $\phi_{2}$. In the spirit of Matzner
& Levin (2005) and Kratter & Matzner (2006), we define the disc radius as
$r_{d}=\phi_{d}r_{s}$. (Here $\phi_{d}$ is significantly smaller than unity,
and reflects the typical turbulent angular momentum of a region on the brink
of collapse.)

The disc density can be written as
$\rho_{d}=\phi_{3}\Sigma_{d}\Omega/c_{s}=\phi_{3}\Omega^{2}/(\pi GQ)$, where
$\Sigma_{d}$ is the surface density, $\Omega=\phi_{4}Gm_{s}/r_{d}^{3}$ is the
orbital speed and $Q=c_{s}\Omega/(\pi G\Sigma_{d})$ is the Toomre (1964)
stability parameter. Using our definition for the virial parameter (eq. 5),
and adopting $R_{c}=L/2$, we obtain

$e^{s_{\rm d}}=\phi_{w}\frac{\alpha_{\rm
vir,0}\mathcal{M}^{\frac{2}{\beta}-2}}{Q}$ (22)

where $\phi_{w}$ is a collection of the other coefficients, defined as
$24\phi_{1}\phi_{2}\phi_{3}\phi_{4}/(5\pi\phi_{d}^{3})$. We expect $Q\approx
1-2$ so that gravitational instability stimulates the disc angular momentum
transport.

While every one of these dimensionless quantities may depend on the
dimensionless physical and numerical parameters of our simulations (and
presumably does), we hypothesize that the dependence is weak so long as the
disc scale is reasonably well resolved. In this spirit, we perform a one
parameter fit of our model to the estimates of $s_{\rm d}$ from the
simulations. For this, we fix $\beta=1/2$ and allow $\phi_{w}$ to be a free
parameter. We find that a global fit to all our simulations results in
$\phi_{w}=62.0\pm 7.6$. This would mean that the best fit value for the
product $\phi_{1}\phi_{2}\phi_{3}\phi_{4}\approx 0.13$, which is reasonable.
We show the predictions using this best-fitting model in Figure 7 and Figure
8. While there is a significant amount of scatter (at least in part because
$s_{\rm d}$ is not always well-determined in our simulations), we see that the
model captures both the absolute value of $s_{\rm d}$ and its dependence on
$\alpha_{\rm vir}$ and $\mathcal{M}$ reasonably well.

##  5 Conclusions

Self-gravitating turbulence is a common description of the flow in star-
forming molecular clouds, and a common paradigm for theories of star
formation. However, we still lack a complete theory for the density
distributions produced in such media. In this work we characterise the density
distribution in self-gravitating, non-magnetised, turbulent media by
performing numerical simulations with a range of Mach numbers and virial
parameters. Using these simulations, we show that:

  1. 1.

The volumetric density PDF can be divided into three regimes: turbulence
dominated, gravity dominated, disc/rotation dominated. These regimes are
characterised by different functional forms: the PDF transitions from being a
lognormal (LN) at low densities (turbulence-dominated regime) to a first
powerlaw (PL; gravity-dominated) and a second (typically shallower) powerlaw
(disc regime) at higher densities. Our LN+2PL model includes elements of the
LN+PL model by Burkhart (2018) and Burkhart & Mocz (2019), in particular, the
importance of the sonic scale in setting a PL transition. However, contrary to
their results, we find for the most part that the PDF does not evolve with
star formation efficiency (SFE), and instead remains nearly steady once the
SFE reaches a few percent.

  2. 2.

While the break in the powerlaw (PL) slope can be difficult to see directly in
the PDF, statistical tests confirm that a fit containing two powerlaws is
almost always preferred over one with a single powerlaw, and the presence of a
second powerlaw becomes obvious when we examine the star formation efficiency
$\epsilon_{\rm ff}(s)$ as a function of log density $s$ (Khullar et al.,
2019). In this diagnostic, the break between the gravity-dominated and disc-
dominated regimes manifests as a characteristic oscillatory pattern whereby
$\epsilon_{\rm ff}(s)$ reaches a local maximum in the gravity-dominated regime
and then a local minimum at the transition from this regime to the disc
regime. This oscillatory pattern is present in all of our simulations at
almost all times once star formation begins.

  3. 3.

The gravity dominated part of the PDF is characterised by gas with local
virial parameter $\alpha_{\rm vir}<1$ (in agreement with observational results
of Chen et al. 2018, 2019). This is a result of stellar gravity being more
important as we go from $s_{\rm g}$, the log density at which gravity becomes
dominant and where gas self-gravity and stellar gravity are about equally
important, to $s_{\rm d}$, the log density at which rotational support becomes
dominant, where only stellar gravity is relevant. At densities above $s_{\rm
d}$ rotational support kicks in, leaving $\alpha_{\rm vir}\sim 1$, and the
morphology becomes disc-like.

  4. 4.

The width of the LN part of the PDF obeys increases with Mach number, and
scale at which the LN gives way to a gravity-dominated PL increases with both
Mach number and virial parameter. The rates of increase are consistent with
prior theoretical predictions, provided we make one substantial modification
to these models: self-gravitating turbulence is significantly more compressive
than turbulence driven with a natural mix of modes, which has the effect of
increasing both the LN width and the transition density. Models predicting
these two parameters yield reasonably accurate results only if we set the
turbulence-driving parameter $b$ that appears in them to a value close to that
expected for purely compressive turbulence, rather than one that reflects the
properties of the turbulence that is intially present at the onset of star
formation.

  5. 5.

We show that the disc formation scale, and its variation with Mach number and
virial parameter, can be understood as arising from the characteristic amount
of angular momentum at the sonic scale of the turbulence.

We therefore conclude that in a turbulent self-gravitating medium, the nature
of gas flows are captured in the shape of the density PDF and that a LN+2PL
model can be used to describe this density PDF well. With the help of this
model, we can determine the effects of the molecular cloud properties on the
density PDF, a necessary step toward formulating a predictive theory of star
formation. However, we caution that we have not yet included the effects of
stellar feedback or magnetic fields, which would likely modify the functional
form of the PDF. We leave these steps to future work, and caution against
applying the conclusions of our study directly to observed molecular clouds
until we are able to carry them out.

## Acknowledgements

We would like to thank the anonymous referee for their comments. SK would like
to thank Norman Murray and Piyush Sharda for valuable discussions. The
research of SK and CDM is supported by an NSERC Discovery Grant. MRK
acknowledges support from the Australian Research Council through its Future
Fellowship and Discovery Projects funding schemes, awards FT180100375 and
DP190101258. CF acknowledges funding provided by the Australian Research
Council (Future Fellowship FT180100495), and the Australia-Germany Joint
Research Cooperation Scheme (UA-DAAD). Computations were performed on the
Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the
Canada Foundation for Innovation; the Government of Ontario; Ontario Research
Fund - Research Excellence; and the University of Toronto. We further
acknowledge high-performance computing resources provided by the Australian
National Computational Infrastructure (grant ek9) in the framework of the
National Computational Merit Allocation Scheme and the ANU Merit Allocation
Scheme, and by the Leibniz Rechenzentrum and the Gauss Centre for
Supercomputing (grant pr32lo). The simulation software FLASH was in part
developed by the DOE-supported Flash Center for Computational Science at the
University of Chicago.

## Data availability

The simulation data underlying this article will be shared on reasonable
request to the corresponding author.

## References

  * Alves et al. (2017) Alves J., Lombardi M., Lada C. J., 2017, A&A, 606, L2
  * Ballesteros-Paredes et al. (1998) Ballesteros-Paredes J., Vazquez-Semadeni E., Scalo J., 1998, Astrophys. J., 515, 286
  * Ballesteros-Paredes et al. (2011) Ballesteros-Paredes J., Hartmann L. W., Vázquez-Semadeni E., Heitsch F., Zamora-Avilés M. A., 2011, Mon. Not. R. Astron. Soc., 411, 65
  * Bate & Burkert (1997) Bate M. R., Burkert A., 1997, MNRAS, 288, 1060
  * Berger & Colella (1989) Berger M. J., Colella P., 1989, Journal of Computational Physics, 82, 64
  * Berkhuijsen & Fletcher (2008) Berkhuijsen E. M., Fletcher A., 2008, Mon. Not. R. Astron. Soc. Lett., 390, 19
  * Brunt et al. (2010a) Brunt C. M., Federrath C., Price D. J., 2010a, MNRAS, 403, 1507
  * Brunt et al. (2010b) Brunt C. M., Federrath C., Price D. J., 2010b, MNRAS, 405, L56
  * Burkert & Bodenheimer (2000) Burkert A., Bodenheimer P., 2000, ApJ, 543, 822
  * Burkhart (2018) Burkhart B., 2018, ApJ, 863, 118
  * Burkhart & Lazarian (2012) Burkhart B., Lazarian A., 2012, Astrophys. J. Lett., 755, 19
  * Burkhart & Mocz (2019) Burkhart B., Mocz P., 2019, ApJ, 879, 129
  * Burkhart et al. (2010) Burkhart B., Stanimirovic S., Lazarian A., Kowal G., 2010, Astrophys. J., 708, 1204
  * Burkhart et al. (2015) Burkhart B., Collins D. C., Lazarian A., 2015, Astrophys. J., 808, 48
  * Burkhart et al. (2016) Burkhart B., Stalpes K., Collins D. C., 2016, Astrophys. J., 834, L1
  * Burkhart et al. (2017) Burkhart B., Stalpes K., Collins D. C., 2017, ApJ, 834, L1
  * Chen et al. (2018) Chen H. H.-H., Burkhart B., Goodman A., Collins D. C., 2018, The Astrophysical Journal, 859, 162
  * Chen et al. (2019) Chen H. H.-H., et al., 2019, ApJ, 877, 93
  * Collins et al. (2011) Collins D. C., Padoan P., Norman M. L., Xu H., 2011, ApJ, 731, 59
  * Collins et al. (2012) Collins D. C., Kritsuk A. G., Padoan P., Li H., Xu H., Ustyugov S. D., Norman M. L., 2012, Astrophys. J., 750, 13
  * Cunningham et al. (2018) Cunningham A. J., Krumholz M. R., McKee C. F., Klein R. I., 2018, MNRAS, 476, 771
  * Dib et al. (2007) Dib S., Kim J., Vázquez-Semadeni E., Burkert A., Shadmehri M., 2007, ApJ, 661, 262
  * Dib et al. (2020) Dib S., et al., 2020, Astron. Astrophys., 642, 177
  * Donkov & Stefanov (2018) Donkov S., Stefanov I. Z., 2018, MNRAS, 474, 5588
  * Donkov & Stefanov (2019) Donkov S., Stefanov I., 2019, MNRAS, 485, 3224
  * Donkov et al. (2020) Donkov S., Stefanov I. Z., Veltchev T. V., Klessen R. S., 2020, MNRAS, 000, 1
  * Dubey et al. (2008) Dubey A., et al., 2008, Astronomical Society of the Pacific Conference Series, 385, 145
  * Elmegreen (2002a) Elmegreen B. G., 2002a, Astrophys. J., 564, 773
  * Elmegreen (2002b) Elmegreen B. G., 2002b, ApJ, 577, 206
  * Eswaran & Pope (1988) Eswaran V., Pope S. B., 1988, Computers and Fluids, 16, 257
  * Federrath (2013) Federrath C., 2013, MNRAS, 436, 3167
  * Federrath (2015) Federrath C., 2015, MNRAS, 450, 4035
  * Federrath & Banerjee (2015) Federrath C., Banerjee S., 2015, MNRAS, 448, 3297
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Federrath & Klessen (2013) Federrath C., Klessen R. S., 2013, ApJ, 763, 51
  * Federrath et al. (2008a) Federrath C., Klessen R. S., Schmidt W., 2008a, ApJ, 688, L79
  * Federrath et al. (2008b) Federrath C., Klessen R. S., Schmidt W., 2008b, ApJ, 688, L79
  * Federrath et al. (2010a) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010a, A&A, 512, A81
  * Federrath et al. (2010b) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010b, ApJ, 713, 269
  * Federrath et al. (2014) Federrath C., Schrön M., Banerjee R., Klessen R. S., 2014, ApJ, 790, 128
  * Federrath et al. (2016a) Federrath C., et al., 2016a, ApJ, 832, 143
  * Federrath et al. (2016b) Federrath C., et al., 2016b, Astrophys. J., 832, 143
  * Federrath et al. (2021) Federrath C., Klessen R. S., Iapichino L., Beattie J. R., 2021, Nature Astronomy,
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Ginsburg et al. (2013) Ginsburg A., Federrath C., Darling J., 2013, ApJ, 779, 50
  * Girichidis et al. (2014) Girichidis P., Konstandin L., Whitworth A. P., Klessen R. S., 2014, Astrophys. J., 781, 91
  * Grisdale et al. (2019) Grisdale K., Agertz O., Renaud F., Romeo A. B., Devriendt J., Slyz A., 2019, MNRAS, 486, 5482
  * Grudić et al. (2018) Grudić M. Y., Hopkins P. F., Lee E. J., Murray N., Faucher-Giguère C.-A., Johnson L. C., 2018, Monthly Notices of the Royal Astronomical Society, 488, 1501
  * Heiderman et al. (2010) Heiderman A., Evans II N. J., Allen L. E., Huard T., Heyer M., 2010, ApJ, 723, 1019
  * Hennebelle & Chabrier (2008) Hennebelle P., Chabrier G., 2008, Astrophys. J., 684, 395
  * Hennebelle & Chabrier (2009) Hennebelle P., Chabrier G., 2009, Astrophys. J., 702, 1428
  * Hennebelle & Chabrier (2011) Hennebelle P., Chabrier G., 2011, ApJ, 743, L29
  * Heyer et al. (2016) Heyer M., Gutermuth R., Urquhart J. S., Csengeri T., Wienen M., Leurini S., Menten K., Wyrowski F., 2016, A&A, 588, A29
  * Hill et al. (2008) Hill A. S., Benjamin R. A., Kowal G., Reynolds R. J., Haffner L. M., Lazarian A., 2008, Astrophys. J., 686, 363
  * Hopkins (2013a) Hopkins P. F., 2013a, Mon. Not. R. Astron. Soc., 430, 1653
  * Hopkins (2013b) Hopkins P. F., 2013b, MNRAS, 430, 1880
  * Hunter (1977) Hunter C., 1977, ApJ, 218, 834
  * Jaupart & Chabrier (2020) Jaupart E., Chabrier G., 2020, ApJ, 903, L2
  * Kainulainen & Federrath (2017) Kainulainen J., Federrath C., 2017, A&A, 608, L3
  * Kainulainen et al. (2014) Kainulainen J., Federrath C., Henning T., 2014, Science (80-. )., 344, 183
  * Khullar et al. (2019) Khullar S., Krumholz M. R., Federrath C., Cunningham A. J., 2019, MNRAS, 488, 1407
  * Klessen (2000) Klessen R. S., 2000, ApJ, 535, 869
  * Konstandin et al. (2012) Konstandin L., Girichidis P., Federrath C., Klessen R. S., 2012, Astrophys. J., 761, 149
  * Körtgen (2020) Körtgen B., 2020, MNRAS, 497, 1263
  * Kratter & Matzner (2006) Kratter K. M., Matzner C. D., 2006, MNRAS, 373, 1563
  * Kritsuk et al. (2007) Kritsuk A. G., Norman M. L., Padoan P., Wagner R., 2007, Astrophys. J., 665, 416
  * Kritsuk et al. (2011) Kritsuk A. G., Norman M. L., Wagner R., 2011, ApJ, 727, L20
  * Krumholz (2014) Krumholz M. R., 2014, Phys. Rep., 539, 49
  * Krumholz (2015) Krumholz M. R., 2015, arXiv e-prints, p. arXiv:1511.03457
  * Krumholz & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
  * Krumholz et al. (2004) Krumholz M. R., McKee C. F., Klein R. I., 2004, ApJ, 611, 399
  * Krumholz et al. (2012) Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71
  * Lada et al. (2010) Lada C. J., Lombardi M., Alves J. F., 2010, ApJ, 724, 687
  * Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  * Lemaster & Stone (2008) Lemaster M. N., Stone J. M., 2008, Astrophys. J., 682, L97
  * Leroy et al. (2017) Leroy A. K., et al., 2017, ApJ, 846, 71
  * Lin et al. (2016) Lin Y., et al., 2016, ApJ, 828, 32
  * Lombardi et al. (2015) Lombardi M., Alves J., Lada C. J., 2015, A&A, 576, L1
  * Maier et al. (2016) Maier E., Chien L.-H., Hunter D. A., 2016, AJ, 152, 134
  * Maier et al. (2017) Maier E., Elmegreen B. G., Hunter D. A., Chien L.-H., Hollyday G., Simpson C. E., 2017, AJ, 153, 163
  * Marinkova et al. (2020) Marinkova L., Veltchev T. V., Donkov S., 2020, Publications of the Astronomical Society “Rudjer Boskovic”, 20, 51
  * Masunaga & Inutsuka (2000) Masunaga H., Inutsuka S.-i., 2000, ApJ, 531, 350
  * Mathew & Federrath (2020) Mathew S. S., Federrath C., 2020, MNRAS, 496, 5201
  * Matzner & Levin (2005) Matzner C. D., Levin Y., 2005, ApJ, 628, 817
  * Menon et al. (2020) Menon S. H., Federrath C., Kuiper R., 2020, MNRAS, 493, 4643
  * Menon et al. (2021) Menon S. H., Federrath C., Klaassen P., Kuiper R., Reiter M., 2021, MNRAS, 500, 1721
  * Mocz et al. (2017) Mocz P., Burkhart B., Hernquist L., McKee C. F., Springel V., 2017, ApJ, 838, 40
  * Molina et al. (2012) Molina F. Z., Glover S. C., Federrath C., Klessen R. S., 2012, Mon. Not. R. Astron. Soc., 423, 2680
  * Murray & Chang (2015) Murray N., Chang P., 2015, ApJ, 804, 44
  * Murray et al. (2017) Murray D. W., Chang P., Murray N. W., Pittman J., 2017, MNRAS, 465, 1316
  * Myers (2015) Myers P. C., 2015, ApJ, 806, 226
  * Myers et al. (2014) Myers A. T., Klein R. I., Krumholz M. R., McKee C. F., 2014, MNRAS, 439, 3420
  * Nolan et al. (2015) Nolan C. A., Federrath C., Sutherland R. S., 2015, Mon. Not. R. Astron. Soc., 451, 1380
  * Nordlund & Padoan (1999) Nordlund Å. K., Padoan P., 1999, in Franco J., Carraminana A., eds, Interstellar Turbulence. p. 218 (arXiv:astro-ph/9810074)
  * Onus et al. (2018) Onus A., Krumholz M. R., Federrath C., 2018, MNRAS, 479, 1702
  * Ostriker et al. (2001) Ostriker E. C., Stone J. M., Gammie C. F., 2001, ApJ, 546, 980
  * Padoan & Nordlund (2002) Padoan P., Nordlund A., 2002, Astrophys. J., 617, 559
  * Padoan & Nordlund (2004) Padoan P., Nordlund Å., 2004, ApJ, 617, 559
  * Padoan & Nordlund (2011) Padoan P., Nordlund A., 2011, The Astrophysical Journal, 730, 40
  * Padoan et al. (1997) Padoan P., Nordlund A., Jones B. J. T., 1997, MNRAS, 288, 145
  * Padoan et al. (2007) Padoan P., Nordlund A., Kritsuk A. G., Norman M. L., Li P. S., 2007, Astrophys. J., 661, 972
  * Pan et al. (2018) Pan L., Padoan P., Nordlund Å., 2018, ApJ, 866, L17
  * Pan et al. (2019) Pan L., Padoan P., Nordlund Å., 2019, ApJ, 881, 155
  * Passot & Vázquez-Semadeni (1998) Passot T., Vázquez-Semadeni E., 1998, Phys. Rev. E - Stat. Physics, Plasmas, Fluids, Relat. Interdiscip. Top., 58, 4501
  * Penston (1969a) Penston M. V., 1969a, MNRAS, 144, 425
  * Penston (1969b) Penston M. V., 1969b, MNRAS, 145, 457
  * Pokhrel et al. (2016) Pokhrel R., et al., 2016, Mon. Not. R. Astron. Soc., 461, 22
  * Pokhrel et al. (2021) Pokhrel R., et al., 2021, ApJ, 912, L19
  * Pope & Ching (1993) Pope S. B., Ching E. S., 1993, Phys. Fluids A, 5, 1529
  * Price & Federrath (2010) Price D. J., Federrath C., 2010, Mon. Not. R. Astron. Soc., 406, 1659
  * Price et al. (2011) Price D. J., Federrath C., Brunt C. M., 2011, ApJ, 727, L21
  * Ricker (2008) Ricker P. M., 2008, ApJS, 176, 293
  * Salim et al. (2015) Salim D. M., Federrath C., Kewley L. J., 2015, ApJ, 806, L36
  * Schmidt et al. (2009) Schmidt W., Federrath C., Hupp M., Kern S., Niemeyer J. C., 2009, A&A, 494, 127
  * Schneider et al. (2013) Schneider N., et al., 2013, ApJ, 766, L17
  * Schneider et al. (2015a) Schneider N., et al., 2015a, Mon. Not. R. Astron. Soc. Lett., 453, L41
  * Schneider et al. (2015b) Schneider N., et al., 2015b, Astron. Astrophys., 575, 79
  * Schneider et al. (2015c) Schneider N., et al., 2015c, Astron. Astrophys., 578, 29
  * Shu (1977) Shu F. H., 1977, ApJ, 214, 488
  * Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
  * Tremblin et al. (2014) Tremblin P., et al., 2014, A&A, 564, A106
  * Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., 1997, ApJ, 489, L179
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9
  * Vazquez-Semadeni (1994) Vazquez-Semadeni E., 1994, ApJ, 423, 681
  * Veltchev et al. (2019) Veltchev T. V., Girichidis P., Donkov S., Schneider N., Stanchev O., Marinkova L., Seifried D., Klessen R. S., 2019, Mon. Not. R. Astron. Soc., 489, 788
  * Vutisalchavakul et al. (2016) Vutisalchavakul N., Evans Neal J. I., Heyer M., 2016, ApJ, 831, 73
  * Waagan et al. (2011) Waagan K., Federrath C., Klingenberg C., 2011, Journal of Computational Physics, 230, 3331
  * Whitworth & Summers (1985) Whitworth A., Summers D., 1985, MNRAS, 214, 1

##  Appendix A The LN+PL and LN+2PL models

In this section, we describe our LN+PL and LN+2PL models in more detail,
including all the constraints and the relations between the parameters, and
how we use these constraints when we carry out fits.

###  A.1 LN+PL

Starting with the definition of our function, $p(s)$ (see equation 7), we
first require that the integral over all $s$ be normalized to unity. This
constrains the normalization parameter $N$ to be

$N=\frac{1}{\frac{p_{0}e^{-\alpha_{\rm g}s_{\rm g}}}{\alpha_{\rm
g}}+\frac{1}{2}\text{erfc}\left(\frac{s_{0}-s_{\rm g}}{\sqrt{2}\sigma_{\rm
s}}\right)}.$ (23)

Next, we apply the constraint for differentiability at the transition point,
$s_{\rm g}$. The derivatives on either side of $s_{\rm g}$ are

$p^{\prime}(s)=\begin{cases}-\frac{(s-s_{0})e^{-\frac{(s-s_{0})^{2}}{2\sigma_{\rm
s}^{2}}}}{\sqrt{2\pi}\sigma_{\rm s}^{3}\left(\frac{p_{0}e^{-\alpha_{\rm
g}s_{\rm g}}}{\alpha_{\rm g}}+\frac{1}{2}\text{erfc}\left(\frac{s_{0}-s_{\rm
g}}{\sqrt{2}\sigma}\right)\right)}&s-s_{\rm g}<0\\\ -\frac{\alpha_{\rm
g}p_{0}e^{-\alpha_{\rm g}s}}{\frac{p_{0}e^{-\alpha_{\rm g}s_{\rm
g}}}{\alpha_{\rm g}}+\frac{1}{2}\text{erfc}\left(\frac{s_{0}-s_{\rm
g}}{\sqrt{2}\sigma}\right)}&s-s_{\rm g}\geq 0.\end{cases}$ (24)

These two expressions are equal if the parameter $p_{0}$ obeys

$p_{0}=-\frac{(s_{0}-s_{\rm g})e^{\alpha_{\rm g}s_{\rm g}-\frac{(s_{\rm
g}-s_{0})^{2}}{2\sigma_{\rm s}^{2}}}}{\sqrt{2\pi}\alpha_{\rm g}\sigma_{\rm
s}^{3}}.$ (25)

Applying the continuity constraint and substituting for $p_{0}$, we get

$s_{\rm g}=s_{0}+\alpha_{\rm g}\sigma_{\rm s}^{2},$ (26)

which is similar to the expression derived by Burkhart et al. (2017), except
that we cannot use the relation $s_{0}=-1/2\sigma_{\rm s}^{2}$ that holds for
pure LN PDFs, since we impose a constraint on mass conservation in our model.
This final constraint reads $\int_{-\infty}^{\infty}e^{s}p(s)ds=1$. Thus we
can solve for $s_{0}$, under the assumption that $\sigma_{\rm s}>0$ and
$\alpha_{\rm g}>1$, as

$s_{0}=\log\left(\frac{A}{B}\right)$ (27)

where

$\displaystyle A$ $\displaystyle=$ $\displaystyle(\alpha_{\rm
g}-1)\left[\sqrt{\pi}\alpha_{\rm g}\sigma_{\rm s}e^{\frac{\alpha_{\rm
g}^{2}\sigma_{\rm s}^{2}}{2}}\left(\text{erf}\left(\frac{\alpha_{\rm
g}\sigma_{\rm s}}{\sqrt{2}}\right)+1\right)+\sqrt{2}\right],$ (28)
$\displaystyle B$ $\displaystyle=$ $\displaystyle\sqrt{2}\alpha_{\rm
g}e^{\alpha_{\rm g}\sigma_{\rm s}^{2}}$ (29)
$\displaystyle{}+\sqrt{\pi}(\alpha_{\rm g}-1)\alpha_{\rm g}\sigma_{\rm
s}e^{\frac{1}{2}\left[\alpha_{\rm g}^{2}+1\right)\sigma_{\rm
s}^{2}}\left(\text{erf}\left(\frac{(\alpha_{\rm g}-1)\sigma_{\rm
s}}{\sqrt{2}}\right)+1\right].$ (30)

When fitting to LN+PL distributions, our procedure is to leave $\alpha_{\rm
g}$ and $\sigma_{\rm s}$ as free parameters. We then derive $s_{0}$ from
equation 27, $s_{\rm g}$ from equation 26, $p_{0}$ from equation 25, and
finally $N$ from equation 23, which completes specification of the model.

###  A.2 LN+2PL

Starting with the definition of our function $p(s)$, defined in equation 8, we
impose the normalization constraint first. This gives,

$\displaystyle N$ $\displaystyle=$
$\displaystyle\left[\frac{p_{0}}{\alpha_{\rm g}}(e^{-\alpha_{\rm g}s_{\rm
g}}-e^{-\alpha_{\rm g}s_{\rm d}})+\frac{p_{1}}{\alpha_{\rm d}}(e^{-\alpha_{\rm
d}s_{\rm d}}-e^{-\alpha_{\rm d}s_{\rm sink}})\right.$ (31)
$\displaystyle\left.{}+\frac{1}{2}{\rm erfc}\left(\frac{s_{0}-s_{\rm
g}}{\sqrt{2}\sigma_{\rm s}}\right)\right]^{-1}$

Next, we apply the differentiability constraint at the transition points
$s_{\rm g}$, which gives us

$p_{0}=\frac{s_{\rm g}-s_{0}}{\sqrt{2\pi}\sigma_{\rm s}^{3}\alpha_{\rm g}}{\rm
exp}\left(\alpha_{\rm g}s_{\rm g}-\frac{(s_{\rm g}-s_{0})^{2}}{2\sigma_{\rm
s}^{2}}\right).$ (32)

As for the single PL case, we next apply the continuity constraint at $s_{\rm
g}$, which leads once again to

$s_{\rm g}=s_{0}+\alpha_{\rm g}\sigma_{\rm s}^{2}.$ (33)

We cannot require differentiability at the second transition point $s_{\rm d}$
as that would imply $\alpha_{\rm g}=\alpha_{\rm d}$. However, we do require
continuity at $s_{\rm d}$, which gives us

$p_{1}=p_{0}e^{-(\alpha_{\rm g}-\alpha_{\rm d})s_{\rm d}}.$ (34)

Finally, we apply the mass conservation constraint. Evaluating the integral
over the various parts of the PDF, this can be expressed as

$\displaystyle\frac{N}{2}(J+K+L)=1$ (35)

where

$\displaystyle J$ $\displaystyle=$ $\displaystyle\frac{2p_{0}}{\alpha_{\rm
g}-1}[e^{(s_{\rm g}-\alpha_{\rm g}s_{\rm g})}-e^{(s_{\rm d}-\alpha_{\rm
g}s_{\rm d})}],$ (36) $\displaystyle K$ $\displaystyle=$
$\displaystyle\frac{2p_{1}}{\alpha_{\rm d}-1}[e^{(s_{\rm d}-\alpha_{\rm
d}s_{\rm d})}-e^{(s_{\rm sink}-\alpha_{\rm d}s_{\rm sink})}],$ (37)
$\displaystyle L$ $\displaystyle=$ $\displaystyle{\rm
exp}\left(s_{0}+\frac{\sigma_{\rm s}^{2}}{2}\right)\times{\rm
erfc}\left(\frac{s_{0}+\sigma_{\rm s}^{2}-s_{\rm g}}{\sqrt{2}\sigma_{\rm
s}}\right).$ (38)

When fitting LN+2PL functions, we leave $\sigma_{\rm s}$, $\alpha_{\rm g}$,
$\alpha_{\rm d}$, and $s_{\rm d}$ as our four free parameters. Since the
equations above constitute a transcendental system for the remaining
parameters, we must solve them numerically. Fortunately, we can reduce the
problem to a one-dimensional root-finding exercise as follows: for any set of
the four fitting parameters and a trial value of $s_{0}$, we can use equation
33 to obtain $s_{\rm g}$, equation 32 to obtain $p_{0}$, equation 34 to obtain
$p_{1}$, and then equation 31 to solve for $N$. In general the value of $N$
obtained thereby will not satisfy the mass constraint given by equation 35. We
therefore iteratively adjust $s_{0}$ until this condition is satisfied, at
which point we will have found the correct, converged values for $s_{0}$,
$s_{\rm g}$, $p_{0}$, $p_{1}$, and $N$.

##  Appendix B The disc sizes and resolution criterion

Figure 9: Each row represents simulations with 3 different Mach numbers,
$\mathcal{M}\approx 2.5$ (top), $\mathcal{M}\approx 5$ (middle),
$\mathcal{M}\approx 10$ (bottom), but the same virial parameter, $\alpha_{\rm
vir,0}=1$. Left column: The density PDF for a single snapshot shown in blue,
alongside our LN+2PL model fit (red). The vertical lines and corresponding
shaded regions (errors in the fitted parameters) indicate the transition
points from the LN part to the first PL, $s_{\rm g}$ (purple), first PL to
second PL, $s_{\rm d}$ (green), and sink particle creation threshold density,
$s_{\rm sink}$ (black). The change in the color-bar occurs at $s=s_{\rm g}$ to
indicate different regions in the gas flow. Centre column: The volume-weighted
column density, $\rho_{\rm mw}$ (equation 14), for the same simulation
snapshots. The colorbar change occurs at $s=s_{\rm g}$. Right column: A 0.2 pc
box ($\mathcal{M}\approx 2.5$) and 0.1 pc boxes ($\mathcal{M}\approx
5,\mathcal{M}\approx 10$) centred around the maximum density in the simulation
boxes (green squares in the middle panel). We only show gas that is denser
than the second transition point, $s_{\rm d}$. The gas denser than $s_{\rm g}$
consists of dense cores and filamentary-like structures, whereas gas denser
than $s_{\rm d}$ is almost always found in rotating discs. The disc sizes
decrease as the $\mathcal{M}$ number increases. Figure 10: Top: $\epsilon_{\rm
ff}(s)$ curves for the simulations with $\alpha_{\rm
vir,0}=2,\mathcal{M}\approx 5$ (case v2M5), but differing effective
resolutions: 1611 AU (blue), 805 AU (red), 402 AU (green), 201 AU (orange),
101 AU (cyan). The curves are shown for an SFE=5%, but are time-averaged over
the range 4-6%. The vertical dashed lines indicate the sink particle density
threshold for the corresponding simulations. The slope of the $\epsilon_{\rm
ff}(s)$ curve depends on the PL tail slope in the density PDF and therefore
can be used as a proxy for estimating the resolution required to resolve the
second PL tail and disc-like structures. Bottom: $\epsilon_{\rm ff}(s)$ curves
at different SFE for the simulations with $\alpha_{\rm
vir,0}=1,\mathcal{M}\approx 10$ (v1M10) but differing effective resolutions:
201 AU (orange and cyan), and 50 AU (red and teal). Once again, SFE=5%
represents a time-average over the SFE range 4-6%, and SFE=13% is a time-
average over the range 12-14%. The $\Delta x\approx 50$ AU simulation curves
start to show indications of the rise and fall structure that is indicative of
a second PL tail, but do not quite reach the resolution that would be required
to capture it.

In this appendix we investigate the resolution required for a simulation to
capture the formation of discs and the transition from the first to second PL.
Figure 9 is similar to Figure 4, and is constructed following the same
procedure; however, in Figure 9 each of the three rows shows a simulation with
a different Mach number (top: $\mathcal{M}\approx$ 2.5, middle:
$\mathcal{M}\approx$ 5, bottom: $\mathcal{M}\approx$ 10) but the same virial
parameter, $\alpha_{\rm vir,0}=1$. The right-most column of Figure 9 shows the
gas denser than the second transition point $s_{\rm d}$, which lies almost
entirely in discs. There is a clear dependence of the disc size on the Mach
number, which is not unexpected, since in equation 20 we show that disc sizes
are a fraction of the sonic scale, $r_{\rm d}\propto\mathcal{M}^{-2}$. This
implies that a disc in a $\mathcal{M}\approx 5$ simulation is about 4 times
larger than in a $\mathcal{M}\approx 10$ simulation, completely consistent
with the ratio of disc sizes visible in Figure 9.

We next look at the resolution required to resolve the discs and as a result
the second PL tail. As shown in Section 2.3, the $\epsilon_{\rm ff}(s)$ curve
is very sensitive to the slope of the PL tails of the density PDF. As a
result, a second PL tail in $p(s)$ having a PL slope less than 1.5 can be seen
in the $\epsilon_{\rm ff}(s)$ curve as the range where the slope,
$d\epsilon_{\rm ff}(s)/ds$, is negative (see e.g., Figure 1). In the absence
of a second PL tail, we do not expect a change in slope at high $s$ in the
$\epsilon_{\rm ff}(s)$ curve. Since gas in the second PL tail is present in
discs, we do not expect the second PL tail to be resolvable when the discs are
unresolved. To test this, we repeat simulation v2M5 with five different
maximum resolutions: 1611 AU, 805 AU, 402 AU, 201 AU, and 101 AU; the 201 AU
case is the standard one that we use elsewhere in the paper. In Figure 10
(top), we show $\epsilon_{\rm ff}(s)$ for these 5 simulations; the curves are
averaged over the SFE range 4% and 6% for each individual simulation. The
blue, red and green curves, which belong to simulations with a minimum cell
size $\Delta x\approx 1611,805$ and $402$ AU, respectively, are not able to
resolve the discs, a failure that manifests as the absence of the
characteristic local maximum and minimum that we have shown is associated with
the development of the second PL tail. By contrast, there is a clear
indication of a second PL tail in the orange and cyan curves, which correspond
to resolutions of 201 and 101 AU, respectively. These simulations are able to
resolve the discs and capture the gas dynamics more accurately. The disc size
seen in this simulation has a radius $r_{\rm d}\sim$ 2000 AU. Therefore we
estimate that capturing the disc / second PL transition requires that
$\mathcal{H}\equiv r_{\rm d}/\Delta x\gtrsim 10$.

To check that this condition matches other simulations, and also whether the
resolution requirement is time-dependent, we also repeat the v1M10 simulation
at a resolution of $\Delta x=201$ AU, a factor of 4 lower than our fiducial
resolution of 50 AU for this case. The bottom panel of Figure 10 shows the
resulting $\epsilon_{\rm ff}(s)$ curves at SFEs of 5% (averaged over SFEs of
4-6%) and 13% (averaged over 12-14%). As expected, a minimum cell size of
$\Delta x\approx 201$ AU is not sufficient to resolve the discs in this case.
For $\Delta x\approx 50$ AU resolution, we begin to see a flattening of the
$\epsilon_{\rm ff}(s)$ curve at $s>5$, indicating the first hints of discs and
the formation of a PL slope $\alpha_{\rm d}$ close to the asymptotic value of
1.5. This is consistent with our resolution criterion of $\mathcal{H}\gtrsim
10$, because in this case we find disc radius $\sim 300$ AU. Thus our maximum
resolution reaches $\mathcal{H}\approx 6$, close to but not quite reaching our
requirement $\mathcal{H}\gtrsim 10$.