The Effects of Magnetic Fields and Outflow Feedback on the Shape and Evolution of the Density Probability Distribution Function in Turbulent Star-forming Clouds

Star formation takes place in dense and cold giant molecular clouds (GMCs) that are subject to magnetized supersonic turbulent motions (e.g., Padoan et al. 1997; McKee & Ostriker 2007; Lazarian 2009; Kennicutt & Evans 2012; Myers et al. 2014; Krumholz et al. 2018). Star formation is, in part, controlled by nonlinear fluid dynamics, and models of star formation focus on understanding the roles that self-gravity, magnetohydrodynamic (MHD) turbulence, and galactic environment play in the formation of dense collapsing gas that may form stars (e.g., Collins et al. 2012; Padoan et al. 2012; Burkhart et al. 2017). For individual star-forming clouds, a common approach to modeling the relevant physics is to investigate the distribution of gas via a density probability distribution function (PDF) analysis (e.g., Krumholz & McKee 2005; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart 2018; see Padoan et al. 2014 for a review). Models for star formation that use a density PDF have been used to explain a range of phenomena including the mass distribution of cores and the stellar initial mass function (e.g., Hennebelle & Chabrier 2008, 2009, 2011; Hopkins 2012). Such models can also be used as subgrid models in galaxy formation simulations to model the dependence of local star formation on the dynamical state of the gas (Braun & Schmidt 2015; Semenov et al. 2016, 2021; Li et al. 2017; Trebitsch et al. 2017; Lupi et al. 2018; Gensior et al. 2020; Kretschmer & Teyssier 2020; Kretschmer et al. 2022; Olsen et al. 2021).

The shape of the density PDF in GMCs is expected to be lognormal when isothermal supersonic turbulence dominates the gas dynamics (see, e.g., Passot & Vázquez-Semadeni 1998, 2003; Klessen 2000) and to transition to a power law as gravitational contraction takes over at high densities (Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2013; Jaupart & Chabrier 2020; Khullar et al. 2021). This shape has been seen in both observations (i.e., in column density) and in numerical simulations (both in density and projected column density; see, e.g., Vazquez-Semadeni & Garcia 2001; Wada & Norman 2007; Ossenkopf-Okada et al. 2016; Veltchev et al. 2019). The lognormal form of the column density PDF describes the behavior of diffuse H i and ionized gas as well as some molecular clouds that are not forming massive stars (Hill et al. 2008; Burkhart et al. 2010; Kainulainen & Tan 2013; Schneider et al. 2015a). The dense gas in molecular clouds (i.e., extinctions greater than A_v > 1, corresponding to n > 10³ cm⁻³) is predominantly found to have a power-law PDF rather than a lognormal PDF (Collins et al. 2012; Girichidis et al. 2014; Myers 2015; Schneider et al. 2015b; Alves et al. 2017; Burkhart et al. 2017; Chen et al. 2018, 2019; Kainulainen & Federrath 2017; Mocz et al. 2017; Myers 2017). We note that when the medium is not isothermal, a bimodal density PDF distribution is observed (see, e.g., Vázquez-Semadeni et al. 2000; de Avillez & Breitschwerdt 2004; Gazol et al. 2005; Ballesteros-Paredes et al. 2011); however, for studies of molecular clouds, an isothermal approximation is reasonable given the cloud temperature does not strongly vary beyond T = 10–15 K.

Most analytic models of star formation have used only a lognormal form for the density PDF, where the width of the distribution is set by the properties of MHD turbulence, including the sonic Mach number and the type of driving (e.g., Krumholz & McKee 2005; Federrath & Klessen 2012; Salim et al. 2015). The star formation rates (SFRs) are obtained by integrating the PDF past a predefined critical density for collapse. ⁵ However, models that use a lognormal density PDF lack a time-varying treatment of stellar feedback (e.g., stellar winds, ionizing radiation, radiation pressure, and collimated outflows or jets) and neglect the evolution of the PDF shape under the influence of self-gravity.

To remedy this, such models modify the SFR by an efficiency factor, such as with a constant parameter ₀ or ${\epsilon }_{\mathrm{core}}$ that is due to stellar feedback, since feedback processes, along with magnetic fields, are now understood to be key ingredients for the low efficiency of star formation (see, e.g., Rosen et al. 2014, 2020; Federrath 2015; Grudić et al. 2018; Rosen & Krumholz 2020). As such, analytical models of star formation that assume a lognormal form for the density PDF are only able to take into account the important aspects of feedback in simple ad hoc ways that reduce the overall SFRs and efficiencies down to the observed values. Additionally, because the lognormal form is set only by turbulence, it neglects the effects of self-gravity on the density distribution, e.g., the power-law tail. ⁶

Motivated by these studies, Burkhart (2018) presented an analytic model for star formation in a gravo-turbulent medium based on a piecewise lognormal and power-law density PDF that can allow for mass cycling. The model presented in Burkhart (2018) and Burkhart & Mocz (2019) accounts for self-gravity and stellar feedback by including a power-law tail on the high-density end of the density PDF and allowing the form of the PDF to be time variable. This model estimates the star formation efficiency (SFE) by taking the ratio of the mass of star-forming gas in the power-law tail to the total gas in the cloud. The model has thus far been tested on MHD gravo-turbulent simulations that lacked stellar feedback (Burkhart & Mocz 2019; Khullar et al. 2021).

Few studies have focused on understanding how mass transfers between different parts of the density PDF and how the overall shape of the PDF is affected by the inclusion of stellar feedback. In this work, we analyze the density PDFs of a suite of hydrodynamical simulations of star formation within a molecular cloud that progressively adds physics starting from only self-gravity, then adding turbulence, magnetic fields, and finally including stellar feedback from protostellar jets and slower but wider disk winds. ⁷ Collimated protostellar jets are magnetically launched by the winding-up of the magnetic field in the accretion disk, following a magneto-centrifugal mechanism (Blandford & Payne 1982; Shu et al. 1988; Pelletier & Pudritz 1992; Bontemps et al. 1996; Lynden-Bell 2003; Maud et al. 2015; Kölligan & Kuiper 2018; Rosen et al. 2020). Jets and disk winds are well known to play a significant role in setting the structure of their host molecular cloud and in setting SFRs (see, e.g., Federrath et al. 2014; Rosen et al. 2020). In simulations, protostellar outflows are often included via subgrid models, such as the model described in Federrath et al. (2014), and used in this paper.

Using these simulations, we investigate how the density PDFs of star-forming molecular clouds evolve when we separately study the effects of gravity, turbulence, magnetic fields, and protostellar outflows. Our goal is to understand how different physical processes influence the density PDF and the movement of material from diffuse gas into stars. We also analyze how different portions of the density PDF change over time and how the mass flow of the gas connects to the change in the SFR with time.

