- 3

A Markov model for non-lognormal density distributions in compressive isothermal turbulence

Philip Mocz ∗ 1 and Blakesley Burkhart 2,3

- 1 Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ, 08544, USA ∗ 2 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Rd, Piscataway, NJ 08854, USA

(Received xxx xx, 2019; Revised xxx xx, 2019; Accepted xxx xx, 2019)

Submitted to ApJL

## ABSTRACT

Compressive isothermal turbulence is known to have a near lognormal density probability distribution function (PDF) with a width that scales with the sonic Mach number and nature of the turbulent driving (solenoidal vs compressive). However, the physical processes that mold the extreme high and low density structures in a turbulent medium can be different, with the densest structures being composed of strong shocks that evolve on shorter timescales than the low density fluid. The density PDF in a turbulent medium exhibits deviations from lognormal due to shocks, that increases with the sonic Mach number, which is often ignored in analytic models for turbulence and star formation. We develop a simple model for turbulence by treating it as a continuous Markov process, which explains both the density PDF and the transient timescales of structures as a function of density, using a framework developed in Scannapieco &amp; Safarzadeh (2018). Our analytic model depends on only a single parameter, the effective compressive sonic Mach number, and successfully describes the non-lognormal behavior seen in both 1D and 3D simulations of supersonic and subsonic compressive isothermal turbulence. The model quantifies the non-lognormal distribution of density structures in turbulent environments, and has application to star forming molecular clouds and star formation efficiencies.

Keywords: methods: statistical - stars: formation - turbulence

## 1. INTRODUCTION

The interstellar medium (ISM) is a supersonic, highly turbulent, magnetized medium (Larson 1981; Shu et al. 1987). The turbulent density distribution plays a crucial role in the process of star formation, including setting the initial mass function, star formation rates and star formation efficiency. There is a well established numerical and analytic result that the density distribution of a supersonic, isothermal medium is approximately lognormal, with a variance that grows with the sonic Mach number M (Passot &amp; V´zquez-Semadeni 1998; Padoan a &amp; Nordlund 2011a; Price et al. 2011; Molina et al. 2012; Burkhart &amp; Lazarian 2012). Past a critical density threshold, the high density tail of the distribution may experience collapse under self-gravity, and form stars. It is thus important to understand the distribution and timescales of high density structures from a theoretical perspective, to understand the initial conditions of star formation and predict star formation efficiencies.

Corresponding author: Philip Mocz pmocz@astro.princeton.edu

- ∗ Einstein Fellow

Despite the wide usage of the lognormal model (Vazquez-Semadeni 1994; Federrath et al. 2008; Price et al. 2011; Hopkins 2013), supersonic turbulence is known to also exhibit deviations from a lognormal density PDF deviations from a lognormal density PDF due to the presence of strong shocks and intermittency (Kowal et al. 2007; Federrath 2013; Squire &amp; Hopkins 2017; Pan et al. 2019). The deviation becomes stronger at higher Mach numbers (i.e., for stronger density compressions and larger degree of density contrast). The non-lognormal density PDF can be described well with a simple spatial model that considers the density to be

arranged as a collection of strong shocks of width M -2 (Squire &amp; Hopkins 2017).

Recently, Robertson &amp; Goldreich (2018) put forth the theoretical picture (confirmed by simulations) that a supersonic, turbulent, isothermal medium is a network of quasi-1D strong shock sheet fragments with exponential post-shock profiles 1 . These strong shocks make up most of the mass fraction (and high density tail of the density PDF), but fill just a tiny fraction of the volume. The high density shocks are transient structures, that evolve and dissipate on faster timescales than the cloud crossing time (but see also Falceta-Gon¸alves et al. 2008). c

A complimentary picture to describing the density PDF as a spatial composition as was done in Robertson &amp; Goldreich (2018), is to consider the statistical timeevolution of the densities of fluid elements. The volume elements that make up the density PDF, for example may be considered a Markov process (similarly, the velocity field in a turbulent Navier-Stokes velocity field can be described by a Markov process Novikov 1989; Pedrizzetti &amp; Novikov 1994; Pedrizzetti 1999; Renner et al. 2001). Indeed, a lognormal distribution is suggestive of an Ornstein-Uhlenbeck stochastic process, i.e., a mean-reverting random walk (Uhlenbeck &amp; Ornstein 1930). The density PDF in supersonic turbulence has been modelled as such in Konstandin et al. (2016). More recently, Scannapieco &amp; Safarzadeh (2018) considered a Markov model for explaining not only the turbulent PDF, but also the rate of change of the log density ( s ) over a finite time interval, given as a function of final value of s . Here we adopt a similar approach to study the non-lognormality of the density PDF. The important differences between our model and theirs are: (1) we seek to describe the non-lognormality of the density PDF whereas Scannapieco &amp; Safarzadeh (2018) was designed to predict near lognormal distribution, (2) our model is fully analytic, (3) we describe the behavior in the Eulerian frame as opposed to Lagrangian, and (4) we ensure exact mass conservation with proper renormalization.

