Gravo-Turbulent Fragmentation

# Turbulent Disks are Never Stable: Fragmentation and Turbulence-Promoted Planet Formation

Philip F. Hopkins11affiliation: TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA; E-mail:phopkins@caltech.edu 22affiliation: Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720 & Jessie L. Christiansen33affiliation: SETI Institute/NASA Ames Research Center, M/S 244-30, Moffett Field, CA 94035
###### Abstract

A fundamental assumption in our understanding of disks is that when the Toomre parameter , the disk is stable against fragmentation into smaller, self-gravitating objects (and so cannot, for example, form planets via direct collapse). However, if disks are turbulent, this criterion neglects a broad spectrum of stochastic density fluctuations that can produce rare but extremely high-density local mass concentrations that will easily collapse. We have recently developed an analytic framework to predict the statistics of these fluctuations. Here, we use these models to consider the rate of fragmentation and mass spectrum of fragments formed in a turbulent, Keplerian disk (e.g. a proto-planetary or proto-stellar disk). Turbulent disks are never completely stable: we calculate the (always finite) probability of the formation of self-gravitating structures via stochastic turbulent density fluctuations (compressions, shocks, etc.) in such a disk. Modest sub-sonic turbulence above a minimum Mach number is sufficient to produce a few stochastic fragmentation or “direct collapse” events over Myr timescales, even if and cooling is relatively “slow” (). In trans-sonic turbulence () this extends to as large as . We derive the “true” criterion needed to suppress such events (over some timescale of interest), which scales exponentially with Mach number. We specify this to cases where the turbulence is powered by MRI, convection, or density/spiral waves, and derive equivalent criteria in terms of and/or the disk cooling time. In the latter case, cooling times as long as may be required to completely suppress this channel for collapse. These gravo-turbulent events produce a mass spectrum concentrated near the Toomre mass (spanning rocky-to-giant planet masses, and increasing with distance from the star), with a wider mass spectrum as increases. We apply this to plausible models of proto-planetary disks with self-consistent cooling rates and temperatures, and show that, beyond radii au, no disk temperature can fully suppress stochastic events. For theoretically expected temperature profiles, even a minimum mass solar nebula could experience stochastic collapse events, provided a source of turbulence with modest sub-sonic Mach numbers.

protoplanetary discs — planets and satellites: formation — accretion, accretion disks — hydrodynamics — instabilities — turbulence — star formation: general — galaxies: formation
slugcomment: Submitted to ApJ, May, 2013

## 1. Introduction

### 1.1. The Problem: When Do Disks Fragment?

Fragmentation and collapse of self-gravitating gas in turbulent disks is a process central to a wide range of astrophysics, including planet, star, supermassive black hole, and galaxy formation. The specific case of a “marginally stable,” modestly turbulent Keplerian disk is particularly important, since this is expected in proto-planetary and proto-stellar disks as well as AGN accretion disks. There has been considerable debate, especially, regarding whether planets could form via “direct collapse” – fragmentation of a self-gravitating region in a proto-planetary disk – as opposed to accretion onto planetesimals (see e.g. Boley et al., 2006; Boley, 2009; Dodson-Robinson et al., 2009; Cai et al., 2008, 2010; Vorobyov & Basu, 2010; Johnson et al., 2010; Boss, 2011; Stamatellos et al., 2011; Seager et al., 2011; Forgan & Rice, 2013, and references therein). Even if this is not a significant channel for planet formation, it is clearly critical to understand the conditions needed to avoid fragmentation. This is especially demanding because such fragments need only form an order-unity number of times over the millions of dynamical times most proto-planetary disks survive, to account for a significant fraction of planets.

Yet there is still no consensus in the literature regarding the criteria for “stability” versus fragmentation in any disk, especially in nearly-Keplerian proto-planetary disks. The classic Toomre criterion:

 Q≡σgκπGΣgas>1 (1)

where is the gas velocity dispersion, the epicyclic frequency, and the gas surface density (Toomre, 1964, 1977; Goldreich & Lynden-Bell, 1965) appears to offer some guidance. And indeed, for , most of the mass in disks fragments into self-gravitating clumps in roughly a single crossing time. But this was derived for a smooth, homogeneous disk, dominated by thermal pressure with no cooling, so does not necessarily imply stability in any turbulent system. Gammie (2001) studied a more realistic case of a turbulent disk with some idealized cooling, and showed that if the cooling time exceeded a couple times the dynamical time, the disk could maintain its thermal energy (via dissipation of the turbulent cascade), with a steady-state and transsonic Mach numbers (powered by local spiral density waves), thus avoiding runway catastrophic fragmentation.

However, many subsequent numerical simulations have shown that, although catastrophic, rapid fragmentation may be avoided when these criteria are met, simulations with larger volumes (at the same resolution) and/or those run for longer timescales still eventually form self-gravitating fragments, even with cooling times as long as times the dynamical time (see references above and e.g. Rice et al., 2005, 2012; Meru & Bate, 2011b, a; Meru & Bate, 2012; Paardekooper et al., 2011; Paardekooper, 2012). Increasing spatial resolution and better resolution of the turbulent cascade (higher Reynolds numbers) appear to exacerbate this – leading to the question of whether there is “true” stability even at infinitely long cooling times (Meru & Bate, 2012). And these simulations are still only typically run for a small fraction of the lifetime of such a disk (or include only a small fraction of the total disk mass, in shearing-sheet models), and can only survey a tiny subsample of the complete set of parameter space describing realistic disks.

### 1.2. The Role of Turbulent Density Fluctuations

Clearly, the theory of disk fragmentation requires revision. But almost all analytic theoretical work has assumed that the media of interest are homogeneous and steady-state (though see Kratter & Murray-Clay, 2011, and references above), despite the fact that perhaps the most important property of turbulent systems is their inhomogeneity. In contrast to the “classical” homogenous models, Paardekooper (2012) went so far as to suggest that fragmentation when may be a fundamentally stochastic process driven by random turbulent density fluctuations – and so can never “converge” in the formal sense described by the simulations above. That said, simulations of idealized turbulent systems in recent years have led to important breakthroughs; in particular, the realization that compressible, astrophysical turbulence obeys simple inertial-range velocity scalings (e.g. Ossenkopf & Mac Low, 2002; Federrath et al., 2010; Block et al., 2010; Bournaud et al., 2010), and that – at least in isothermal turbulence – the density distribution driven by stochastic turbulent fluctuations develops a simple shape, with a dispersion that scales in a predictable manner with the compressive Mach number (Vazquez-Semadeni, 1994; Padoan et al., 1997; Scalo et al., 1998; Ostriker et al., 1999).