This paper is organized as follows: In Section 2, we describe the setup of the simulations used in this paper, and in Section 3, we review analytic models for density PDFs. In Section 4, we discuss our fits to the simulated density PDFs and detail how mass flows between different sections of the PDF and into stars. In Section 5, we discuss the implications of our results, and in Section 6, we summarize our work.

For our analysis, we use a suite of four simulations that model a star-forming region within a molecular cloud. These four simulations include runs drawn from Federrath (2015, 2016a), and an entirely new simulation that includes protostellar outflows. Our simulations include varying physical processes, such as self-gravity, turbulence, magnetic fields, and protostellar outflows. All four simulations were performed using the FLASH hydrodynamics code, a publicly available adaptive mesh hydrodynamics code (Fryxell et al. 2000). FLASH solves the fully compressible MHD equations using adaptive mesh refinement (AMR) and can include many inter-operable modules. The simulations presented here use a Godunov-type method with the second-order, five-wave approximate HLL5R Riemann solver (Waagan et al. 2011).

The first simulation (Gravity) includes only self-gravity with an initial Gaussian random density distribution and zero initial velocities. This is the gravity only (Gaussian initial conditions) simulation from Federrath (2016a). The second simulation (Turbulence) initializes with uniform density and drives turbulence for two eddy-turnover times (t_turnover = L/(2σ_v ) ≈ 0.98 Myr, where L = 2 pc is the box size and σ_v = 1 km s⁻¹ is the velocity dispersion on the driving-scale L/2), before turning on self-gravity. Turbulence is continually driven at half the box size throughout the run to maintain the sonic Mach number close to ${{ \mathcal M }}_{s}=5$ and with a natural mix of forcing modes, such that the resulting turbulence driving parameter is b ∼ 0.4 (see the turbulence driving method developed in Federrath et al. 2010b). The third simulation (B-Fields) is identical to the Turbulence run, but adds magnetic fields with an Alfvén Mach number of ${{ \mathcal M }}_{{\rm{A}}}=2$. The magnetic field is initialized with uniform strength and oriented along the z-direction, as described in Federrath (2015). However, the turbulent stirring during the driving phrase randomizes the magnetic field orientation so that when we turn on gravity, the magnetic field is no longer uniform. Both the Turbulence and B-Fields runs were first presented in Federrath (2015). The fourth simulation (All + Outflows) is first presented in this paper. The All + Outflows simulation is similar to the B-Fields run but adds stellar feedback in the form of protostellar outflows and protostellar heating feedback, as described in Federrath et al. (2017). Table 1 lists the physical parameters of each simulation, and Figure 1 shows density projections for each of these four simulations.

Zoom In Zoom Out Reset image size

Table 1. Summary of Key Simulation Parameters

Simulation	Gravity?	σv (km s−1)	b	${ \mathcal M }$ s	B (μG)	${{ \mathcal M }}_{{\rm{A}}}$	Outflows?	ρt (g cm−3)	${N}_{\mathrm{res}}^{3}$	References
Gravity	Yes	N/A	N/A	N/A	N/A	N/A	No	N/A	10243	Federrath (2016a)
Turbulence	Yes	1	0.4	5	N/A	N/A	No	1.64 × 10−20	10243	Federrath (2015)
B-Fields	Yes	1	0.4	5	10	2.0	No	1.64 × 10−20	10243	Federrath (2015)
All + Outflows	Yes	1	0.4	5	10	2.0	Yes	1.64 × 10−20	20483	First presented here
No Gravity	No	1	0.4	5	N/A	N/A	No	1.64 × 10−20	10243	First presented here

Note. A summary of the simulations and key parameters: σ_v is the velocity dispersion, b is the turbulent driving parameter, ${ \mathcal M }$ _s is the sonic Mach number, B is the strength of the magnetic fields, ${{ \mathcal M }}_{{\rm{A}}}$ is the Alfvénic Mach number, and ρ_t is the calculated reference transition density. ${N}_{\mathrm{res}}^{3}$ shows the maximum effective grid resolution, and we list the reference for the paper in which each simulation was first presented. Note that the simulation with outflows also includes protostellar heating feedback. The initial mean volume density of each simulation is ρ₀ = 3.28 × 10⁻²¹ g cm ⁻³.

Download table as:  ASCIITypeset image

We ran an additional simulation without gravity (No Gravity, not shown in Table 1) to demonstrate the effect of self-gravity on the density PDF. This simulation has the same initial conditions as the Turbulence and B-Fields simulations in Table 1 and includes both turbulence and magnetic fields, but does not include gravity and therefore does not produce sink particles. We use one snapshot from this simulation only to compare the density PDF with the other four simulations and do not use it in the rest of the analysis.

All of the simulations have a box size of 2 pc and a total mass of M_init = 388 M_⊙, resulting in a mean density of ρ₀ = 3.28 × 10⁻²¹ g cm ⁻³. Periodic boundary conditions are used for all of the simulations. The virial parameter (the ratio of twice the kinetic energy to the gravitational energy) is α_vir = 1.0. The simulations that include magnetic fields have a field strength of B = 10 μG, an Alfvén Mach number of ${{ \mathcal M }}_{{\rm{A}}}=2.0$, and a plasma beta parameter (representing the ratio of thermal and magnetic pressure) of β = 0.33 (Federrath 2015). All of the simulations are isothermal except for the All + Outflows simulation, which is initially isothermal but allows for protostellar heating as implemented in Federrath et al. (2017). No cooling is modeled.

To account for the formation of stars, the simulations use sink particles as described in Federrath et al. (2010a, 2014). Sink particles are formed when a local region undergoes gravitational collapse and exceeds a threshold density, ρ_sink. The Gravity, Turbulence, and B-Fields simulations form sink particles at ${s}_{\mathrm{sink}}=\mathrm{ln}({\rho }_{\mathrm{sink}}/{\rho }_{0})=8.75$. The All + Outflows simulation, which has a higher maximum resolution than the other simulations, forms sink particles at ${s}_{\mathrm{sink}}=\mathrm{ln}({\rho }_{\mathrm{sink}}/{\rho }_{0})=10.14$. Sink particles only form in regions that are maximally refined (ρ_sink is greater than all of the grid level refinement density thresholds). The sink radius is set to 2.5 grid cell lengths (of the maximally refined cells) to avoid artificial fragmentation following the Jeans criterion from Truelove et al. (1998). On all other AMR levels, the refinement criterion is set to refine the local Jeans length by at least 32 grid cells, except for the All + Outflow model, which resolves the Jeans length by at least 16 cells (Federrath et al. 2011). Sink particles continue to accrete gas from any cells within the accretion radius of the sink particle that have exceeded ρ_sink, and are gravitationally bound and converging toward the sink (Federrath et al. 2010a).