Thus our focus is to construct a simple, predictive Markov model to describe the non-lognormal distribution of density in compressive turbulence. We do so by applying simple assumptions for the timescales of evolving structures as a function of density. We then measure the timescales in 1D simulations of isothermal 'turbulence' - such simulations are entirely shock-dominated, as there is no true turbulence in 1D. Thus, the simple set up we work with here is akin to Burgers turbulence.

1 see also Mocz &amp; Burkhart 2018 for the magnetohydrodynamic version of the picture

These 1D simulations allow us to fit the free parameters of our model and study the dependence on the sonic Mach number. We then compare our model to 3D simulations, which exhibit an interplay between strong shocks and a turbulent cascade.

We lay out our notation and a brief overview of the isothermal fluid equations and the lognormal model in § 2 and 3. We describe out analytic Markov model for transient shocked-structures turbulence in § 4, as well as fitting the model parameters with 1D simulations. We compare our model to previous models in the literature for the PDF, and to results from 3D supersonic isothermal turbulence simulations in § 5. We summarize our main findings in § 6.

## 2. ISOTHERMAL FLUID EQUATIONS

The equations for an isothermal, inviscid fluid can be written compactly as follows. Let ρ be the average density in the medium, and use units where the sound-speed is c = 1. The equation of state of such a fluid is P = ρc 2 = ρ . We define the log-density as: s ≡ log( ρ/ρ ). Then, the isothermal fluid equations are:

$$D _ { t } s = - \nabla \cdot \mathbf v$$

$$D _ { t } \mathbf v = - \nabla s$$

where D t = ∂ t + v · ∇ is the convective derivative. In the Lagrangian frame, density changes are entirely due to the divergence of the velocity field (compressibility of the fluid).

In the Eulerian frame, the corresponding equations are:

$$\partial _ { t } s = - \mathbf v \cdot \nabla s - \nabla \cdot \mathbf v \quad \quad \quad ( 3 )$$

$$\partial _ { t } \mathbf v = - \mathbf v \cdot \nabla \mathbf v - \nabla s$$

Changes in density at a fixed location are due to a combination of advection and compression.

## 3. LOG-DENSITY PDF

Our object of study is the density PDF. Consider the volume-weighted PDF of the log-density, P V ( s ), which is the probability of a given unit volume V having an average density ρ = exp( ), per unit s s .

The PDF must obey 2 constraints, the conservation of probability,

$$\int _ { - \infty } ^ { \infty } P _ { V } ( s ) \, d s = 1,$$

and the conservation of mass,

$$\int _ { - \infty } ^ { \infty } \rho P _ { V } ( s ) \, d s = \overline { \rho } = 1.$$

We can also define the related mass-weighted density PDF, as P M ( s ) = ρP V ( s ).

A common model for the density PDF is the lognormal:

$$\omega \colon \\ P _ { V } ^ { \text{LN} } ( s ) = \frac { 1 } { \sqrt { 2 \pi \sigma ^ { 2 } } } \exp \left [ - \frac { ( s - s _ { 0 } ) ^ { 2 } } { 2 \sigma ^ { 2 } } \right ] \quad \ ( 7 )$$

where s 0 = -σ / 2 2 is the center of the normal distribution, and σ 2 is its variance. The variance has been empirically determined from numerical simulations of driven isothermal turbulence to be a function of the sonic Mach number and the turbulent driving (b):

$$\sigma ^ { 2 } = \log ( 1 + b ^ { 2 } \mathcal { M } ^ { 2 } )$$

where b glyph[similarequal] 1 for compressive driving and b glyph[similarequal] 1 / 3 for solenoidal (Federrath et al. 2008; Price et al. 2011). If the volume-averaged PDF is lognormal, the massweighted PDF is too, just with shifted mean: s 0 → σ / 2 2. A well-known motivation for the origin of a lognormal PDF comes from considering the statistics of the divergences in the velocity field in the right-hand-side of Eqn. 1 responsible for evolving s and applying the Central Limit Theorem (e.g. Vazquez-Semadeni 1994). If a Lagrangian fluid element of log-density s encounters instances of diverging velocity fields in the fluid which behave like independent random variables, then by the Central Limit Theorem, the mass-weighted distribution of the log-density s is normal, and thus so is the volumeweighted distribution. Any deviation from the assumption of independent random variables will lead to deviations in the lognormal distribution.