Recently, in a series of papers (Hopkins, 2012b, 2013b, c, d, 2013a), we showed that the excursion-set formalism could be applied to extend these insights from idealized simulations, and analytically calculate the statistics of bound objects formed in the turbulent density field of the ISM. This is a mathematical formulation for random-field statistics (i.e. a means to use the power spectra of turbulence to predict the statistical real-space structure of the density field), well known from cosmological applications as a means to calculate halo mass functions and clustering in the “extended Press-Schechter” approach from Bond et al. (1991). This is a well-known theoretical tool in the study of large scale structure and galaxy formation, and underpins much of our analytic understanding of halo mass functions, clustering, mergers, and accretion histories (for a review, see Zentner, 2007). The application to turbulent gas therefore represents a means to calculate many interesting quantities analytically that normally would require numerical simulations. In Hopkins (2012b) (hereafter Paper I), we focused on the specific question of giant molecular clouds (GMCs) forming in the ISM, and considered the simple case of isothermal gas with an exactly lognormal density distribution. We used this to predict quantities such as the rate of GMC formation and collapse, their mass function, size-mass relations, and correlation functions/clustering, and showed that these agreed well with observations. In Hopkins (2013a) (hereafter Paper II), we generalized the models to allow arbitrary turbulent power spectra, different degrees of rotational support, non-isothermal gas equations of state, magnetic fields, intermittent turbulence, and non-Gaussian density distributions; we also developed a time-dependent version of the theory, to calculate the rate of collapse of self-gravitating “fragments.”

### 1.3. Paper Overview

In this paper, we use the theory developed in Paper I & Paper II to calculate the statistics of fragmentation events in Keplerian, sub and trans-sonically turbulent disks, with a particular focus on the question of fragmentation in proto-planetary disks. We develop a fully analytic prediction for the probability, per unit mass and time, of the formation of self-gravitating fragments (of a given mass) in a turbulent disk. In addition to providing critical analytic insights, this formulation allows us to simultaneously consider an enormous dynamic range in spatial, mass, and timescale, and to consider extremely rare fluctuations (e.g. fluctuations that might occur only once over millions of disk dynamical times), which is impossible in current numerical simulations.111It is important to clarify, when we refer to “large fluctuations,” are not referring to extremely large, single-structure “forcing” events (e.g. a very strong shock on large scales). In fact, we assume that the probability of such “positive intermittency events” (isolated, large amplitude Fourier modes) is vanishingly small (, in the language of the intermittency models discussed in § 4 and Hopkins 2013a, 2012a). What we calculate is the (rare) probability of many small, independent fluctuations on different scales acting, by random chance, sufficiently “in phase” so as to produce a large density perturbation. If “positive intermittency events” do occur, it may significantly increase the probability of rare collapse events. We will show that this predicts “statistical” instability, fragmentation, and the formation of candidate “direct collapse” planets/stars even in the “classical” stability regime when and cooling is very slow. However, these events occur stochastically, separated by much larger average timescales than in the catastrophic fragmentation which occurs when . We will show that this naturally explains the apparently discrepant simulation results above, and may be important over a wide dynamic range of realistic disk properties for formation of both rocky and giant planets.

## 2. “Classical” vs. “Statistical” Stability

Motivated by the discussion above, and the models developed in this paper, we will introduce a distinction between two types of (in)stability: “classical” and “statistical.”

By “classical” instability, we refer to the traditional regime where , in which small perturbations to the disk generically grow rapidly. A “classically unstable” disk will fragment catastrophically, with a large fraction of the gas mass collapsing into self-gravitating objects on a few dynamical times. A “classically stable” disk will survive many dynamical times in quasi-steady state, and gas at the mean density will not be self-gravitating.

By “statistical” instability, we refer to the regime where an inhomogeneous disk experiences sufficiently large fluctuations in density such that there is an order-unity probability (integrated over the entire volume and lifetime of the disk) of the formation of some region which is so over-dense that it can successfully collapse under self-gravity. A “statistically stable” disk has a probability much less than unity of such an event occurring, even once in its lifetime.

Classically unstable disks are always statistically unstable, and statistically stable disks are always classically stable. However, we argue in this paper that there is a large regime of parameter space in which disks can be classically stable, but statistically unstable. In this regime, disks can, in principle, evolve for millions of dynamical times with , and nearly all the disk mass will be stable, but rare density fluctuations might form an order-unity number of isolated self-gravitating “fragments.”

## 3. Model Outline

Here, we present an order-of-magnitude, qualitative version of the calculation which we will develop rigorously below. This serves to illustrate some important scalings and physical processes.

Consider an inhomogeneous disk (around a star of mass ) with scale radius , scale-height , mean surface density and mid-plane density , and Toomre . Although the disk is classically stable (and so not self-gravitating at the mean density), if a local “patch” of the disk exceeds some sufficiently large density then that region will collapse under self-gravity. The most unstable wavelengths to self-gravity are of order the disk scale height : roughly speaking, a gas parcel of this size will collapse under self-gravity if it exceeds the Roche criterion (overcomes tidal forces): .222To derive this, we assume a Keplerian disk, , and vertical equilibrium, for a disk supported by thermal pressure.

So for we see that gas at the mean density will not collapse. However, turbulence produces a broad spectrum of density fluctuations. For simple, isothermal turbulence with characteristic (compressive) Mach number , the (volumetric) distribution of densities is approximately log-normal, with dispersion ( for ).

Since density fluctuations exceeding () can collapse, we integrate the tail of the log-normal density probability distribution function (PDF) above this critical density to estimate the probability , per unit volume, of a region being self-gravitating at a given instant. Since we are considering regions of size , the disk contains approximately independent volumes, and so (assuming is small) the probability of any one such volume having a coherent volume-average density above the critical threshold is .