The All + Outflow simulation uses a custom module for implementing two-component jet feedback. This two-component jet module was first described in Federrath et al. (2014) and is used in Federrath (2015), as well as for the All + Outflow simulation shown here. The All + Outflow simulation additionally implements protostellar heating as described in Federrath et al. (2017). This two-component jet module consists of a wide-angle low-speed component and a collimated high-speed component, and has been physically calibrated based on theoretical models of jet launching, dedicated high-resolution numerical simulations of single accretion disks with jets (see, e.g., Blandford & Payne 1982; Lynden-Bell 2003), and observations. This jet module uses a velocity profile as shown in Figure 2 of Federrath et al. (2014). For more details about the sink particles used in these simulations and the protostellar outflow prescription, see Table 1 of Federrath et al. (2014). See Federrath et al. (2017) for further details about the protostellar heating feedback.

Projection plots for each run are shown in Figure 1. Examination of these projections indicates that the Gravity and Turbulence runs form large filamentary structures of dense gas with large regions of diffuse gas in between filaments. However, while the density distribution for the Gravity run is smooth, the inclusion of turbulence results in more defined filaments, somewhat more small-scale structure, and greater contrast between high- and low-density regions. In contrast, with the inclusion of magnetic fields and outflows, these large-scale filamentary structures become progressively less prominent, and smaller-scale fluctuations in the density start to dominate.

Furthermore, as demonstrated in Federrath (2015), the differing physics produces vastly different instantaneous SFRs. In addition, the first sink particle is produced at different times in each simulation. The output files are saved based on the integrated star formation efficiency (how much of the gas has been converted into stars); thus, each of the simulation snapshots corresponds to a different physical time. Figure 2 shows the integrated SFE (M_*/M_init) as a function of time for each of the simulations. Note the similarity of the SFE evolution between the Turbulence run and Gravity run is somewhat misleading as these two simulations have very different initial conditions for their density fields. This does not imply that turbulence does not play a role in setting the SFR. In fact, for the same initial condition density field, Federrath (2015) showed that the addition of turbulence reduces the SFR per freefall time by a factor of about two to three. In addition, Figure 2 suggests that outflows play a significant role in slowing down star formation relative to the B-Fields simulation after t/t_ff ∼ 1.3. Protostellar feedback significantly slows star formation only at later times because it takes time for outflows to develop, both because it takes time for the sink particles to form and because time is needed for the outflows to propagate through the gas and affect the gas at larger scales (see, e.g., Gallegos-Garcia et al. 2020).

Zoom In Zoom Out Reset image size

The first snapshot is from just before the first sink particle is formed for each simulation—this is the SFE = 0% snapshot, and we define this to be the t = 0 point in our analysis (see, e.g., Figure 2). Note that the amount of time that passes between when turbulence is fully developed and self-gravity is turned on, and when the first sink particle is formed is different for each simulation. Since the latter point is where we set t = 0, we expect gravity to have already established a power-law tail by t = 0. We choose this as our t = 0 point since we are interested in the evolution of the cloud during star formation and are not currently considering the preceding evolution. Subsequent snapshots may be referred to by their integrated SFE, which can be seen as a proxy for time as a higher SFE refers to a more evolved snapshot. For example, the 4% snapshot is when 4% of the gas mass has been converted into stellar mass (this is the snapshot shown in Figure 1).

For each of the simulations, we consider the SFE = 0%, 1%, 2%, 3%, 4%, 5%, and 10% snapshots. Due to a much longer computational run time, the SFE = 10% snapshot for the All + Outflows simulation was not available, and we instead consider an SFE = 6% snapshot.

3.1. Pure Lognormal Density PDFs

A standard approach for modeling star formation at the cloud scale is to assume a time-independent, lognormal density PDF, which is the expected density distribution for supersonic turbulent molecular clouds (e.g., Vazquez-Semadeni 1994; Vazquez-Semadeni et al. 1995; Padoan et al. 1997; Scalo et al. 1998; Kravtsov 2003; Krumholz & McKee 2005; Hennebelle & Chabrier 2008; Robertson & Kravtsov 2008; Price & Federrath 2010; Burkhart & Lazarian 2012; Collins et al. 2012; Federrath & Klessen 2012; Molina et al. 2012; Hopkins 2013; Walch et al. 2013). In these analytic star formation theories, supersonic turbulence produces gravitationally unstable density fluctuations and sets the overall fraction of dense star-forming gas. The lognormal volume-weighted density PDF is described by