## 4. A MARKOV MODEL FOR NON-LOGNORMAL DENSITY DISTRIBUTION IN COMPRESSIBLE TURBULENCE

We model the volume-weighted density PDF in the Eulerian frame as a Markov process:

$$s ( t + d t ) = s ( t ) + A ( s ) d t + \mathcal { N } ( 0, 1 ) \sqrt { D ( s ) d t } \quad ( 9 ) \quad \text{dic}$$

where

$$A ( s ) = \frac { s _ { * } - s } { \tau _ { A } ( s ) } \quad \quad ( 1 0 )$$

is a drift term that returns s towards the mean s ∗ of the distribution on a timescale of τ A ( s ), and

$$D ( s ) = \frac { 2 \sigma _ { s } ^ { 2 } } { \tau _ { D } ( s ) } \quad \quad ( 1 1 )$$

is is a diffusion term, with timescale τ D ( s ) and N (0 , 1) is a random number drawn from a normal distribution of mean zero and variance one.

In the case that τ A ( s ) = τ D ( s ) = const, the Markov process is the well-known Ornstein-Uhlenbeck process

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00544v2-figure-ddda493a-ab21-434b-b7da-3c2053c3ead2.jpg)

4

0.1

τ

3

2

1

0

1

2

3

-

-

-

-

4

-

2

0

2

4

Figure 1. Top panel : A plot of our model for the density PDF (Eqn. 15). f = 0 (solid) is a simple Gaussian, while f &gt; 0 introduces negative skewness (shown in dotted lines is f = 0 15). . Bottom panel : Measure of the density-dependent transient timescales in 1D isothermal shock simulations used to tune our model. ∆ s reverts back to the straight gray line after ∼ τ eddy , meaning that the structure at a given location is now uncorrelated with what used be there. Shown is our model fit (very thick, transparent line) and simulation data ( M = 2 (thin), 4, 8, 16 (thick)) at 4 different times. High density structures tend to revert back to mean distribution (gray line) faster than structure at low densities.

(mean-reverting random walk), and results in a lognormal distribution for s (Eqn. 7). Having constant drift and diffusion timescales is the unique solution that predicts a lognormal distribution, although modifications to the timescales can yield distributions that are very close to lognormal (Scannapieco &amp; Safarzadeh 2018).

Despite the turbulent density PDF being near lognormal, the types of structures at high and low densities are quite different. High density gas is composed of strong shocks arranged in quasi-1D sheet fragments that can interact. They have a large mass fraction but small volume-filling fraction. Low density regions are created by the interaction of 3-D expansion waves, structured as large volume-filling irregular voids. The high density part of the PDF establishes itself first very quickly when the medium is turbulently stirred, in just a fraction of the dynamical time (Nordlund &amp; Padoan 1999). These shocks are known to operate on faster timescales than

0.2

eddy eddy

τ

τ

0.8

τ

0.4

eddy eddy

the rest of the fluid (Scannapieco &amp; Safarzadeh 2018; Robertson &amp; Goldreich 2018).

Thus we consider a drift timescale that is constant at low densities s &lt; s ∗ , but decreases at high densities s &gt; s ∗ :

$$\tau _ { A } ( s ) = \frac { \tau _ { A, 0 } } { 1 + \mathrm H ( s - s _ { * } ) 3 f / 2 } \quad \quad ( 1 2 )$$

where H( ) is the Heaviside step function, s τ A, 0 is a characteristic timescale, and f ≥ 0 is a free-parameter setting how strongly the timescale decreases at high densities. This is a phenomenological model meant to capture the fact that strong shocks can form, expand, and revert back to the mean on faster timescales than the typical dynamical time of the system, but it is rather quite general to describe a two-state behaviour which ends up yielding tractable analytic solutions. We also consider a constant dispersion term:

$$D ( s ) = \frac { 2 \sigma _ { s } ^ { 2 } } { \tau _ { D, 0 } } \quad \quad \quad ( 1 3 )$$

where σ s is the size of the dispersion and τ D, 0 is the associated dispersion timescale.

This Markov process has a steady-state analytic solution for the density PDF, which can be obtained from solving the time-steady Fokker-Planck equation:

$$0 = - \frac { \partial } { \partial s } \left [ A ( s ) P ( s ) \right ] + \frac { \partial ^ { 2 } } { \partial s ^ { 2 } } \left [ \frac { D ( s ) } { 2 } P ( s ) \right ] \quad ( 1 4 ) \quad \text{f} \quad \text{i}$$

The solution is:

$$P _ { V } ( s ) \in \exp \left [ - \frac { ( s - s _ { * } ) ^ { 2 } \, ( 1 + f ( s - s _ { * } ) \mathrm H ( s - s _ { * } ) ) } { 2 \tilde { S } }. \right ] _ { ( 1 5 ) } \quad \begin{array} { c } \sigma _ { y s } \end{array} \quad \text{$\sigma _ { y s } \cdot}$$

where we have defined

$$\begin{smallmatrix} \cdot \cdot \cdot \\ \tilde { S } \equiv \sigma _ { s } ^ { 2 } \frac { \tau _ { A, 0 } } { \tau _ { D, 0 } } \end{smallmatrix}$$

The solution depends only on the drift-to-dispersion timescale ratio, τ A, 0 /τ D, 0 , not their absolute values.

The normalization of Eqn. 15 is computed from the probability constraint (Eqn. 5), and s ∗ from the mass constraint (Eqn. 6). The analytic expressions are a bit cumbersome but it is straightforward to compute these values numerically. Note that the lognormal distribution (Eqn. 7) is recovered when f = 0 and τ A, 0 = τ D, 0 .

The model, at this stage, has two free parameters that need to be fit to data: f , a knob for the degree of nonlognormality, and ˜ which sets the width of the distriS bution. In the top panel of Fig. 1 we plot the shapes of the predicted PDFs for f = 0 15 and compare it to . f = 0 (lognormal), for various widths ˜ . S The f &gt; 0 parameter introduces a degree of negative skewness. The deviation from lognormal grows stronger as the PDF is made wider. We will fit the free parameters to 1D 'turbulence' simulations, described in sections 4.1 and 4.1.1.

Figure 2. Space-time diagram for 1D isothermal supersonic M = 4 'turbulence' ( top : density, bottom : velocity). The system quickly establishes a network of self-similar interacting and expanding shocks.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00544v2-figure-6dc577e9-4ecf-4e88-9e76-629fd6946019.jpg)

## 4.1. 1D isothermal 'turbulence' simulations

We carry out 1D isothermal 'turbulence' simulations to which we will tune the Markov model. This section is dedicated to the numerical implementation details. We solve Eqns. 1 and 2 using a 2nd order finite volume (Godunov) method, following Springel (2010) on a Cartesian grid. We add a random driving term to the right-hand-side of Eqn. 2 to drive the turbulence.

Shocks are driven on a scale of L stir = 1. The box size is much larger, L box = 1000, to capture the statistics of shocks. The resolution is ∆ x = 0 0625 to resolve shocks. .

The shocks are driven by a random driving term, which has to be compressive in 1D. The details of the driving do not matter significantly in 1D for the statistical properties of the shocks, unlike in 3D where they can be a mix of compressive and solenoidal. We drive the turbulence by adding 5 × L box Gaussian pulses and uniform random locations x rand L , box in the box, in time

intervals of ∆ t kick = 0 1: .

$$v _ { \text{source} } & = \sum _ { \text{shock=1} } ^ { 5 \times L _ { \text{box} } } \mathcal { A } \times \mathcal { N } ( 0, 1 ) \exp \left ( - x _ { \text{rand}, L _ { \text{box} } } / \left ( 2 L _ { \text{stir} } ^ { 2 } \right ) \right ) \\ \text{throws} \quad \mathcal { A } \colon = \mathcal { H } \, \sum _ { \text{threshold} } \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \, \mathcal { H } \colon }$$

where A is the amplitude of the kicks, which sets the resulting turbulent Mach number. We simulate turbulence with mach numbers M = 2 4 8 16. , , ,

We find the driving establishes a quasi-steady state solution (approximately constant Mach number) where the energy injected is dissipated by the shocks. The result is 'burgulence': a network of 1D shocks (saw-tooth shaped as a function of s ), with a velocity power spectrum of k -2 as in Burgers' turbulence, as opposed to the k -5 / 3 Kolmogorov result from 3D incompressible turbulence. Fig. 2 shows a space-time diagram of the full solution from t = 0 to t = 100.

## 4.1.1. Fitting the Markov model timescales from 1D simulations

Our Markov model is best fit with f = 0 15 and . ˜ = S 1 4 . M with the 1D simulations (bottom panel Fig. 1).

We have carried out a suite of 1D driven isothermal 'turbulence' simulations, described in § 4.1, with sonic Mach numbers of M = 2 4 8 16. The simulations have , , , a characteristic timescale τ eddy = L stir M c associated with the stirring length scale.

We use the simulations to fit τ A, 0 and f , and found τ A, 0 = τ eddy and f = 0 15. The fit was performed as fol-. lows, using an approach similar to that in Scannapieco &amp; Safarzadeh (2018) (that work instead treated Lagrangian tracer particles rather than volume elements). We considered the average change ∆ s of an Eulerian volume element over a time interval, given a final value of s 2 . The average change of our Markov process should evolve as