Now consider that a typical proto-planetary disk (at a few AU) might have and , with . In this case, and are extremely small!

However, this analysis applies to the disk viewed at a single instant. Turbulent density fluctuations evolve stochastically in time, with a coherence time (on a given scale) about equal to the bulk flow crossing time (since this is the time for fluid flows to cross and interact with a new independent region of size ). So the density PDF is re-sampled or “refreshed” on the timescale . But this is just a dynamical time, which is very short relative to a typical disk lifetime ( yr at  AU). If the disk survives for a total timescale , then, the entire volume is re-sampled independent times. So the probability, integrated over time, of any one of these volumes, at any one time, exceeding the self-gravity criterion, is .

If a typical disk with parameters above survives for  Myr, or crossing times, we then obtain a time-and-volume-integrated probability of at least a single stochastic “fragmentation event” driven by turbulent density fluctuations. The mass of the self-gravitating “fragment” will be (, for these parameters in a minimum mass solar nebula). And, despite the fact that the average timescale between such events may be long ( Myr), if a fragment forms, it forms rapidly (in  yr) on the turbulent crossing time  yr and has a short collapse/free-fall time yr. Despite having , then, we estimate that this disk could be statistically unstable!

For otherwise fixed and , declines exponentially with increasing . So for , and such a disk is both classically and statistically stable. For , the disk is classically unstable, and even material at the mean density (i.e. an order-unity fraction of the mass) begins collapsing on a single dynamical time. But for , such a disk is classically stable, but statistically unstable.

Below, we present a more formal and rigorous derivation of statistical (in)stability properties, and consider a range of disk properties. But the simple order-of-magnitude arguments above capture the most important qualitative behaviors we will discuss.

## 4. The Model: Turbulent Density Fluctuations

As mentioned above, turbulence in approximately isothermal gas (neglecting self-gravity) generically drives the density PDF to a approximate lognormal (a normal distribution in ). This follows from the central limit theorem (see Passot & Vazquez-Semadeni, 1998; Nordlund & Padoan, 1999), although there can be some corrections due to intermittency and mass conservation (Klessen, 2000; Hopkins, 2012a). And in the simplest case of an ideal box of driven turbulence, the variance simply scales with the driving-scale Mach number as (where is a constant discussed below). These scalings have been confirmed by a huge number of numerical experiments with , sampling the PDF down to values as low as in the “tails” (Federrath et al., 2008, 2010; Federrath & Klessen, 2012; Schmidt et al., 2009; Price & Federrath, 2010; Konstandin et al., 2012).

This is true in both sub-sonic and super-sonic turbulence (see references above), as well as highly magnetized media (with the magnetic fields just modifying or the effectively compressible component of ; see Kowal et al. 2007; Lemaster & Stone 2009; Kritsuk et al. 2011; Molina et al. 2012), and even multi-fluid media with de-coupled electrons, ions, neutral species, and dust (see Downes, 2012).333It is a common misconception that log-normal density PDFs apply only to super-sonic, non-magnetized turbulence. In fact, while they apply strictly to only isothermal (non-intermittent) turbulence (as discussed in § 4), the analytic derivation of lognormal PDFs actually assumes small (local) Mach numbers (see Nordlund & Padoan 1999). The log-normal model (with the higher-order intermittency corrections we include in § 5) and the simple assumptions we use for the scaling of the density power spectrum with velocity power spectrum have been tested in the sub-sonic limit in both simulations (see Shaikh & Zank, 2007; Kowal et al., 2007; Burkhart et al., 2009; Schmidt et al., 2009; Konstandin et al., 2012), and also in experimental data from the solar wind (Burlaga, 1992; Forman & Burlaga, 2003; Leubner & Vörös, 2005) and laboratory MHD plasmas (Budaev, 2008) as well as jet experiments (Ruiz-Chavarria et al., 1996; Warhaft, 2000; Zhou & Xia, 2010). The power spectrum predictions in this limit, in fact, generically follow the well-known and tested weakly-compressible Kolmogorov-like scaling (Montgomery et al., 1987; Zank & Matthaeus, 1990). And in Hopkins (2012a), we show that deviations from lognormal statistics in the high-density wing of the distribution are less significant at lower Mach numbers. Likewise, because of the lower compressive Mach numbers, these assumptions are more accurate in solenoidally-forced and/or magnetized turbulence (see Kowal et al., 2007; Hopkins, 2012a; Federrath, 2013). And although the density PDF is not exactly lognormal when the gas is no longer isothermal, the inviscid Navier-Stokes equations show that the local response is invariant under the substitution , allowing an appropriately modified PDF which provides an excellent approximation to the results in simulations (see Scalo et al., 1998; Passot & Vazquez-Semadeni, 1998).

In Paper I-Paper II, we use these basic results to show how excursion set theory can be used to analytically predict the statistical structure of turbulent density fluctuations. The details are given therein; for the sake of completeness, we include a summary of the important equations in Appendix B. Here, we briefly describe the calculation.

Consider the field : the density field about the (random) coordinate x in space, smoothed with some window function of characteristic radius (e.g. the average density in a sphere). If the gas is isothermal (we discuss more complicated cases below), this is distributed as a log-normal:

 P(ln[ρx(R)])=1√2πS(R)exp(−[lnρx(R)+S(R)/2]22S(R)) (2)

where is the variance on each scale . is just the Fourier transform of the density power spectrum, which itself follows from the (well-defined) velocity power spectrum (see Eq. B2). Essentially, is determined by integrating the contribution to the variance on each scale from the velocity field, using the relation ( is the scale-dependent Mach number of the velocity field).

“Interesting” regions are those where exceeds some critical value, above which the region is sufficiently dense so as to be self-gravitating. Including the effects of support from angular momentum/shear, thermal and magnetic pressure, and turbulence, this is given by

 ρcrit(R)ρ0≡Q2~κ(1+hR)[σ2g(R,ρcrit)σ2g(h,ρ0)hR+~κ2Rh] (3)