$$\begin{eqnarray}&&{p}_{\mathrm{LN}}(s)=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{s}^{2}}}\exp \left(-\displaystyle \frac{{\left(s-{s}_{0}\right)}^{2}}{2{\sigma }_{s}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

expressed in terms of the logarithmic density,

$$\begin{eqnarray}&&s\equiv \mathrm{ln}(\rho /{\rho }_{0}),\end{eqnarray} \tag{ 2 }$$

where σ_s is the standard deviation of the lognormal distribution. The quantities ρ₀ and s₀ denote the mean density and mean logarithmic density, respectively; the latter of which is related to σ_s by

$$\begin{eqnarray}&&{s}_{0}=-\displaystyle \frac{1}{2}\,{\sigma }_{s}^{2}.\end{eqnarray} \tag{ 3 }$$

For an isothermal equation of state, the width of the lognormal is determined by the turbulent sonic Mach number (Padoan et al. 1997; Federrath et al. 2008; Burkhart et al. 2009; Price et al. 2011; Konstandin et al. 2012) and the turbulence driving parameter b:

$$\begin{eqnarray}&&{\sigma }_{s}^{2}=\mathrm{ln}\left[1+{b}^{2}{{ \mathcal M }}_{s}^{2}\right],\end{eqnarray} \tag{ 4 }$$

where the sonic Mach number depends on the rms velocity dispersion (v_rms) and the sound speed (c_s):

$$\begin{eqnarray}&&{{ \mathcal M }}_{s}=\displaystyle \frac{{v}_{\mathrm{rms},3{\rm{D}}}}{{c}_{{\rm{s}}}}.\end{eqnarray} \tag{ 5 }$$

The turbulent driving parameter, b, describes the mix of solenoidal and compressive modes of the turbulent driving and ranges from b ≈ 1/3 for purely solenoidal driving, to b ≈ 1 for purely compressive driving (Federrath et al. 2008, 2010b). Values of b for real clouds occupy the full range from b ∼ 0.3–1 (see, e.g., Brunt 2010; Price et al. 2011; Ginsburg et al. 2013; Kainulainen et al. 2013; Federrath et al. 2016; Menon et al. 2021), depending on the physical conditions and location of the clouds. For the simulations used in this work, which are driven with a natural mix of forcing modes, resulting in b ∼ 0.4 (Federrath et al. 2010b) and ${{ \mathcal M }}_{s}=5$, Equation (4) predicts a width of σ_s = 1.27.

In the presence of magnetic fields, the width of the lognormal distribution is expected to change. Padoan & Nordlund (2011) and Molina et al. (2012) derive an expression for the width of the lognormal that is modified by the ratio of the thermal to magnetic pressure or the plasma β₀ (Molina et al. 2012):

$$\begin{eqnarray}&&{\sigma }_{m}^{2}=\mathrm{ln}\left[1+{b}^{2}{M}_{s}^{2}\left(\displaystyle \frac{{\beta }_{0}}{1+{\beta }_{0}}\right)\right].\end{eqnarray} \tag{ 6 }$$

For the simulations used in this work, the plasma β₀ ∼ 0.33, resulting in a predicted value of σ_m = 0.83. In case of a strong guide field, Beattie et al. (2021) derived an anisotropic version of Equation (6). In this case, Equation (6) still provides a reasonable approximation, although for ${{ \mathcal M }}_{{\rm{A}}}\lt 2$, anisotropies start to play a significant role (Federrath 2016b).

3.2. A Piecewise Density PDF

Burkhart (2018) proposed using a piecewise lognormal plus power-law density PDF in order to account for the effects of early stellar feedback (e.g., winds, jets, radiation) and self-gravity when calculating the SFE (see also Burkhart & Mocz 2019). In particular, this model accounts for self-gravity by introducing a power-law tail and allows the slope of this power law, as well as the transition density, to vary in response to the influence of, for example, stellar feedback. This change was motivated by observational and numerical work suggesting that the density distribution takes the form of a power law, rather than a lognormal distribution, at the high-density end of the PDF (see, e.g., Klessen 2000; Kritsuk et al. 2011; Collins et al. 2012; Federrath & Klessen 2013; Girichidis et al. 2014; Lombardi et al. 2015; Burkhart et al. 2017; Leroy et al. 2017; Chen et al. 2019; Jaupart & Chabrier 2020; Khullar et al. 2021).

The density PDF proposed by Burkhart & Mocz (2019) consists of a lognormal distribution at low densities and a power-law distribution at high densities, with a transition point denoted as ${s}_{t}=\mathrm{ln}({\rho }_{t}/{\rho }_{0})$, and is described by

$$\begin{eqnarray}{p}_{\mathrm{LN}+\mathrm{PL}}(s)=\left\{\begin{array}{ll}N\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{s}}{e}^{-\tfrac{{\left(s-{s}_{0}\right)}^{2}}{2{\sigma }_{s}^{2}}} & s\lt {s}_{t}\\ {{NCe}}^{-\alpha s} & s\gt {s}_{t},\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where α is the slope of the power-law tail, and σ_s is the width of the lognormal distribution (e.g., given by Equation (4)). The point where the density PDF transitions between the lognormal component and the power-law component can be derived from considerations of continuity and differentiability (Burkhart et al. 2017):

$$\begin{eqnarray}&&{s}_{t}=(\alpha -1/2){\sigma }_{s}^{2}.\end{eqnarray} \tag{ 8 }$$

By combining Equations (8) and (4), the transition point can be further related to the Mach number and driving parameter (Burkhart et al. 2017):

$$\begin{eqnarray}&&{s}_{t}=\displaystyle \frac{1}{2}(2\alpha -1){\rm{ln}}[1+{b}^{2}{{ \mathcal M }}_{s}^{2}].\end{eqnarray} \tag{ 9 }$$

A density PDF with both a lognormal and a power-law component is compatible with a gravo-turbulent model of star formation since the development of a power-law density distribution is consistent with the gas undergoing gravitational contraction while the lognormal distribution still retains some memory of its turbulent initial conditions. In the initial stages of collapse, the power-law slope will be very steep. The slope becomes shallower on a timescale of a freefall time at the critical density (Federrath & Klessen 2013), and may reach −1.5 in the limit of uniform pressure-less spherical collapse (Girichidis et al. 2014; Jaupart & Chabrier 2020; Khullar et al. 2021). Since different values of the power-law slope are expected early in the cloud's evolution, the value of the transition density and fraction of dense self-gravitating gas is also expected to change. In addition, gas may cycle through the power law to lognormal portions of the PDF due to stellar feedback, again adjusting the amount of dense gas available for star formation.

4.1. Overview of the Density PDFs

We now consider the shapes of the simulated density PDFs, which each include different physical processes. Figure 3 shows the volume density PDFs for selected snapshots of each of the simulations described in Section 2. The SFE = 0% snapshot is just before the first sink particle forms in each simulation, the SFE = 2% snapshot is an approximate midpoint in the evolution of each simulation, and the SFE = 5% snapshot is the latest snapshot where the same SFE was available for all of the simulations.

Zoom In Zoom Out Reset image size

All simulations that include turbulence (i.e., Turbulence, B-Fields, and All + Outflows) show a prominent lognormal peak in the PDF, in agreement with expectations for supersonic turbulence (see Section 3.1). For reference, we overplot a single snapshot of the density PDF for the No Gravity run (dashed black line), which includes only turbulence and magnetic fields but no gravity, and a theoretical lognormal shape with a width of σ_s = 1.27 (red dotted line), as predicted by Equation (4), using the driving parameters of the simulations: a sonic Mach number of 5 and a driving parameter b = 0.4. Both lines match the lognormal peak in all three simulations reasonably well. The minor discrepancy with the theoretical curve is likely due to statistical variations of the PDF that occur between individual snapshots (see, e.g., Federrath et al. 2008). We do not make this comparison for the Gravity simulation, as it lacks any turbulent driving, and the narrow lognormal peak of its PDF is simply an imprint of the Gaussian random density field used to initialize this simulation.

At high densities, the simulated PDFs quickly depart from a pure lognormal shape as a power-law tail forms due to the influence of self-gravity, in agreement with the theoretical expectation discussed earlier. The power-law tails have slopes close to −1.9 on average, although there is also some evidence of multiple power-law tails with different slopes. Two or more power-law slopes have previously been observed in molecular clouds, for example, in Schneider et al. (2015a). This may be because of the time-dependence of the emergence of different power-law slopes due to gravitational collapse (Burkhart et al. 2017; Jaupart & Chabrier 2020) and/or due to the formation of rotationally supported structures, i.e., formation of accretion disks (Khullar et al. 2021). The latter may only be seen in the simulation with outflows, as this simulation has sufficient resolution to capture hints of accretion disk formation.

To gauge the transition point between the lognormal peak and power-law tail, we overplot the reference transition density of s_t = 2.25 (vertical black line) calculated from Equations (8) and (4), using a power-law tail slope of −1.9 and the lognormal width of σ_s = 1.27. As the comparison with the No Gravity run shows, this value of s_t is close to the density at which the power-law tail forms and the PDF deviates from a lognormal shape. Note that in the Gravity simulation, which lacks any turbulent driving, we do not expect a lognormal distribution to form, and the power-law tail forms at a significantly lower density. Given that our subsequent analysis mainly focuses on the effects of magnetic fields and protostellar feedback, we choose to use the same s_t in the Gravity simulation as in the runs with turbulence driving. The reference transition density is shown in the Gravity panel of Figure 3 for the sake of comparison with the other panels.

Finally, in addition to the high-density power-law tail, the All + Outflows simulation also deviates from a lognormal distribution at the low-density end of the PDF. When protostellar outflows and heating are included, the low-density end of the PDF exhibits significant fluctuations as the simulation progresses. Figure 4 shows all available snapshots of the All + Outflows simulation, demonstrating that the significant fluctuations in the low-density tail of the PDF are time-varying.

Zoom In Zoom Out Reset image size

Deviations from lognormality at both low and high densities for each of the runs are discussed further below in Section 4.2.

4.2. Fits to the Density PDFs

To compare our density PDFs with the models discussed in Section 3, we fit a lognormal curve to the peak of the density PDF and a power-law relation to the high-density tail of the PDF. We use a separate lognormal fit and power-law fit so that we can focus separately on the behavior of the PDF in the lognormal portion of the PDF and in the power-law portion. Recent work has begun to implement a piecewise lognormal plus power-law fit (Khullar et al. 2021), but for this work, we choose to use separate fits in accordance with the approach used in, e.g., Kritsuk et al. (2011), Schneider et al. (2015a, 2015b), Alves et al. (2017), and Soler (2019).

We fit a Gaussian curve (i.e., a parabola in log space) to the peak of the density PDF using curve_fit from the scipy optimize package, which performs a nonlinear least-squares fit to the data. For each snapshot, we performed a single Gaussian fit to a specific range of densities selected to exclude regions of the PDF that appeared during visual inspection to deviate significantly from a lognormal distribution. The upper limit for all snapshots was set at s = 2.0 in order to exclude regions of the PDF that visually appear to not be purely lognormal. For the Gravity, Turbulence, and B-Fields simulations, no lower limit was set, and all bins below s = 2.0 were used for the fits. For the All + Outflows simulation, we use a lower limit of s = −4.0 in order to exclude the excess of low-density gas in the PDF, produced by protostellar outflows, that deviates significantly from a lognormal PDF (see, for example, Figure 4). These upper and lower limits ensure that each fit minimized the influence of deviations from lognormality at high and low densities. The curve_fit method produces a covariance matrix, from which we derived 1σ uncertainty error bars for each fit.

We also used curve_fit to fit a power-law function at the high-density ends of the PDFs. For each snapshot, we performed a single linear fit ($\mathrm{log}(\mathrm{PDF})=\alpha s+\beta$, in log space, where α is the slope and β is the y-intercept), with a fitting range between s = 3.0–8.0. We selected the upper limit of the fitting range to capture the behavior of the high-density end of the PDF while minimizing the impact of the PDF fluctuations at high densities. For the SFE = 1% snapshot of the B-Fields simulation, we adjusted the fitting range to be between s = 3.0–6.0 because the PDF for this snapshot drops off dramatically above s = 6.0 and skews the fit to a much steeper slope. This is likely due to the formation of a single 3.9 M_⊙ sink particle just before the time of the snapshot. For all of the snapshots, we selected s = 3.0 as the lower limit of the fitting range to exclude portions of the PDFs that are clearly no longer purely linear under visual inspection. We derived the 1σ uncertainties for each fit from the covariance matrix produced by the curve_fit procedure.

The widths and slopes of the fits described above are shown in Figure 5. The top panel of Figure 5 shows the widths of the Gaussian fits for all of the snapshots, with the error bars derived from the covariance matrix, for all runs except the Gravity simulation. The Gravity simulation is not shown because it does not include turbulence, and the narrow lognormal peak apparent in the Gravity panel of Figure 3 is due to the initial Gaussian random field. The bottom panel of Figure 5 shows the fitted slopes for each snapshot of each simulation, with the error bars showing the 1σ uncertainty of the fit. We find that a slope of α = −1.9 (shown with the dotted line) is a reasonable value for the power-law slopes, although we do see values ranging from α = −1.4 to α = −2.4. This range of values agrees with the range of slopes reported for observed GMCs in Burkhart & Mocz (2019). It is also well within the theoretically expected range of slopes (−1 to −3) discussed in Burkhart & Mocz (2019). Thus, we will use the value of α = −1.9 for determining a constant transition density in our later analysis. This value of the slope is shown in Figures 3 and 4 for reference. The value of the transition density calculated using a slope of −1.9 is also shown in Figures 3 and 4 and aligns well with where the PDF deviates from a lognormal distribution.

Zoom In Zoom Out Reset image size

Interestingly, the top panel of Figure 5 shows that the fitted widths fall above the predicted value of σ_s = 1.27 that we find using Equation (4) if we assume that ${{ \mathcal M }}_{s}=5$ and b = 0.4 are fixed (black dashed line in the figure). However, when physics beyond turbulence is included, the values of both ${ \mathcal M }$ _s and b can change (Jaupart & Chabrier 2020; Khullar et al. 2021). For example, gravity will possibly drive more compressive modes of turbulence (Jaupart & Chabrier 2020; Khullar et al. 2021). The formation of the power-law tail at high densities may also bias the fits of the lognormal portion of the density PDF toward a wider value. Likewise, outflows may induce more solenoidal and compressive modes of turbulence (see also Offner & Chaban 2017; Rosen & Krumholz 2020).

We check whether an increase in the sonic Mach number may be the reason for this increase in the fitted PDF widths by plotting the evolution of the sonic Mach numbers versus time in Figure 6. The solid lines show the rms sonic Mach number for each snapshot as calculated using the magnitude of the gas velocity and the sound speed. The dashed black line shows the sonic Mach number that is continuously driven during the simulation run (${{ \mathcal M }}_{s}=5$).

Zoom In Zoom Out Reset image size

The rms sonic Mach numbers for the Turbulence and B-Fields simulations remain fairly close to the input value. In contrast, the rms sonic Mach number for the All + Outflows simulation rises significantly over time, due to the added kinetic energy from the protostellar outflows. In fact, protostellar outflows significantly increase the rms sonic Mach number only after about a freefall time. As discussed with respect to Figure 2, outflows have a more significant impact on the simulation at later times because of the gradual increase of the SFR and the time required for the outflows to propagate through the gas.

Because the rms sonic Mach numbers do not diverge from the driving value of ${{ \mathcal M }}_{s}=5$ with the inclusion of magnetic fields and gravity, the fitted widths in Figure 5 for those simulations cannot be explained solely by a change in sonic Mach number. In the case of the All + Outflows simulation, a sonic Mach number of ≈7 could explain the fitted widths. However, at later snapshots, the rms Mach number increases dramatically, whereas the fitted width does not. Furthermore, at early snapshots, the rms Mach number is too low to reflect the PDF width. These results suggest that either b must be changing or the fitting procedure overestimates the width due to the influence of the power-law tail and/or due to the low-density fluctuations from the outflows.

4.3. Transition Density between the Lognormal and the Power-law Tail

Using our fitted slopes and the width predicted by Equation (4), we calculate the transition density (Equation (8)) for each snapshot. As discussed above, the fitted widths (shown in Figure 5) are wider than the value predicted by Equation (8). This may be due to a bias in the fits from the influence of the power-law tail or due to a changing value of b from additional compressive motions. Furthermore, in Burkhart (2018), the width used to calculate the transition density is determined only by the properties of the large-scale turbulent flow (i.e., the driving-scale sonic Mach number and the turbulent driving parameter b). In the context of these simulations, this suggests the use of the known parameters of the driven turbulence, ${{ \mathcal M }}_{s}=5$ and the driving-scale b = 0.4, to infer the width. Thus, we calculate the transition density using the width predicted by Equation (4), which depends on only the large-scale, driven turbulence, in accordance with Burkhart (2018).

We plot the computed values of s_t versus time in Figure 7. We show a reference transition density of s_t = 2.25 (dashed black line), computed using a reference power-law slope of α = −1.9 and the width predicted by Equation (4). This reference transition density s_t = 2.25 is close to a median value around which the fluctuations occur.

Zoom In Zoom Out Reset image size

There is significant variation of s_t with time, which reflects the variation in the fitted slopes, as shown in the bottom panel of Figure 5. The inferred values of s_t for the Turbulence and All + Outflows simulations generally correspond to the reference value across the snapshots. The B-Fields simulation shows significantly more variation of s_t with time. This is because the power-law slope is significantly more time variable, for the reasons discussed above.

We overplot the reference transition density (s_t = 2.25) as a vertical line in Figures 3 and 4. The s_t = 2.25 transition density provides a compelling point of divergence between the lognormal and power-law portions of the PDF. It closely matches the density at which the simulation without gravity (black dashed lines) and the simulations with gravity and driven turbulence no longer track each other. We therefore use the reference transition density of s_t = 2.25 for all of our subsequent calculations (Section 4.4) in order to provide a single reference density for comparison of all simulations. Future work will further investigate a variable transition density and different methods for fitting the transition density.

4.4. Mass Flow from Diffuse to Collapsing Gas as Traced by the PDF

In this section we analyze the total mass contained in the lognormal and power-law portions of the density PDF and study how that mass moves into the sink particles. We considered the total mass contained in the sink particles, in the density PDF overall, in the power-law tail, and in the lognormal part of the PDF. The power-law tail is defined as any gas that is denser than the reference transition density, s_t = 2.25. The lognormal peak is defined as any gas that is less dense than the reference transition density.

Figure 8 shows the mass contained within each component of the density PDF as a function of time for each simulation. The horizontal red line shows the initial mass of the simulation for reference. The total gas mass (orange line) decreases in time as sink particles form and accrete mass (blue line). Figure 8 also separately shows the mass of the diffuse gas in the lognormal portion of the PDF (i.e., at s < s_t ; purple line) and the dense, self-gravitating mass in the power-law tail (i.e., at s > s_t ; green line). The power-law tail mass is generally stable in time for all of the simulations, while the lognormal peak mass steadily decreases in step with the total gas mass decrease. Overall, the Gravity and Turbulence simulations exchange gas mass into stellar mass on a significantly faster timescale than the other two simulations. Because these simulations evolve much faster, the x-axes for these two simulations are different than the other panels in order to make the evolution easier to see.

Zoom In Zoom Out Reset image size

Figure 8 shows that the diffuse gas reservoir is replenishing the power-law tail mass at a rate that is equal to the SFR, such that it maintains a stable power-law tail. The main difference between the four simulations is that the inclusion of magnetic fields and outflow feedback results in slower accretion of gas from the lognormal portion of the PDF to the power-law tail, in agreement with previous work (see, e.g., Cunningham et al. 2012; Federrath 2015; Offner & Chaban 2017; Rosen & Krumholz 2020).

To illustrate this difference in the accretion rates, Figure 9 compares the time-averaged rates of change for each of the PDF components. To produce this plot, we fit a linear function to each of the trends in Figure 8 using the LinearRegression method from scikit-learn's linear_model, which performs an ordinary least-squares linear regression. Figure 9 shows the slopes, or ΔM/Δt, of each of these linear regression fits. The total gas mass ΔM/Δt values are not shown in Figure 9 because they are identical and opposite to the stellar mass slopes. Note when considering Figure 9 that the linear regressions fit the overall trend and do not capture variations in the slope. For example, a fit to only the final two or three points of the Gravity case would result in a much steeper slope than is measured when fitting all of the points.

Zoom In Zoom Out Reset image size

For all four simulations, the change in the power-law tail mass with time is small, as expected from Figure 8. We also find that the slopes for the lognormal and stellar mass components are similar to each other in all four cases. However, for the B-Fields and All + Outflows simulations, the change in the lognormal mass with time is comparatively much smaller than in the other simulations.

We also consider the cumulative change in mass for the power-law tail mass, the lognormal peak mass, and the stellar mass of each simulation in Figure 10. We plot the mass lost from the power-law tail (green line) and the lognormal peak (purple line), as well as the mass gained in the sinks (blue line). In other words, Figure 10 shows the negative, cumulative change in mass for the power-law tail and the lognormal peak but the positive, cumulative change in the sink mass. This allows us to track which part of the density PDF is ultimately contributing to the growth of the sink particles.

Zoom In Zoom Out Reset image size

We see a similar situation in all four simulations. The lognormal mass loss and the total stellar mass gain track each other closely, while the power-law tail mass loss remains close to zero or slightly negative (indicating a small mass gain). This suggests, in agreement with the discussion of Figures 8 and 9, that the mass in the power-law portion of the PDF is replenished from the lognormal portion at the same rate that it loses mass to the sink particles. However, the timescale on which this occurs is different for each simulation. As seen above, the two simulations with magnetic fields (B-Fields and All + Outflows) evolve much more slowly. Similar to Figure 8, the x-axes are different between the runs with and without magnetic fields.

Furthermore, for the All + Outflows simulation, the cumulative mass change remains shallower at all time steps than in the other simulations, suggesting a slower SFR for this simulation. The lognormal peak mass loss drops slightly below the stellar mass gain in the last few time steps of the All + Outflows simulation, and the power-law tail mass loss becomes slightly positive. This suggests that, in these time steps, protostellar outflows slow down the rate at which the lognormal replenishes the power-law tail, and the power-law tail is slightly depleted. This may be due to the additional kinetic energy injected by protostellar outflows into nearby gas, making it difficult for diffuse gas to collapse.

Interestingly, the Turbulence simulation converts all of the power-law tail mass from t = 0 (M_PL = 11.3 M_☉ at t = 0) into stars within a freefall time of the transition density, i.e., by t = t_ff(ρ_t). This suggests that large-scale driven turbulence does not provide any additional support in the collapsing dense gas. This agrees with the fact that turbulence can lead to increased fragmentation as well as longer overall collapse timescales due to the production of additional diffuse gas (i.e., density fluctuations via the creation of a lognormal distribution). Thus, turbulence also slows the global rate of mass transfer from the lognormal to the power-law tail (Offner et al. 2009; Federrath 2015; Rosen et al. 2019) on timescales longer than t = t_ff(ρ₀).

When magnetic fields are included, all of the mass rates of change are smaller, indicating that the mass transfer from the diffuse, lognormal portion of the PDF to the dense gas and into the stars is stalled and that the depletion times in these simulations are significantly longer than the Gravity and Turbulence runs. Comparison of the timescales of stellar mass growth (Figure 10) indicates that inclusion of magnetic fields results in ∼4–5 times longer depletion times. This is consistent with the effects on the star formation efficiency seen in Figure 2 and discussed in further detail in Federrath (2015). Magnetic fields are able to provide support to the gas against collapse and converging motions at all scales and densities. This is in agreement with the findings of Rosen & Krumholz (2020), who found that magnetic fields strongly suppress fragmentation. In addition, the mass rates of change are even smaller for the All + Outflows case. This effect may be due in part to mass lost to the outflows or may be a result of the increased kinetic energy seen earlier. Our results suggest that both magnetic fields and protostellar outflows are extremely important for inhibiting the accretion of diffuse gas into collapsing gas and of collapsing gas into stars, and that turbulence alone does very little to prevent collapse in the dense gas of the power-law tail.

In this paper, we investigated the density PDFs of a suite of four simulations of a 2 pc region of a star-forming molecular cloud. We find that the inclusion of physical processes beyond turbulence results in deviations from a pure lognormal density PDF. Specifically, self-gravity produces a power-law tail in the high-density, collapsing gas, and protostellar feedback produces other deviations from lognormallity, particularly at the low-density end of the PDF. Thus, the inclusion of self-gravity and protostellar feedback in realistic cloud environments implies non-lognormal density PDFs.

Our results agree with recent observed density PDFs, which have suggested evidence of non-lognormality. For example, some recent observational and numerical work shows that the density PDF of the dense gas in molecular clouds is found to have a power-law distribution (for observational work see, e.g., Collins et al. 2012; Myers 2015, 2017; Schneider et al. 2015b; Alves et al. 2017; Kainulainen & Federrath 2017; Chen et al. 2019; Ma et al. 2021; for numerical and theoretical work, see, e.g., Klessen 2000; Kritsuk et al. 2011; Federrath & Klessen 2013; Girichidis et al. 2014; Lombardi et al. 2015; Burkhart et al. 2017; Chen et al. 2018; Mocz et al. 2017; Jaupart & Chabrier 2020; Khullar et al. 2021). However, with the inclusion of stellar feedback, we no longer observe a lognormal density PDF at the lowest densities and instead observe significant time-dependent fluctuations in the low-density part of the PDF. We only see evidence of a stable lognormal distribution near the mean density. Future work with these simulations will include constructing the column density PDFs and comparing them to observations such as those in Lombardi et al. (2015).

We also find that the measured lognormal widths from these simulations do not match the width predicted by the variance-sonic Mach number relationship (Equation (4)), if we assume that the large-scale turbulence driving parameter b = 0.4 and the sonic Mach number ${{ \mathcal M }}_{s}=5$ are fixed throughout the simulations (e.g., Figure 5, top panel). We find that the rms sonic Mach number cannot account for the wider PDFs, although we do find that outflows inject additional kinetic energy and increase the value of ${ \mathcal M }$ _s (see Figure 6), in agreement with the findings of Rosen et al. (2020). Our PDF fit may be biased to wider values due to the influence of the power-law tail. It is also likely that gravity changes the turbulent driving parameter, b, by increasing the production of compressive modes (e.g., Jaupart & Chabrier 2020; Khullar et al. 2021). However, Pan et al. (2016) found in large-scale (∼200 pc) simulations that b stays relatively low, though their setup differs significantly from ours in that gravitational collapse on small scales is not well resolved; hence, there is only relatively little effect from gravity on the density distributions in their simulations. Further work is needed to understand whether and how the turbulent driving parameter b changes with the inclusion of physics beyond supersonic turbulence. In summary, we find that the inclusion of physics beyond turbulence induces significant changes in the density PDF (the power-law tail and the low-density fluctuations discussed above).

Our results point to the importance of the overall rate at which mass flows from diffuse to dense states in setting the SFR. In these simulations, the amount of dense gas in the power law is similar and remains roughly constant in time (i.e., Figures 8 and 9). Figure 10 shows that mass moves from the lognormal to the power-law tail at the same rate that mass moves into sink particles, resulting in a minimal change in the total mass contained in the power-law tail. This is true for all four simulations studied here, despite the different physics included.

We find that turbulence does not provide significant supportive pressure in the power-law tail. Indeed, the presence of shocks in a supersonic flow will enhance collapse locally around the shock. Furthermore, turbulence decays toward smaller scales, corresponding to less support for higher-density gas in collapsing environments. However, shocks also produce rarefied gas and hence a broad lognormal distribution, which helps to slow collapse in the diffuse gas, as the less-dense gas has a longer freefall time. Thus, over a freefall time, turbulence overall acts to reduce the SFR.

Our results provide evidence for the importance of mass cycling on cloud scales in setting lower SFRs when magnetic fields and stellar feedback are considered. The time-varying fluctuations in the low-density end of the volume density PDF for the All + Outflows simulation provide evidence of mass cycling between high and low densities. The simulated clouds collapse under the influence of gravity, which leads to the local accretion of mass onto the sink particles (e.g., Figures 8–10). However, Figure 4 demonstrates that the amount of low-density gas also changes with time, initially growing rapidly and then falling again. This suggests that gas is cycling between high- and low-density states.

Semenov et al. (2017, 2018) explored an analogous process in galaxy simulations and provided a physical framework for connecting gas depletion times and SFRs with the timescale of this cycling. This framework relies only on mass conservation of gas as the gas cycles between different interstellar medium states, and therefore this framework can also be used to describe depletion times and SFRs of individual star-forming regions like the ones that we study here.

To apply this framework in our simulations, we can consider the power-law tail of the PDF as the actively star-forming gas reservoir that can exchange mass with a more diffuse and non-star-forming lognormal part. The depletion time of the entire region, ${\tau }_{\mathrm{dep}}={M}_{\mathrm{gas}}/{\dot{M}}_{\star }$ can be expressed as τ_dep = τ_dep,PL/f_PL, where ${\tau }_{\mathrm{dep},\mathrm{PL}}={M}_{\mathrm{PL}}/{\dot{M}}_{\star }$ is the depletion time of the power-law tail, and f_PL = M_PL/M_gas is its mass fraction. The latter can be expressed via the timescales of dense gas supply from the lognormal part, ${t}_{\mathrm{LN}}$, and the lifetime of dense gas before it is converted into stars or removed back to the diffuse state by, e.g., turbulent shear or stellar feedback, t_PL: ${f}_{\mathrm{PL}}\sim {t}_{\mathrm{PL}}/({t}_{\mathrm{LN}}+{t}_{\mathrm{PL}})$.

Figures 8–10 show that the power-law mass fraction f_PL is nearly constant in time and does not change with the inclusion of magnetic fields and protostellar feedback. This indicates that the timescales ${t}_{\mathrm{LN}}$ and t_PL are set predominantly by turbulence and gravity. At the same time, the timescale on which the power-law tail is converted into stars, τ_dep,PL, increases when magnetic fields and protostellar outflows are included.

Burkert (2017) put forward an analogous analytic model, also based on mass conservation, that assumes that the evolutionary timescales of the diffuse and dense molecular states are set by the freefall time at the characteristic densities of these states, ⁸ t_ff,diff and t_ff,dense (see Equation (1) and the definition of the accretion rate in Burkert 2017). After substituting these timescales in the above expression for the dense mass fraction and assuming t_ff,diff ≫ t_ff,dense, we get f_dense ∼ t_ff,dense/(t_ff,diff + t_ff,dense) ∼ t_ff,dense/t_ff,diff, which is also the result obtained by the Burkert (2017) model. Thus, this model produces a constant f_dense if the average densities of dense and diffuse gas reservoirs are constant and the timescales of gas exchange between these reservoirs are set by corresponding freefall times. However, in our preliminary study, we find that, while the gas in the power-law tail indeed evolves on its local freefall time, the more diffuse lognormal component evolves on shorter-than-freefall timescales, indicating the importance of turbulence in controlling this evolution.

Our future work will focus on a more detailed analysis of the evolutionary timescales and their dependence on the parameters of turbulence and other physical processes included in our simulations.

In this paper, we analyze a suite of four hydrodynamical simulations in order to gain a deeper understanding of how individual components of the density PDF are affected by turbulence, magnetic fields, gravity, and protostellar outflow feedback. We find the following:

• 1.  
  The inclusion of protostellar outflows produces time-varying, non-lognormal features in the low-density end of the density PDF, which results in an excess of diffuse gas relative to simulations that do not include feedback.
• 2.  
  A power-law tail is evident at the high-density end of the density PDF for all simulations that include self-gravity. The power-law tail mass and slope do not significantly vary in time despite the presence of outflow feedback.
• 3.  
  The density distribution retains a lognormal shape only around the mean density, reflecting the presence of large-scale turbulence. However, we measure wider lognormal fits than would be predicted by the variance-sonic Mach number relationship using the driven Mach number and the driving-scale forcing parameter b. This may be because the individual lognormal fits are biased by the growth of the power-law tail and/or because gravity modifies the value of b.
• 4.  
  Mass flows from the lognormal portion of the density PDF to the power-law tail at a rate that is nearly equal to the sink particle mass accretion. This implies that stars form at the rate at which the diffuse gas contracts and replenishes the power-law tail of the density PDF.
• 5.  
  Unlike magnetic fields, we find that driven turbulence does not provide significant support in the dense gas of the power-law portion of the PDF. Over a freefall time at the transition density, our supersonic hydrodynamic turbulence simulation converts all of the mass in the power-law tail into stars.
• 6.  
  The inclusion of magnetic fields and protostellar outflows slows the transference of mass from gas into stars, as well as from the diffuse gas to the dense collapsing gas, by approximately a factor of ∼4–5.

Our results shed light on the way that the gas mass flows between different parts of the density PDF and highlights the importance of such flows in setting slow SFRs in molecular clouds. Low SFRs result from the slowdown of the net mass transfer due to the combined effects of turbulence, gravity, magnetic fields, and protostellar feedback. Our future work will take a closer look at how mass flow occurs, how it connects to analytic star formation models that use the density PDF, and how this process depends on a wider parameter spread of magnetic field strength and turbulent driving.

S.A. thanks Matthew E. Orr and Charlotte Olsen for their valuable comments on the manuscript. S.A. is grateful for the support of NASA grant No. 80NSSC20K0500. The analysis and one of the simulations for this work was performed using computing resources provided by the Flatiron Institute. B.B. is grateful for funded support by the Simons Foundation, Sloan Foundation, and the Packard Foundation. Support for V.S. was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51445.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. C.F. acknowledges funding provided by the Australian Research Council (Future Fellowship FT180100495), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). 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). A.L.R. acknowledges support from NASA through Einstein Postdoctoral Fellowship grant No. PF7-180166 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060; and support from Harvard University through the ITC Post-doctoral Fellowship. The software used in this work was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago. The authors acknowledge the use of the following software: yt (Turk et al. 2011), flash (Fryxell et al. 2000), SciPy (Virtanen et al. 2020), scikit-learn (Pedregosa et al. 2012), Matplotlib (Hunter 2007), and astropy (Astropy Collaboration et al. 2013).