$$\Delta s ( s _ { 2 }, \Delta t ) = ( 1 - \exp \left [ - \Delta t / \tau _ { A } ( s _ { 2 } ) \right ] ) \left ( s _ { 2 } - s _ { * } \right ) \ \ ( 1 8 ) \quad \text{rep}$$

The relation approaches s 2 -s ∗ at long times; that is, after an eddy time the fluid of a given density s 2 is expected to lose memory of its past and on average revert back to the mean log-density. We measured ∆ s as a function of time and density in our 1D simulations of varying Mach number, and fit it with Eqn. 18, see Fig. 1 (bottom panel). This idea of measuring ∆ s comes from Scannapieco &amp; Safarzadeh (2018), where in that work it was measured in a Lagrangian frame.

Then, we fit the term controlling the diffusion, ˜ , S to match the resulting PDFs of the simulations. We find ˜ = 1 4 S . M for our 1D simulations.

In 1D, our model reduces to a single parameter model (the sonic Mach number M ). We shall see in § 5.1 that in 3D our model still works remarkably well, if M is replaced with the effective compressive Mach number.

## 5. RESULTS

## 5.1. Applying the Markov model to turbulent simulations

We investigate how well our Markov model can predict the non-lognormal density PDF of 1D and 3D turbulent simulations. In Fig. 3, we verify that our model, calibrated to the 1D simulations, also predicts the resulting density PDF of the 3D simulations accurately.

We plot the model and 1D simulation data for Mach numbers M = 2 4 8 16. As the Mach number increases, , , , the density distribution becomes wider, has a more negative mean, and a more negative skewness (due to strong shocks). In addition, we have carried out 1D simulations with Mach numbers M = 0 5 1 to show that our model . , is also applicable to the compressive subsonic and transonic case, where our model reduces to a near lognormal.

We also show the PDFs of very high resolution (4096 3 ) 3D isothermal turbulent simulations with M∼ 17 from Federrath (2013). The tails of the PDF can be tricky to measure due poor statistics at lower resolution, which is why we rely on this high resolution run here for our comparison . 2 There are two different turbulent driving modes simulated: fully compressive, and fully solenoidal. The compressive driving exhibits a much wider and more non-lognormal density distribution than does the solenoidal driving. In the compressive case, the PDF looks similar to a 1D simulation with the a similar Mach number ( M ∼ 16). In the solenoidal case, the PDF looks like that of a 1D simulation with the Mach number reduced by a factor of 4. These are rather remarkable results. That is, the 3D simulations can be modeled with our Markov model (with f = 0 15 and . ˜ = 1 4 S . M as in the 1D case) if the Mach number M is replaced with an effective compressive Mach number describing the system, which we point out is a function of the sonic Mach number and the details of the turbulence driving. In our 1D modeling, there are only compressive (shock) motions. In 3D, even if the large-scale driving is compressive, shock collisions will generate vorticity. If the driving is solenoidal, the vorticity in the fluid is larger (by a factor of ∼ 2; Federrath 2013), because solenoidal motions are directly injected. But in the turbulent cascade, some of those motions are converted into compressive motions and shocks form. The amount of compressive motion will be smaller in solenoidally driven

2 see Federrath (2013) for a discussion of how finite resolution biases the PDF

Figure 3. Volume-average density PDFs of our Markov model and simulations. We compare our Markov model prediction for the PDF (gray lines) to 1D 'turbulence' shock simulations (purple dots) with Mach numbers M = 0 5 1 2 4 8 16. . , , , , , The distributions are close to lognormal, but strong shocks lead to negative skewness that increases with Mach number. Our models also do a good job at explaining the PDFs of 3D turbulent simulations. Shown here are two M∼ 17 high-resolution (4096 3 ) simulations from Federrath (2013), one driven solenoidally, the other compressively. In the compressive case, the PDF looks similar to a 1D simulation with a similar Mach number. In the solenoidal case, the PDF looks like that of a 1D simulation with the Mach number reduced by a factor of 4.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00544v2-figure-7e7c590d-29ca-448b-abeb-9e3be63a5f96.jpg)

turbulence, and thus the strength of shocks will be reduced.

While we have not compared the timescale ∆ s (which is a function of τ A ) in the 3D simulations of Federrath (2013) directly to our model, we point out that we do not expect it to be different from the 1D simulations, up to constant factor, otherwise the resulting density PDF would also change significantly unless τ D somehow conspires to give the same shape.

## 5.2. Non-lognormal PDF statistics