(Vandervoort, 1970). Here is the mean midplane density of the disk, is the disk scale height, for a Keplerian disk ( is the epicylic frequency), and is the Toomre parameter. The effective gas dispersion (Eq. B5) includes the thermal (), turbulent (), and magnetic support (Alfvén speed ). A full derivation is given in Paper I; but we stress that this not only implies/requires that a region with is gravitationally self-bound, but also that such regions will not be unbound by tidal shear (i.e. have sizes within the Hill radius) and that they will not be unbound/destroyed by energy input from the turbulent cascade (shocks and viscous heating).444Note that in Eq. 3 we consider the rms turbulent velocity (i.e. do not explicitly treat different velocity fluctuations), even though we consider density fluctuations from the turbulence. In Paper I & Paper II we show that this is a very small source of error (see also e.g. Sheth & Tormen 2002, who show that even for pressure-free flows with velocities typically above escape velocities, this introduces a correction to the predicted mass function). But particularly here, since the turbulence of interest is sub-sonic, is never the dominant component of . Moreover, even though velocity fluctuations do drive the density fluctuations, because these are built up hierarchically in the cascade, the magnitude of the coherent velocity fluctuations on a given scale is instantaneously nearly uncorrelated with the density fluctations (see Federrath et al. 2010, who find ). Inserting such a scaling into our model is straightforward, but has no detectable effect in any figure shown. Finally, the most important effects of coupled velocity-density fluctuations are already included in the models for intermittency we consider, since this coupling makes the PDF non-lognormal. Knowing the size and critical density of a region, it is trivial to translate this to the total collapsing gas mass .

We desire the mass and initial size spectrum of regions which exceed (so are self-gravitating) on the scale specifically defined as the largest scale on which the region is self-gravitating (i.e. excluding “sub-units” so that we do not double-count “clouds within clouds”). In Paper I we show this reduces to a derivation of the “first-crossing” distribution for the field . The mass function of collapsing objects can be written generally as

 dndM=ρcrit(M)Mff(M)∣∣dSdM∣∣ (4)

where is a function given in Eq. B8 that is somewhat cumbersome to derive (see Paper I for details), but depends only on how the dimensionless quantities and “run” as a function of scale .

In Paper II, we further generalize this to fully time-dependent fields. In statistical equilibrium, the density field obeys a modified Fokker-Planck equation with different modes evolving stochastically – they follow a damped random walk with a correlation time equal to the turbulent crossing time for each spatial scale/wavenumber. This allows us to directly calculate the probability per unit time of the formation of any bound object in the mass function above. The exact solution requires a numerical approach described in Paper I (§ 7) and Paper II (§ 9-10), but to very good approximation (for ) this is just

 ΔP(M,Δt)≈dndMΔtτ(M) (5)

where is the turbulent crossing time on the scale corresponding to the mass (Eq. B6). This is because the coherence time of density fluctuations on a given scale is just the crossing time.

Now consider a disk, or disk element (if we consider a series of cylindrical annuli at different disk-centric radii) with total mass . By the definitions used in Paper I-Paper II, this means the total “effective volume” is , i.e. that the total (integrated) number of objects per unit mass is

 (6)

For any polytropic gas, we can factor out all dimensional parameters, and work in units of and . Mass then has units of . But since, for a disk in vertical equilibrium , and for a Keplerian disk (where is the central mass and is the disk-centric radius), we have ). So we can write where

 μ≡πΣgasr2∗M∗≈MdM∗ (7)

where is the disk mass within . For an exponential vertical profile used to define our dispersion relation, , so we can use this and to link . We can then remove the local quantities and and define mass in the much simpler global units of and . The units of defined above are , so this can be similarly re-written. And since the disk crossing time scales as , provides a natural time unit.

Having defined units, the model is completely specified by dimensionless parameters. These are the spectral index of the turbulent velocity spectrum, (), and its normalization, which we define by the Mach number on large scales , as well as the Toomre parameter . We must also specify the parameter , i.e. the mean fraction of the velocity in compressive modes. This is for purely compressively forced turbulence, for purely solenoidally forced turbulence, and for random forcing. But we can almost completely factor out the dependence on this parameter if we simply define the compressive component of the turbulence, i.e. work in units of the compressive Mach number555At sufficiently small , the compressive-to-total Mach number ratio can scale steeply with Mach number depending on the turbulent forcing (see Zank & Matthaeus, 1990; Zank et al., 1990). Because we work specifically with the compressive Mach number, this is largely irrelevant to our calculation (what does matter is that the relation between and density fluctuations remains intact; see Konstandin et al. 2012). Moreover the “steepening” becomes significant only below the minimum values we will identify as interesting (compare e.g. Fig. 6 of Konstandin et al. 2012), especially for magnetized turbulence (where other effects have the opposite sense; see Ostriker et al. 2001; Shaikh 2007; Price et al. 2011; Molina et al. 2012). .

In Paper II, we show how these equations generalize for the cases with non-isothermal gas, intermittent turbulence, and non-isotropic magnetic fields. The qualitative scalings are similar, but the math becomes considerably more complicated (and the “first-crossing distribution” and its evolution in time must be solved via a numerical Monte Carlo method rather than analytically). The presence of intermittency makes the density PDF non-Gaussian, and introduces explicitly correlated fluctuation structures, but this can be entirely encapsulated in a modified form of Eq. B8 above (see Paper II, Appendix D), leaving the rest of our derivation intact. We will consider such non-Gaussian, correlated statistics and show they give very similar results. For non-isothermal cases, we must replace (in calculating both the variance and critical density), and again allow for the density PDF to be non-Gaussian (with a skew towards lower/higher densities as the equation of state is made more or less “stiff,” respectively; see § 3-4 in Paper II; this also introduces higher-order correlations in the fluctuation statistics, distinct from those associated with intermittency). Modulo these changes, however, our derivations (and the changes made to specify to the Keplerian disk case) are identical for the simple case of a polytropic gas where ( is the polytropic index). Magnetic fields can produce global anisotropy, but this can be simply absorbed into the form of as well; their largest effect in suppressing fluctuations manifests as a field-strength-dependent value (Paper II, § 6, and e.g. Kowal et al. 2007; Lemaster & Stone 2009; Molina et al. 2012). For the strong-field limit, however, this is just similar to the pure-solenoidal turbulence case (in both cases, there is only one spatial dimension along which compression is possible); so is within the range of variations we will consider.666We do caution that some of our simple assumptions (for example, how we assume the density and velocity power spectra “turn over” at large scale heights ) remain to be tested in simulations of fully nonlinear, shearing, vertically stratified, MHD turbulence. However, within the disk scale height, preliminary comparisons suggest these corrections have little effect. In § 6.2 below, we explicitly compare the predictions from our model to the variance in the mid-plane density field calculated in such simulations (vertically stratified, shearing MRI boxes; Bai & Stone 2012), and find our predictions agree remarkably well over a range of field strengths. Moreover, we can repeat this comparison at various vertical heights up to , and find agreement. Recently, Arena & Gonzalez (2013) have performed numerical experiments of three-dimensional (non-magnetized, isothermal) disks with effective Mach numbers from to (); the results agree well with our assumed form for the density distribution and power spectrum shape. In future work (J. Lynn et al., in preparation), we will study the effect of non-isothermal, stratified, rotating, and self-gravitating flows on turbulent density fluctuations.

## 5. Results: Fragmentation Rates and Mass Spectra in the General Case

In Fig. 1, we now use this to predict the mass spectrum – specifically the probability per unit time, per unit mass – of the formation of self-gravitating regions, as a function of various properties of the system.777For simplicity, we will refer to these calculations as if the disk is constant surface-density out to some maximum radius, with total mass and (the orbital velocity) defined at that radius, and constant Toomre parameter. However, the results shown could be applied to any radius in a disk with any mass profile (provided it is still Keplerian), with defined as the disk mass inside a radial annulus, and and evaluated at that same radius.

We define our “reference” model to be isothermal (), non-intermittent (i.e. log-normal density PDF), with turbulent spectral slope (appropriate for compressible turbulence), Toomre , and have (random forcing) with rms three-dimensional compressive Mach number on the largest scales . But we then vary all of these parameters. For now, we treat each as free – in other words, we make no assumptions about the specific mass profile, temperature structure, cooling chemistry, or other microphysics of the disk. These microphysics are, of course, what ultimately determine the value of the model parameters (and may build in some intrinsic correlations between them); but for any given set of parameters in Fig. 1, the prediction is independent of how the microphysics produce those parameters.

In general, the shape of the mass spectrum is similar regardless of these variations. This is peaked (very approximately lognormal-like) around a characteristic mass . For a relatively massive disk , then, we expect characteristic masses in such events of , corresponding to gas giants. If however the disk mass is lower, , then this becomes , typical of rocky (Earth-mass) planets.

The dependence on the slope of the turbulent spectrum is quite weak: Kolmogorov-like (, appropriate for incompressible turbulence) spectra give nearly identical results; shallower spectra slightly broaden the mass range predicted (since declines more slowly at small scales), but these are not seen in realistic astrophysical contexts. Likewise, at fixed , there is some dependence on (because of how enters into the critical density for collapse), but this is small in this regime, because turbulence is not the dominant source of “support” resisting collapse. The effects of intermittency are very weak, since they only subtly modify the density PDF shape (see Paper II; the important point here is that our predictions do not much change for reasonable departures from Gaussian or log-normal statistics). Even changing the equation of state has surprisingly mild effects. A “stiffer” (higher-) equation of state is more resistive to fragmentation on small scales and leads to a more sharply peaked spectrum. But for fixed , the variance near the “core” of the density PDF is similar independent of , so this does not have much effect on the normalization of the probability distribution. We stress that the medium could have a perfectly adiabatic equation of state with no cooling, and if it had the plotted Mach number and Toomre , our result would be identical (and the probability of fragmentation would still be finite). Very soft equations of state, on the other hand, can lead to a runaway tail of small-scale fragmentation, but this is not likely to be relevant. Clearly, the largest effects come from varying and ; at or the mass distribution rapidly becomes more broad, while at or the characteristic mass remains fixed but the normalization of the probability becomes exponentially suppressed.

This is summarized in Fig. 2. Since the mass range of expected “events” is relatively narrow, we integrate over mass to obtain the total probability of an event per unit time,

 dNfragdt=∫dlogMdNdlogMdt (8)

and plot this as a function of for different parameter choices. As before, most parameters make a surprisingly small difference; and dominate.

### 5.1. A General Statistical Stability Criterion

Formally is always non-zero for ; but we see that for for , the probability per unit time of forming a self-gravitating fluctuation drops rapidly. However, recall that the total lifetime of e.g. a proto-planetary disk is many, many disk dynamical times , and our “time unit” is . Consider a typical lifetime of ; then the disk at  AU around a solar-mass star (where ) experiences dynamical times. For a disk-to-total mass ratio of , this is “time units,” so if we integrate the probability of a fragmentation event over the lifetime of the disk we obtain an order-unity probability even for in the units here (i.e.  as small as ). Of course, following such an integration in detail requires knowing the evolution of the disk mass, , , etc. But we can obtain some estimate of the value of required for statistical stability (ensuring the probability of fragmentation events is negligible) by simply assuming all quantities are constant and integrating over an approximate timescale, also shown in Fig. 2. Here we take our standard model, and consider three timescales in units of (a factor of shorter and longer than the value motivated above), and consider the minimum at each needed to ensure that the time-integrated probability of a fragmentation event is .

This minimum is at ; equivalently, disks with and have an order-unity probability of at least one stochastic fragmentation event over their lifetime (). At larger , is required for ; and by sonic Mach numbers , is required for .

We can approximate the scaling of by the following: recall that the critical density near the Toomre scale scales approximately as (Eq. 3), while the density dispersion scales as (Eq. B2). As noted above, if the system evolves for a total timescale (time in our dimensionless units), then an event with probability per unit time has an order-unity probability of occurring. If the probabilities are approximately normally distributed then this is just , where is the barrier and the variance. Since the mass function is peaked near the Toomre scale we can approximate both by their values near the “driving scale” , and . Thus, statistical stability over some timescale of interest requires a in Fig. 2 of

 Qmin ≈0.5exp[√2ln(t0μ−2Ω)ln(1+M2c)] (9)

For typical values of , , and – i.e. the typical number of independent realizations (in both time and space) of the turbulent field in a protoplanetary disk, this becomes

 Qmin∼0.5exp(6√ln(1+M2c)) (10)