Higher-order moments of the turbulent density distribution can reveal its underlying physical conditions and the nature of turbulent driving (Burkhart et al. 2009). We calculate the resulting mean, variance, and skewness of our model as a function of 1D (compressive) Mach number, and compare it to the lognormal model, shown in Fig. 4. Strong shocks cause non-lognormal density distributions, and evolve on faster timescales than the rest of the fluid, leading to a negative skewness that grows with Mach number.

The variance of the PDF as a function of Mach number is a highly discussed quantity in the study of turbulent statistics (see, e.g., the discussion around Fig. 6 of Hopkins 2013). Across various turbulent simulations in the literature, it has been seen that the skewness as a function of Mach number does not vary simply as log(1 + M 2 compressive ), and it can in fact be larger, and there is large variance at high Mach numbers. The variance is a function of the degree of non-lognormality in the simulation, which can be affected by the details of the turbulent driving. The measured variance in simulations may also be affected due to limitations of resolution and finite box size (which bias the variance to slightly smaller values). Interestingly, our Markov model predicts that when the turbulence is significantly non-lognormal, the variance scales more strongly than log(1 + M 2 compressive ) at high Mach numbers. In fact, it scales linearly with the Mach number (Fig. 4). This linear scaling for variance in 1D simulations is a theoretical upper limit for 3D simulations, where rotational motions in turbulence may reduce the variance.

Figure 4. Mean, variance, and skewness of the density PDF as a function of 1D (compressive) Mach number (dotted), and comparison to the lognormal model (thin). With increased Mach number, the PDF becomes more negatively skewed as strong shocks make the turbulence more nonlognormal, and variance varies linearly with Mach number.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00544v2-figure-10317b41-be5d-4569-a3bb-a6b0ffc0968e.jpg)

The mean of the PDF is simply set by mass conservation. Our non-lognormal Markov model has decreased mean log-density relative to the lognormal model.

## 5.3. Comparison to other PDF models

Scannapieco &amp; Safarzadeh (2018) have also recently considered a Markov model for explaining the PDF in turbulence, which serves as a general framework for our approach. There are a few important differences between our model and theirs. We note that both models choose time scales that are phenomenological, based on simulation data. First, our model has an analytic solution, and describes the strong deviations from lognormal of the turbulence, while Scannapieco &amp; Safarzadeh (2018) choose τ D ( s ) = τ A ( s ) purposefully to yield a near lognormal PDF, and their model has to be evaluated numerically through a Markov simulation. Second, we modeled the volume-averaged PDF with the timescales being measured from simulations in an Eulerian frame, while Scannapieco &amp; Safarzadeh (2018) modeled the mass-averaged PDF with timescales measured for Lagrangian volume elements, although the the PDFs can easily be related to each other. We note however that our timescales cannot directly be compared to Scannapieco &amp; Safarzadeh (2018) because they are measured in different frames and Eulerian quantities include the contribution of advection in changes in s , and hence will always measure larger changes over a fixed time interval. Third, we ensure exact mass conservation through proper renormalization of the PDF, whereas it is only approximate in Scannapieco &amp; Safarzadeh (2018).

Hopkins (2013); Squire &amp; Hopkins (2017) provide an analytic model for non-lognormal PDFs, motivated by a phenomenological scaling of structure functions in intermittent turbulence, characterized by a parameter T , where T = 0 is lognormal, and T = ∞ is maximally intermittent. It was found that the parameter T increases systematically with the compressive Mach number of turbulent simulations. As such, the PDF is described by a single parameter (the compressive Mach number), just as we have found. The model is complimentary to ours, where we have obtained the PDF by considering the timescales of transient structures, rather than their spatial distributions. The T -models do have a time interpretation as well, discussed in Hopkins (2013): they represent the most general steady-state result of a class of multiplicative-random-relaxation processes (with e.g. a Poisson-distributed 'waiting time' distribution between discrete events that produce random multiplicative changes to the local value then damp exponentially). It is related to our model, except for the assumptions about the discrete nature of 'events' and damping. The T -model also has the nice mathematical feature that when the density field is convolved on different scales with a filter, the field statistics of log-density remain exactly a T -model, as opposed to log-normal distributions where convolution can lead to small biases in the statistics.

## 5.4. Applications of the PDF to star formation

In the theoretical picture of a turbulent ISM, gas becomes gravitationally unstable beyond a critical density. Such a threshold has been investigated analytically and with simulations (Krumholz &amp; McKee 2005; Hennebelle &amp; Chabrier 2011; Padoan &amp; Nordlund 2011b; Federrath &amp; Klessen 2012; Mocz et al. 2017; Burkhart 2018; Burkhart &amp; Mocz 2019). The star formation efficiency per freefall time is the fraction of gas above this threshold.