In other words, a event has order-unity chance of occuring once over the disk lifetime, so for any this implies a minimum needed to ensure statistical stability in such an extreme event.

## 6. How Is the Turbulence Powered? Statistical Stability in Specific Models for Turbulence

Thus far, we have considered the general case, varying the Mach numbers independent of other disk properties such as , , and . However, in a realistic physical model, the mechanisms that drive turbulence may be specifically tied to these properties. Moreover, there may be certain characteristic Mach numbers expected or ruled out. In this section, we therefore consider some well-studied physical scenarios for the driving of turbulence in Keplerian disks, and examine their implications for the “statistical stability” we have described above.

### 6.1. The “Gravito-turbulent” Regime (Gravity-Driven Turbulence)

Much of the work studying fragmentation in Keplerian disks has considered disks with a (locally) constant-cooling rate (“” disks, with locally fixed). In particular, this includes the scenario of a “gravito-turbulent” steady-state from Gammie (2001), with local instabilities (spiral waves) powering turbulence which contributes an effective viscosity and maintains a steady temperature and . The theory we present above is more general than this: we make no assumption about the detailed cooling physics, or that the disk is an -disk, and allow the various parameters , , , etc. to freely vary, whereas many of these are explicitly linked in these models. But in the theory above we cannot predict these quantities (or their co-dependencies); therefore, this model provides a simple and useful way to relate and predict some of the otherwise independently free parameters of the more general case, and is worth considering in detail.

#### 6.1.1 General Scalings

Consider a cooling rate which is uniform over a region (annulus) of the disk,

 tcool=ζΩ−1 (11)

If dissipation of gravitational instabilities (e.g. spiral waves) provides a source of heating balancing cooling, and , the system can develop a quasi steady-state angular momentum transport and Toomre parameter, as in a Shakura & Sunyaev (1973) -disk. Pringle (1981) and Gammie (2001) showed that in this equilibrium, the “effective” viscosity parameter

 α≈491γ(γ−1)ζ (12)

is approximately constant. This corresponds to the amplitude of density waves , leading to a maximum above (hence minimum below) which the local and catastrophic fragmentation will occur (see Gammie, 2001; Rice et al., 2005; Cossins et al., 2009). In an -disk, the inflow rate is also determined, as

 ˙M=3παc2sΣgasΩ−1 (13)

so this also corresponds to a maximum “classically stable” inflow rate below which catastrophic fragmentation will not occur.

Implicitly, the relations above also define a steady-state Mach number. Recall, the dissipation of spiral instabilities is ultimately governed by the turbulent cascade. Since the turbulent dissipation rate is constant over scale in a Kolmogorov cascade, we can take it at the top level, (where here we consider the rate per unit area) and equate it to the cooling rate , giving . Equivalently we could have equated the Reynolds stress that leads to , with ; we obtain

 M2c≈32α≈231γ(γ−1)ζ (14)

again.888One subtlety here is that the hydrodynamic Reynolds stress, and/or the dissipation on small scales, is dominated by the longitudinal (compressive) component. So this relation actually determines , not necessarily . But for our purposes, this is particularly convenient, as it allows us to drop the term from our earlier derivation. We now see how this relates to the theory developed in this paper.

In the language here, increasing enters our theory by – as discussed in the text – changing the equilibrium balance between thermal and turbulent energy, i.e. the Mach number. As cooling becomes less efficient, maintaining the same needed to power the turbulence (since spiral waves are still generated) requires a smaller turbulent dispersion, hence makes the system “more stable.”

#### 6.1.2 A Statistical Stability Criterion

But this also suggests an improved statistical stability criterion, accounting not just for regions where but also for stochastic local turbulent density fluctuations. For a given , we have calculated (Fig. 2) the Mach number above which the system will be probabilistically likely to fragment on a given timescale. Eq. 14 allows us to translate this to a minimum . We can do this exactly by simply reading off the numerically calculated values, but we can also obtain an accurate analytic approximation by the following (see § 8). Recall that for a system which evolves for a total timescale (time in the dimensionless units we adopt), we obtain the approximate Eq. 9 for the needed to ensure :

 Q≥12exp[√2ln(τ0)ln(1+M2c)] (15)

We can invert this to find the maximum for statistical stability at a given :

 M2c≲exp[[ln(2Q)]22lnτ0]−1≈[ln(2Q)]22lnτ0 (16)

Where the second equality follows from the fact that (for the systems of interest) is almost always true.

Combining this with Eq. 14, we obtain:

 ζmin ≈43γ(γ−1)ln(τ0)[ln(2Q)]2 (17) =43γ(γ−1)ln[t0/(μ2Ω−1)][ln(2Q)]2 (18)

We can immediately see some important consequences. Because of the stochastic nature of turbulent density fluctuations, will never converge in time integration (assuming the disk can maintain steady-state mean parameters) – there is always a finite (although possibly extremely small) probability of a strong shock or convergent flow forming a region which will collapse rapidly. However, the divergence in time is slow (logarithmic). The critical also scales with just as in the Gammie (2001) case; as the equation of state is made “stiffer,” the Mach numbers and density fluctuations are suppressed so faster cooling (lower ) can be allowed without fragmentation. And scales inversely with , so indeed higher- disks are “more stable,” but there is no “hard” cutoff at a specific value.

We can turn this around, and estimate the typical timescale for the formation of an order-unity number of fragments at a given , obtaining

 (19)

As expected, this quickly becomes large for modest and/or . We stress, however, that this is a probabilistic statement. Although the mean timescale between fragment formation events might be millions of dynamical times, if and when individual fragments form meeting the criteria in the text, they do so rapidly – on of order a single crossing time.

Given our derivation of , what do we expect in realistic systems such as proto-planetary disks? For a physical disk with we should expect (depending on the gas phase), and in equilibrium ; and we should integrate over the entire lifetime of the disk, (with values motivated in § 5). We then expect

 ζintmin≈51γ(γ−1)1+0.05ln(τint0/108)(1+1.4lnQ)2 (20)

This is fairly sensitive to – note if instead – and weakly sensitive to for reasonable variations. But this implies that only disks with extremely slow cooling, , (corresponding to steady-state Mach numbers ) are truly statistically stable with over such a long lifetime.

#### 6.1.3 Comparison with Simulations

Now consider the parameter choices in some examples that have been simulated. Gammie (2001) considered the case with , and a steady-state , evolving their simulations for typical timescales (though they consider some longer-scale runs). Because these were two-dimensional shearing-sheet simulations, the appropriate is somewhat ambiguous, but recall by our definitions, and for the assumptions in Gammie (2001) their “standard” simulation corresponds to an (where we equate the “full disk size” to the area of the box simulated). Plugging in these values, then, we predict , in excellent agreement with the value found by trial of several values therein. Paardekooper (2012) considered very similar simulations but with and all sheets run for ; for this system we predict – again almost exactly their estimated “fragmentation boundary.” Meru & Bate (2012) and Rice et al. (2012) consider three-dimensional global simulations; here is well-defined, a more realistic is adopted, and the disks self-regulate at ; the simulations are run for a shorter time (where for convenience we defined at the effective radius of the disk, since it is radius-dependent, but this is where the mass is concentrated), giving , in very good agreement with where both simulations appear to converge (using either SPH or grid-based methods). This is also in good agreement with the earlier simulations in Rice et al. (2005), for (predicting , vs. their estimated ) and (predicting , vs. their estimated ). Of course, we should naturally expect some variation with respect to the predictions, since this is a stochastic process, but we do not find any highly discrepant results.

We should also note that convergence in the total fragmentation rate in simulations – over any timescale – requires resolving the full fragmentation mass distribution in Fig. 1. Unlike time-resolution above, this is possible because there is clearly a lower “cutoff” in the mass functions (they are not divergent to small mass), but requires a mass resolution of (depending on the exact parameters). This is equivalent to a spatial resolution of , i.e. a small fraction of the disk scale-height . This also agrees quite well with the spatial/numerical resolution where (at fixed time evolution) many of the studies above begin to see some convergence (e.g. Meru & Bate, 2012; Rice et al., 2012), but it is an extremely demanding criterion.

### 6.2. The Magneto-Rotational Regime

In the regime where the disk is magnetized and ionized, the magneto-rotational instability (MRI) can develop, driving turbulence even if the cooling rate is low and . We therefore next consider the simple case where there is no gravo-turbulent instability (), but the MRI is present.

#### 6.2.1 General Scalings

Given MRI and no other driver of turbulence, Alfvén waves will drive turbulence in the gas to a similar power spectrum to the hydrodynamic case (within the range we examine where the power spectrum shape makes little difference), with driving-scale rms . In terms of the traditional parameter (ratio of thermal pressure to magnetic energy density; as magnetic field strengths increase), , so the rms driving-scale Mach number is . Magnetically-driven turbulence is close to purely solenoidal, so and .

We stress, though, that what is important is the saturated local plasma , which can be very different from the initial mean field threading the disk. As the MRI develops, the plasma field strength increases until it saturates in the fully nonlinear mode. Direct simulations have shown that for initial fields , (see e.g. Bai & Stone, 2012; Fromang et al., 2012, and references therein). The saturation occurs in rough equipartition with the thermal and kinetic energy densities – i.e. the turbulence is trans-sonic or even super-sonic (). In the weak-field limit, however, with , the saturation is much weaker, with , so .

These same simulations allow us to directly check our simple scaling with ; the authors directly measure the rms standard deviation in the (linear) density , which for a lognormal density distribution is (by our definitions) identical to . For four simulations with , , (but higher-resolution), , they see (midplane) saturation and (compared to a predicted , respectively). Moreover, these and a number of additional simulations have explicitly confirmed that our lognormal assumption (in the isothermal case) remains a good approximation for the shape of the density PDF (see Kowal et al., 2007; Lemaster & Stone, 2009; Kritsuk et al., 2011; Molina et al., 2012). So for a given and , our previous derivations remain valid.

#### 6.2.2 A Statistical Stability Criterion

The strong-field limit therefore leads to large density fluctuations. However, strong magnetic fields will also provide support against gravity, modifying the collapse criterion; this appears in Eq. B5. But for a given , this simply amounts (to lowest order) to the replacement . Because , near the Toomre scale, this is approximately equivalent to raising the stability parameter as (where is the including only thermal support). The energy and momentum of the bulk flows in the gas turbulence also provides support against collapse, so the “effective dispersion” in Eq. B5 includes all three effects; however this is already explicitly accounted for in our previous calculations for any . But while the effective increases in the strong-field limit with , so does , and the needed for statistical stability on long timescales (Fig. 2) increases exponentially with – so the net effect of MRI is always to increase the probability of stochastic collapse.

Putting this into our general criterion Eq. 9, we can write the statistical stability requirement

 Q(vA=0) ≡csκπGΣgas≈csΩπGΣgas ≥12√1+2β−1satexp[√2ln(τ0)ln(1+2/9β−1sat)] ≳12√1+2β−1satexp[(2/3)β−1/2sat√lnτ0] (21)

where as before and the latter uses the fact that is not extremely small in the cases of interest. Integrated over the lifetime of the disk, this becomes

 csΩπGΣgas ≳12√1+2β−1satexp(2.9β−1/2sat) (22)

This increases rapidly with increasing magnetic field strength: for .

So MRI with saturation (“seed” ) will make even disks statistically unstable, without the need for any other source of turbulence. On the other hand, weak-field MRI with produces only small corrections to statistical stability.

### 6.3. Convective Disks

A number of calculations have also shown that proto-planetary disks are convectively unstable over a range of radii (Boss, 2004; Boley et al., 2006; Mayer et al., 2007). Most simulations which see convection have also seen fragmentation, which has been interpreted as a consequence of convection enhancing the cooling rates until they satisfy the Gammie (2001) criterion for fragmentation. But Rafikov (2007) and others (Cai et al., 2006, 2008) have argued that while convection can and should develop in these circumstances, the radiative timescales at the photosphere push the cooling time above this threshold. However, as we have discussed above, that would not rule out rarer, stochastic direct collapse events.

Consider a polytropic thin disk; this is convectively unstable when it satisfies the Schwarzschild criterion

 dlnTdlnP>γ−1γ (23)