For illustrative purposes, we take a transition density 3 of s crit = log( M 2 ). This is the post-shock density in the isothermal limit. The post-shock density can be taken as a critical density for collapse in virialized cores without strong magnetic fields (Krumholz &amp; McKee 2005; Padoan &amp; Nordlund 2011b; Burkhart &amp; Mocz 2019). The star formation efficiency (mass-weighted) per unit freefall time is estimated as:

$$\phi = \int _ { s _ { \text{crit} } } ^ { \infty } P _ { M } \, d s \, \text{ \quad \ \ } \quad \ \ ( 1 9 )$$

In Eqn. 19 there is no fudge-factor to account for the removal of gas by stellar winds or feedback. We will examine here how non-lognormality affects star formation efficiency, which are observed to be a few percent (Krumholz et al. 2012).

Our definition of the SFE in Eqn. 19 is similar to that of Burkhart &amp; Mocz (2019), where the efficiency is estimated as the fraction of gas above the critical density, relative to the total gas. The primary difference between our approach and that of Burkhart &amp; Mocz (2019) is that the form of the density PDF is different. Burkhart &amp; Mocz (2019) consider Eqn. 19 using a power law form of the density PDF, which is the expected form considering gas that is influence primarily by its own self-gravity. In our model, there is no self-gravity and we only take into account the non-lognormality of the PDF. We also do not include the 'multi-free fall' density inside the integral as is done in studies such as (Federrath &amp; Klessen 2012), which give a very strong increase of the star formation efficiency per free fall with sonic Mach number. The 'multi-free fall' scenario re-weights the integrand of Eqn. 19 by t ff ( ρ /t ) ff ( ρ ).

Fig. 5 shows the star formation efficiency of our Markov model compared to the prediction for a lognormal, as a function of sonic Mach number. The star formation efficiency is increased in our model for turbulence at high Mach numbers, due to the negative skew of the PDF. If the degree of non-lognormality of star forming gas varies as a function of time, it means that the average star formation efficiency can change by a factor of ∼ 2. The role of this non-lognormality may be incorporated into more detailed models of star formation efficiency.

## 6. CONCLUSIONS

3 Various treatments of the critical density in the literature may vary with order unity factors, which we ignore. Alternatively one may forego the use of a transition density and compare timescales τ ( s ) and t ff ( s ) in a Lagrangian frame. But in our Eulerian frame τ ( s ) also includes advection which does not oppose gravitational collapse. So here we just use an 'illustrative' transition density, rather than one computed directly from the model.

ϵ

Figure 5. Star formation efficiency (calculated as fraction of gas above critical density s crit ) vs. sonic Mach number. The efficiency is increased in our non-lognormal model for turbulence at high Mach number, due to the density PDF being negatively skewed. For comparison, we show efficiencies calculated from a lognormal model and the 'lognormal plus ( α = 2) powerlaw' model of Burkhart &amp; Mocz (2019). The multi-free fall time factor (Federrath &amp; Klessen 2012) is not included in this comparison, which can reverse the trend with Mach number.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00544v2-figure-a9023e87-622d-403f-98a1-ea65fbe12433.jpg)

We have developed a simple Markov model with an analytic solution to describe the non-lognormal density distribution in supersonic, isothermal turbulence. The model was calibrated against 1D 'turbulent' shock simulations, and shown to fit PDFs of 3D supersonic turbulence as well. Our main conclusions are as follows:

- · High density structures in supersonic turbulence are produced by strong shocks that have shorter transient timescales than low density structures.
- · These strong shocks lead to deviations from Gaussianity in the natural log of the density PDF, particularly an increasing negative skewness as a function of turbulent Mach number.
- · We have characterized these features of supersonic turbulence with a simple Markov model.
- · Our results can be applied to both 1D (shocks only) and 3D simulations (shocks plus turbulence) with various driving methods (compressive vs solenoidal). This strongly suggests that the simplified 1D picture of shocks holds in 3D (which is a system of quasi-1D, interacting shocked sheet fragments with turbulence).
- · Star formation efficiencies calculated from density PDFs are expected to increase at high Mach numbers when the PDF exhibits strong deviations from lognormal.

- · Our model also applies to subsonic compressive turbulence, where the model reduces to a near lognormal density distribution.

density PDF is used to predict star formation efficiencies.

Our model provides a simple way to quantify the nonlognormality of density structures in turbulent environments, such as star forming molecular clouds. It can be applied to factor in this feature in studies where the

The authors would like to thank Eve Ostriker, Jim Stone, and Mohammad Safarzadeh for valuable discussions. Support for this work (P.M.) was provided by NASA through Einstein Postdoctoral Fellowship grant number PF7-180164 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS803060. B.B. thanks the Simons Foundation for generous support.