Following Lin & Papaloizou (1980), Bell & Lin (1994), and Rafikov (2007), using the fact that disk opacities can be approximated by ( the midplane pressure), this can also be writted . For the appropriate physical values, this implies strong convective instability in disks with (where is dominated by ice grains) and at higher temperatures K when grains sublimate, and marginal convective instability in between. So this should be a common process.

A convective disk can then accelerate gas via buoyancy at a rate comparable to the gravitational acceleration, implying mach numbers (depending on the driving gradients; recall also this is the three-dimensional ) at the scale height where the disk becomes optically thin.999From mixing-length theory, we can equate the convective energy flux at the scale-height (where at , the acceleration , , and for the relevant parameters in cgs) to the cooling flux . This gives us the approximate estimate . This agrees well with the simulations in the text when , but extrapolates to super-sonic values at low and/or large , so convection could be considerably more important than we estimate if it does not saturate at velocities . Buoyancy-driven turbulence is primarily solenoidal forcing, so while , leading to a “maximal” (assuming the convection cannot become supersonic; this is approximately what is measured in these simulations). If this saturation level is independent of (provided the disk is convectively unstable at all), we then simply need to examine Fig. 2 to determine for statistical stability; from Eq. 9 this is approximately

 Q≥0.5exp[Mc√2lnτ0]∼3 (24)

Thus, while this does not dramatically alter the behavior of the stability criterion , it does systematically increase the threshold for statistical stability by a non-trivial factor. And indeed, in the simulations of Mayer et al. (2007) and Boss (2004), fragmentation occurs when convection is present at radii where .

### 6.4. Additional Sources of Turbulence

There are many additional processes that may drive turbulence in proto-planetary and other Keplerian disks, but under most regimes they are less significant for our calculation here.

In the midplane of a protoplanetary disk, where large grains and boulders settle and are only weakly aerodynamically coupled to the gas, Kelvin-Helmholtz and streaming instabilities generate turbulence. However these only operate in a thin dust layer, and appear to drive rather small Mach numbers in the gas, so are unlikely to be relevant for direct collapse in the gas and we do not consider them further (see e.g. Bai & Stone, 2010; Shi & Chiang, 2012). They may, however, be critical for self-gravity of those grains participating in the instabilities themselves – a more detailed investigation of this possibility is outside the scope of this work (since our derivation does not apply to a weakly coupled, nearly-collisionless grain population), but extremely interesting for future study.

Radiative instabilities should also operate if the disk is supported by radiation pressure, and/or in the surface layer if a wind is being radiatively accelerated off the disk by central illumination. The former case is not expected in the physical disk parameter space we consider; but if it were so, convective, magneto-rotational, and photon-bubble instabilities are also likely to be present (Blaes & Socrates, 2001; Thompson, 2008), which will drive turbulence that saturates in equipartition between magnetic, radiation, and turbulent energy densities, i.e. produce the equivalent of throughout the disk (giving results broadly similar to the strong-field MRI case). In the wind case, the Mach numbers involved can be quite large (since material is accelerated to the escape velocity), but unless the surface layer includes a large fraction of the mass, it is unlikely to be important to the process of direct collapse.

In the case of an AGN accretion disk, local feedback from stars in the disk may also drive turbulence (as it does in galactic disks), and this can certainly be significant in the outer regions of the disk where star formation occurs (see Thompson et al., 2005). In that case the turbulence may even be super-sonic, in which case a more appropriate model is that developed in Paper I-Paper II. In the inner parts of the disk, though, where the turbulence is sub-sonic, we are not necessarily interested in rare single star-formation events.

## 7. Example: ProtoPlanetary Disk and Fragmentation Radii

We now apply the statistical stability criteria derived above to a specific model of a proto-planetary disk. This is highly simplified, but it allows us to estimate physically reasonable sound speeds, cooling times, and other parameters, so allows us to ask whether our revised statistical stability criteria are, in practice, important.

### 7.1. Disk Model Parameters

For convenience, consider a disk with a simple power-law surface density profile

 Σ=Σ0,10001000gcm−2(r∗/au)−α (25)

The minimum mass solar nebula (MMSN) corresponds to and , but we consider a range in these parameters below.

For a passive flared disk irradiated by a central star with radius and temperature , the effective temperature is (Chiang & Goldreich, 1997):

 Teff,∗≈(αTR2∗4r2∗)1/4T∗ (26)

where for a solar-type star , and . If, instead of irradiation, disk heating is dominated by energy from steady-state accretion with some , energy balance requires an effective temperature (i.e. disk flux)

 Teff,acc≈[38π˙MΩ2σSB]1/4 (27)

where is the Stefan-Boltzmann constant. The temperature of interest for our purposes, however, is the midplane temperature , since this is where the disk densities are largest and what provides the resisting collapse; this is related to by a function of opacity and which we detail in Appendix A. But having determined and , it is straightforward to calculate ; the sound speed is where is the Boltzmann constant and is the mean molecular weight.

This is sufficient to specify most of the parameters of interest. In the gravito-turbulent model, however, we also require an estimate of the cooling time to estimate . Rafikov (2007) calculate the approximate cooling time for a convective and radiative disk (depending on whether or not it is convective and, if so, accounting for the rate-limiting of cooling by the disk photosphere). This gives

 tcool =Σc2sσSBT4f(τ) (28) ≈2×104yr(Σ103gcm−2)(T100K)−3f(τ)103

where is a function (shown in Appendix A) of the opacity which interpolates between the optically thin/thick, and convective/radiative regimes.

With these parameters calculated, for a given assumption about what drives the turbulence – e.g. MRI, gravitoturbulence, convection, etc. – the compressive Mach number can be calculated following § 6. We also technically need to assume the details of the turbulent spectral shape, for which we will assume a spectral index and non-intermittent , as well as the gas equation of state, for which we take , appropriate for molecular hydrogen. But these choices have small effects on our results, as shown in Fig. 1.

### 7.2. The Characteristic Initial Fragment Mass

In Fig. 3, we use this model to calculate the expected mass of a self-gravitating “fragment.” Varying , , and , we calculate the expected and , assuming a constant accretion rate of