## REFERENCES

| Burkhart, B. 2018, ApJ, 863, 118                                                        |
|-----------------------------------------------------------------------------------------|
| Burkhart, B., Falceta-Gon¸alves, c D., Kowal, G., & Lazarian, A. 2009, ApJ, 693, 250    |
| Burkhart, B., & Lazarian, A. 2012, ApJL, 755, L19                                       |
| Burkhart, B., & Mocz, P. 2019, ApJ, 879, 129                                            |
| Falceta-Gon¸alves, c D., Lazarian, A., & Kowal, G. 2008, ApJ, 679, 537                  |
| Federrath, C. 2013, MNRAS, 436, 1245                                                    |
| Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156                                     |
| Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79                        |
| Hennebelle, P., & Chabrier, G. 2011, ApJL, 743, L29                                     |
| Hopkins, P. F. 2013, MNRAS, 430, 1880                                                   |
| Konstandin, L., Schmidt, W., Girichidis, P., et al. 2016, MNRAS, 460, 4483              |
| Kowal, G., Lazarian, A., & Beresnyak, A. 2007, ApJ, 658, 423                            |
| Krumholz, M. R., Dekel, A., & McKee, C. F. 2012, ApJ, 745, 69                           |
| Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 Larson, R. B. 1981, MNRAS, 194, 809 |
| Mocz, P., Burkhart, B., Hernquist, L., McKee, C. F., &                                  |
| Mocz, P., & Burkhart, B. 2018, MNRAS, 480, 3916                                         |
| Springel, V. 2017, ApJ, 838, 40                                                         |

| Molina, F. Z., Glover, S. C. O., Federrath, C., & Klessen, R. S. 2012, MNRAS, 423, 2680               |
|-------------------------------------------------------------------------------------------------------|
| Nordlund, ˚. A K., & Padoan, P. 1999, in Interstellar Turbulence, ed. J. Franco & A. Carraminana, 218 |
| Novikov, E. 1989, Physics of Fluids A: Fluid Dynamics, 1, 326                                         |
| Padoan, P., & Nordlund, ˚. A 2011a, ApJ, 730, 40 -. 2011b, ApJL, 741, L22                             |
| Pan, L., Padoan, P., & Nordlund, ˚. A 2019, arXiv e-prints, arXiv:1905.00923                          |
| Passot, T., & V´zquez-Semadeni, a E. 1998, PhRvE, 58, 4501                                            |
| Pedrizzetti, G. 1999, Physics of Fluids, 11, 1694                                                     |
| Pedrizzetti, G., & Novikov, E. A. 1994, Journal of Fluid Mechanics, 280, 69                           |
| Price, D. J., Federrath, C., & Brunt, C. M. 2011, ApJL, 727, L21                                      |
| Renner, C., Peinke, J., & Friedrich, R. 2001, Journal of Fluid Mechanics, 433, 383                    |
| Robertson, B., & Goldreich, P. 2018, ApJ, 854, 88                                                     |
| Scannapieco, E., & Safarzadeh, M. 2018, ApJL, 865, L14                                                |
| Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23                                            |
| Springel, V. 2010, MNRAS, 401, 791                                                                    |
| Squire, J., & Hopkins, P. F. 2017, MNRAS, 471, 3753                                                   |
| Uhlenbeck, G. E., & Ornstein, L. S. 1930, Physical Review, 36, 823                                    |
| Vazquez-Semadeni, E. 1994, ApJ, 423, 681                                                              |
Philip Mocz Blakesley Burkhart
A Markov model for non-lognormal density distributions in compressive isothermal turbulence
Compressive isothermal turbulence is known to have a near lognormal density
probability distribution function (PDF) with a width that scales with the sonic
Mach number and nature of the turbulent driving (solenoidal vs compressive).
However, the physical processes that mold the extreme high and low density
structures in a turbulent medium can be different, with the densest structures
being composed of strong shocks that evolve on shorter timescales than the low
density fluid. The density PDF in a turbulent medium exhibits deviations from
lognormal due to shocks, that increases with the sonic Mach number, which is
often ignored in analytic models for turbulence and star formation. We develop
a simple model for turbulence by treating it as a continuous Markov process,
which explains both the density PDF and the transient timescales of structures
as a function of density, using a framework developed in n Scannapieco &
Safarzadeh (2018). Our analytic model depends on only a single parameter, the
effective compressive sonic Mach number, and successfully describes the
non-lognormal behavior seen in both 1D and 3D simulations of supersonic and
subsonic compressive isothermal turbulence. The model quantifies the
non-lognormal distribution of density structures in turbulent environments, and
has application to star forming molecular clouds and star formation
efficiencies.