TURBULENT DISKS ARE NEVER STABLE: FRAGMENTATION AND TURBULENCE-PROMOTED PLANET FORMATION

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 protoplanetary and proto-stellar disks as well as active galactic nucleus (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 protoplanetary 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; D'Angelo 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 protoplanetary 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 protoplanetary disks. The classic Toomre Q-criterion,

$$\begin{equation} Q \equiv \frac{\sigma _{\rm g}\,\kappa }{\pi \,G\,\Sigma _{\rm gas}} > 1, \end{equation} \tag{ 1 }$$

where σ_g is the gas velocity dispersion, κ ∼ Ω the epicyclic frequency, and Σ_gas the gas surface density (Toomre 1964, 1977; Goldreich & Lynden-Bell 1965), appears to offer some guidance. And indeed, for Q ≪ 1, 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, and 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 Q ≳ 1 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 ∼50 times the dynamical time (see references above and, e.g., Rice et al. 2005, 2012; Meru & Bate 2011a, 2011b, 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

1.3. Paper Overview

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 Q < 1, 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 overdense 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 Q ≫ 1, and nearly all the disk mass will be stable, but rare density fluctuations might form an order-unity number of isolated self-gravitating "fragments."

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 M_*) with scale radius r_*, scale height h, mean surface density Σ and midplane density ρ₀ ≈ Σ/(2 h), and Toomre Q > 1. 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 ∼h; roughly speaking, a gas parcel of this size will collapse under self-gravity if it exceeds the Roche criterion (overcomes tidal forces): $\rho \gtrsim M_{\ast }/r_{\ast }^{3} \sim Q\,\rho _{0}$.⁵

So for Q ≳ 1 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 $\mathcal {M}= \langle v_{\rm turb}^{2} \rangle ^{1/2} / c_{s}$, the (volumetric) distribution of densities is approximately lognormal, with dispersion $\sigma _{\ln {\rho }} \sim \sqrt{\ln [1 + \mathcal {M}^{2}]}$ (${\approx} \mathcal {M}$ for $\mathcal {M}\lesssim 1$).

Since density fluctuations exceeding ρ ≳ Q ρ₀ (ln (ρ/ρ₀) ≳ ln Q) can collapse, we integrate the tail of the lognormal density probability distribution function (PDF) above this critical density to estimate the probability $P_{c}\sim {\rm erfc}[(\ln {Q})/(\sqrt{2}\,\sigma _{\ln {\rho }})]$, per unit volume, of a region being self-gravitating at a given instant. Since we are considering regions of size ∼h, the disk contains approximately ∼(r_*/h)² independent volumes, and so (assuming that P_c is small) the probability of any one such volume having a coherent volume-average density above the critical threshold is P_v ∼ (r_*/h)² P_c.

Now consider that a typical protoplanetary disk (at a few AU) might have Q ∼ 100 and $\mathcal {M}\sim 1$, with h/r_* ∼ 10⁻¹. In this case, P_c ∼ 3 × 10⁻⁸ and P_v ∼ 3 × 10⁻⁶ 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 ∼h/v_t (since this is the time for fluid flows to cross and interact with a new independent region of size ∼h). So the density PDF is re-sampled or "refreshed" on the timescale ${\sim } h/v_{t} \sim 1/(\mathcal {M}\,\Omega)$. But this is just a dynamical time, which is very short relative to a typical disk lifetime (Ω⁻¹ ∼ yr at ∼10 AU). If the disk survives for a total timescale t₀, then the entire volume is re-sampled ${\sim } \mathcal {M}\,t_{0}\,\Omega$ independent times. So the probability, integrated over time, of any one of these volumes, at any one time, exceeding the self-gravity criterion, is $P_{\rm frag}^{\rm int}\sim \mathcal {M}\,t_{0}\,\Omega \,(r_{\ast }/h)^{2}\,{\rm erfc}[\ln {Q}/(\sqrt{2}\,\sigma _{\ln {\rho }})]$.

If a typical disk with parameters above survives for ∼Myr, or ∼10⁶ crossing times, we then obtain a time-and-volume-integrated probability $P_{\rm frag}^{\rm int}\sim 1$ of at least a single stochastic "fragmentation event" driven by turbulent density fluctuations. The mass of the self-gravitating "fragment" will be ∼(4π/3) h³ ρ ∼ 4 Q h³ ρ₀ ∼4 (h/r_*)³ M_* (∼0.1–1 M_Jupiter, 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 ${\sim }1/(\mathcal {M}\,\Omega)\sim$ yr and has a short collapse/free-fall time ${\sim }1/\sqrt{G\rho } \sim \Omega ^{-1}\sim$ yr. Despite having Q ∼ 100, then, we estimate that this disk could be statistically unstable!

For otherwise fixed h/r_* ∼ 0.1 and $\mathcal {M}\sim 1$, $P_{\rm frag}^{\rm int}$ declines exponentially with increasing Q. So for Q ≫ 100, $P_{\rm frag}^{\rm int}\ll 1$ and such a disk is both classically and statistically stable. For Q < 1, 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 1 ≲ Q ≲ 100, 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.

As mentioned above, turbulence in approximately isothermal gas (neglecting self-gravity) generically drives the density PDF to an approximate lognormal (a normal distribution in ln (ρ)). 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 S simply scales with the driving-scale Mach number $\mathcal {M}$ as $S=\ln {[1+b^{2}\,\mathcal {M}^{2}]}$ (where b is a constant discussed below). These scalings have been confirmed by a huge number of numerical experiments with $\mathcal {M}\sim 0.01\hbox{--}100$, sampling the PDF down to values as low as ∼10⁻¹⁰ 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 b or the effectively compressible component of $\mathcal {M}$; 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).⁶ 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 $S(\mathcal {M})\rightarrow S(\mathcal {M}\,|\,\rho) = \ln {[1+b\,v_{t}^{2}/c_{s}^{2}(\rho)]}$, allowing an appropriately modified PDF that provides an excellent approximation to the results in simulations (see Scalo et al. 1998; Passot & Vazquez-Semadeni 1998).

In Paper I and 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 ρ_x(R) = ρ(x | R): the density field about the (random) coordinate x in space, smoothed with some window function of characteristic radius R (e.g., the average density in a sphere). If the gas is isothermal (we discuss more complicated cases below), this is distributed as a lognormal:

$$\begin{equation} P(\ln [\rho _{\bf x}(R)]) = \frac{1}{\sqrt{2\pi \,S(R)}}\,\exp {\left(-\frac{[\ln {\rho _{\bf x}(R)} + S(R)/2 ]^{2}}{2\,S(R)} \right)}, \end{equation} \tag{ 2 }$$

where S(R) is the variance on each scale R. S(R) is just the Fourier transform of the density power spectrum, which itself follows from the (well-defined) velocity power spectrum (see Equation (B2)). Essentially, S(R) is determined by integrating the contribution to the variance on each scale from the velocity field, using the relation ${\Delta }S\approx \ln {(1 + b^{2}\,\mathcal {M}(R)^{2})}$ ($\mathcal {M}(R)$ is the scale-dependent Mach number of the velocity field).

"Interesting" regions are those where ρ_x(R) 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

$$\begin{eqnarray} &&\frac{\rho _{\rm crit}(R)}{\rho _{0}} \equiv \frac{Q}{2\,\tilde{\kappa }}\,\left(1+\frac{h}{R} \right) \left [ \frac{\sigma _{g}^{2}(R,\,\rho _{\rm crit})}{\sigma _{g}^{2}(h,\,\rho _{0})}\,\frac{h}{R} + \tilde{\kappa }^{2}\,\frac{R}{h}\right ] \end{eqnarray} \tag{ 3 }$$

(Vandervoort 1970). Here ρ₀ is the mean midplane density of the disk, h is the disk scale height, $\tilde{\kappa }\equiv \kappa /\Omega =1$ for a Keplerian disk (κ is the epicylic frequency), and Q ≡ (σ_g[h, ρ₀] κ)/(π G Σ_gas) is the Toomre Q parameter. The effective gas dispersion $\sigma _{g}^{2}(R,\,\rho) = c_{s}^{2}(\rho) + \langle v_{t}^{2}(R) \rangle + v_{\rm A}^{2}(\rho,\,R)$ (Equation (B5)) includes the thermal (c_s), turbulent (v_t), and magnetic support (Alfvén speed v_A). A full derivation is given in Paper I, but we stress that this not only implies/requires that a region with ρ(R) > ρ_crit(R) 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).⁷ Knowing the size R and critical density ρ_crit of a region, it is trivial to translate this to the total collapsing gas mass M = M(R).

We desire the mass and initial size spectrum of regions which exceed ρ_crit (so are self-gravitating) on the scale R 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 that this reduces to a derivation of the "first-crossing" distribution for the field ρ_x(R). The mass function of collapsing objects can be written generally as

$$\begin{equation} \frac{{d}n}{{d}M} = \frac{\rho _{\rm crit}(M)}{M}\,f_{f}(M)\,\left |\frac{{d}S}{{d}M} \right |, \end{equation} \tag{ 4 }$$

where f_f(S) is a function given in Equation (B8) that is somewhat cumbersome to derive (see Paper I for details), but depends only on how the dimensionless quantities S(R) and ρ_crit(R)/ρ₀ "run" as a function of scale R.

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 (Section 7) and Paper II (Sections 9 and 10), but to very good approximation (for M dn/dln M ≲ ρ₀) this is just

$$\begin{equation} \Delta P(M,\,\Delta t)\approx \frac{{d}n}{{d}M}\,\frac{\Delta t}{\tau (M)}, \end{equation} \tag{ 5 }$$

where τ(M) = R[M]/v_t(R[M]) is the turbulent crossing time on the scale R[M] corresponding to the mass (Equation (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 M_d. By the definitions used in Paper I and Paper II, this means the total "effective volume" is M_d/ρ₀, i.e., that the total (integrated) number of objects per unit mass is

$$\begin{equation} \frac{{d}N}{{d}\log {M}} = \frac{M_{d}}{\rho _{0}}\,\frac{{d}n}{{d}\log {M}} = \frac{M_{d}}{M}\,\frac{\rho _{\rm crit}}{\rho _{0}}\,f_{f}(M)\,\left |\frac{{d}S}{{d}\log {M}} \right |. \end{equation} \tag{ 6 }$$

For any polytropic gas, we can factor out all dimensional parameters and work in units of h and ρ₀. Mass then has units of ρ₀ h³. But since for a disk in vertical equilibrium h = σ_g/Ω and for a Keplerian disk $\kappa = \Omega = (G\,M_{\ast }\,r_{\ast }^{-3})^{1/2}$ (where M_* is the central mass and r_* is the disk-centric radius), we have $Q =(\sigma _{g}[h]\,\kappa)/(\pi \,G\,\Sigma _{\rm gas}) = h\,\Omega ^{2}/(\pi \,G\,\Sigma _{\rm gas}) = (h/r_{\ast })\,M_{\ast }/(\pi \,\Sigma _{\rm gas}\,r_{\ast }^{2}$). So we can write (h/r_*) = Q μ, where

$$\begin{equation} \mu \equiv \frac{\pi \,\Sigma _{\rm gas}\,r_{\ast }^{2}}{M_{\ast }} \approx \frac{M_{d}}{M_{\ast }}, \end{equation} \tag{ 7 }$$

where M_d is the disk mass within R. For an exponential vertical profile used to define our dispersion relation, Σ_gas = 2 ρ₀ h, so we can use this and (h/r_*) = Q μ to link ρ₀ h³ = (2π)⁻¹ (μ Q)² M_d. We can then remove the local quantities ρ₀ and h and define mass in the much simpler global units of μ and M_d. The units of dN defined above are M_d/(ρ₀ h³), so this can be similarly re-written. And since the disk crossing time scales as ∼h/σ_g[h] = Ω⁻¹, Ω provides a natural time unit.

Having defined units, the model is completely specified by dimensionless parameters. These are the spectral index p of the turbulent velocity spectrum, E(k)∝k^−p ($v_{t}^{2}(R)\propto R^{p-1}$), and its normalization, which we define by the Mach number on large scales $\mathcal {M}_{h}^{2}\equiv \langle v_{t}^{2}(h)\rangle /c_{s}^{2}$, as well as the Toomre parameter Q. We must also specify the parameter b, i.e., the mean fraction of the velocity in compressive modes. This is b = 1 for purely compressively forced turbulence, b = 1/3 for purely solenoidally forced turbulence, and b = 1/2 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 number⁸ $\mathcal {M}_{c}\equiv b\,\mathcal {M}$.

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" f_f(M) 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 Equation (B8) above (see Paper II; Appendix D), leaving the rest of our derivation intact. We will consider such non-Gaussian, correlated statistics and show that they give very similar results. For non-isothermal cases, we must replace c_s → c_s(ρ) (in calculating both the variance and critical density) and again allow for the density PDF to be non-Gaussian (with a skew toward lower/higher densities as the equation of state is made more or less "stiff," respectively; see Sections 3 and 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 $c_{s}^{2} \propto \rho ^{\gamma -1}$ (γ is the polytropic index). Magnetic fields can produce global anisotropy, but this can be simply absorbed into the form of f_f(M) as well; their largest effect in suppressing fluctuations manifests as a field-strength-dependent value b (Paper II, Section 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) and so is within the range of b variations we will consider.⁹

In Figure 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.¹⁰

Zoom In Zoom Out Reset image size

We define our "reference" model to be isothermal (γ = 1), non-intermittent (i.e., lognormal density PDF), with turbulent spectral slope p = 2 (appropriate for compressible turbulence), Toomre Q = 1, and have b = 1/2 (random forcing) with rms three-dimensional compressive Mach number on the largest scales $\mathcal {M}_{c}=1$. 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 Figure 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 ∼0.1–1 μ² M_d. For a relatively massive disk M_d ∼ 0.1 M_*, then, we expect characteristic masses in such events of ∼10⁻³ M_*, corresponding to gas giants. If, however, the disk mass is lower, M_d ∼ 0.01 M_*, then this becomes ∼10⁻⁶ M_*, typical of rocky (Earth-mass) planets.

The dependence on the slope of the turbulent spectrum is quite weak: Kolmogorov-like (p = 5/3, appropriate for incompressible turbulence) spectra give nearly identical results; shallower spectra slightly broaden the mass range predicted (since $\mathcal {M}_{c}$ declines more slowly at small scales), but these are not seen in realistic astrophysical contexts. Likewise, at fixed $\mathcal {M}_{c}$, there is some dependence on b (because of how $\mathcal {M}= \mathcal {M}_{c}/b$ 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 lognormal 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 $\mathcal {M}_{c}$, 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 Q, 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 $\mathcal {M}_{c}$ and Q; at $\mathcal {M}_{c}\gtrsim 1$ or Q ≲ 1 the mass distribution rapidly becomes more broad, while at $\mathcal {M}_{c}\lesssim 1$ or Q ≳ 1 the characteristic mass remains fixed but the normalization of the probability becomes exponentially suppressed.

This is summarized in Figure 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,

$$\begin{equation} \frac{{d}N_{\rm frag}}{{d}t} = \int \,{d}\log {M}\,\frac{{d}N}{{d}\log {M}\,{d}t}, \end{equation} \tag{ 8 }$$

and plot this as a function of $\mathcal {M}_{c}$ for different parameter choices. As before, most parameters make a surprisingly small difference; $\mathcal {M}_{c}$ and Q dominate.

Zoom In Zoom Out Reset image size

5.1. A General Statistical Stability Criterion

Formally dN_frag/dt is always non-zero for $\mathcal {M}_{c}>0$, but we see that for $\mathcal {M}_{c}\lesssim 1/2$, the probability per unit time of forming a self-gravitating fluctuation drops rapidly. However, recall that the total lifetime of, e.g., a protoplanetary disk is many disk dynamical times ∼Ω⁻¹, and our "time unit" is μ² Ω⁻¹. Consider a typical lifetime of τ_Myr ≡ τ_disk/Myr ∼ 1; then the disk at ∼10 AU around a solar-mass star (where Ω ≈ 1 yr⁻¹) experiences ∼10⁶ dynamical times. For a disk-to-total mass ratio of μ ∼ 0.1, this is ∼10⁸ "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 dN_frag/dt ∼ 10⁻⁸ in the units here (i.e., $\mathcal {M}_{c}$ as small as ∼0.15). Of course, following such an integration in detail requires knowing the evolution of the disk mass, Q, $\mathcal {M}_{c}$, etc. But we can obtain some estimate of the value of Q required for statistical stability (ensuring that the probability of fragmentation events is negligible) by simply assuming that all quantities are constant and integrating over an approximate timescale, also shown in Figure 2. Here we take our standard model, consider three timescales in units of μ² Ω⁻¹ (a factor of ∼100 shorter and longer than the value motivated above), and consider the minimum Q at each $\mathcal {M}_{c}$ needed to ensure that the time-integrated probability of a fragmentation event is ≪1.

This minimum Q_min is ∼1 at $\mathcal {M}\sim 0.1$; equivalently, disks with Q ≈ 1 and $\mathcal {M}_{c}\gtrsim 0.1$ have an order-unity probability of at least one stochastic fragmentation event over their lifetime ($P_{\rm frag}^{\rm int}\sim 1$). At larger $\mathcal {M}_{c}\gtrsim 0.3\hbox{--}0.5$, Q ≳ 3–5 is required for $P_{\rm frag}^{\rm int}\ll 1$; and by sonic Mach numbers $\mathcal {M}_{c}\sim 1\hbox{--}3$, Q ≳ 40–1000 is required for $P_{\rm frag}^{\rm int}\ll 1$.

We can approximate the scaling of $Q_{\rm min}(\mathcal {M}_{c})$ by the following: recall that the critical density near the Toomre scale h scales approximately as ln (ρ_crit/ρ₀) ∼ ln (2 Q) (Equation (3)), while the density dispersion scales as $\sigma _{\ln {\rho }}\sim \sqrt{(1+\mathcal {M}_{c}^{2})}$ (Equation (B2)). As noted above, if the system evolves for a total timescale τ₀ = t₀/(μ² Ω⁻¹) (time in our dimensionless units), then an event with probability per unit time P ≈ 1/τ₀ has an order-unity probability of occurring. If the probabilities are approximately normally distributed, then this is just exp (− B²/(2 S)) ≈ 1/τ₀, where B is the barrier and S the variance. Since the mass function is peaked near the Toomre scale, we can approximate both by their values near the "driving scale" ∼h, B ≈ ln (2 Q) and $S\approx \ln {(1+\mathcal {M}_{c}^{2}})$. Thus, statistical stability over some timescale of interest t₀ requires a Q_min in Figure 2 of

$$\begin{eqnarray} Q_{\rm min} &\approx 0.5\,\exp \left[\sqrt{2\,\ln {(t_{0}\,\mu ^{-2}\,\Omega)}\,\ln {\left(1+\mathcal {M}_{c}^{2} \right)}} \right]. \end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray} &&Q_{\rm min}\sim 0.5\,\exp \left(6\,\sqrt{\ln {\left(1+\mathcal {M}_{c}^{2} \right)}}\right). \end{eqnarray} \tag{ 10 }$$

In other words, a ∼5–6 σ_ln ρ event has order-unity chance of occurring once over the disk lifetime, so for any $\mathcal {M}_{c}$ this implies a minimum Q needed to ensure statistical stability in such an extreme event.

Thus far, we have considered the general case, varying the Mach numbers $\mathcal {M}_{c}$ independent of other disk properties such as Q, Σ, 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 t_cool = ζ Ω⁻¹ locally fixed). In particular, this includes the scenario of a "gravito-turbulent" steady-state from Gammie (2001), with local instabilities (spiral waves) powering turbulence that contributes an effective viscosity and maintains a steady temperature and Q ∼ 1. 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 Q, $\mathcal {M}_{c}$, γ, 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,

$$\begin{equation} t_{\rm cool}=\zeta \,\Omega ^{-1}. \end{equation} \tag{ 11 }$$

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

$$\begin{equation} \alpha \approx \frac{4}{9}\,\frac{1}{\gamma \,(\gamma -1)\,\zeta }, \end{equation} \tag{ 12 }$$

is approximately constant. This corresponds to the amplitude of density waves δΣ/Σ∝α^1/2, leading to a maximum α_max ≈ 0.06 above (hence minimum ζ_min below) which the local Q < 1 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

$$\begin{equation} \dot{M}=3\pi \,\alpha \,c_{s}^{2}\,\Sigma _{\rm gas}\,\Omega ^{-1}, \end{equation} \tag{ 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, ${d}E/{d}A\,{d}t = (p-1)^{-1}\,\Sigma _{\rm gas}\,v_{t}^{2}\,\Omega$ (where here we consider the rate per unit area) and equate it to the cooling rate $=[\gamma \,(\gamma -1)]^{-1}\,\Sigma _{\rm gas}\,c_{s}^{2}\,t_{\rm cool}^{-1}$, giving $\mathcal {M}_{c}^{2} \approx (3/2)\,\alpha$. Equivalently, we could have equated the Reynolds stress that leads to α, $\alpha =({d}\ln {\Omega }/{d}\ln {r_{\ast }})^{-1}\,T_{{\rm Rey}}/(\Sigma _{\rm gas}\,c_{s}^{2})$ with T_Rey = 〈Σ_gas δv_rδv_ϕ〉; we obtain

$$\begin{equation} \mathcal {M}_{c}^{2} \approx \frac{3}{2}\,\alpha \approx \frac{2}{3}\,\frac{1}{\gamma \,(\gamma -1)\,\zeta } \end{equation} \tag{ 14 }$$

again.¹¹ 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 Q ∼ 1 needed to power the turbulence (since spiral waves are still generated) requires a smaller turbulent dispersion and 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 Q < 1 but also for stochastic local turbulent density fluctuations. For a given Q, we have calculated (Figure 2) the Mach number $\mathcal {M}_{c}$ above which the system will be probabilistically likely to fragment on a given timescale. Equation (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 Section 8). Recall that for a system which evolves for a total timescale τ₀ = t₀/(μ² Ω⁻¹) (time in the dimensionless units we adopt), we obtain the approximate Equation (9) for the Q_min needed to ensure $P_{\rm frag}^{\rm int}\ll 1$:

$$\begin{equation} Q \ge \frac{1}{2}\,\exp {\left[\sqrt{2\,\ln {(\tau _{0})}\,\ln {\left(1+\mathcal {M}_{c}^{2} \right)}} \right]}. \end{equation} \tag{ 15 }$$

We can invert this to find the maximum $\mathcal {M}_{c}$ for statistical stability at a given Q:

$$\begin{equation} \mathcal {M}_{c}^{2} \lesssim \exp \left [ \frac{[\ln {(2\,Q)}]^{2}}{2\,\ln {\tau _{0}}} \right ] - 1 \approx \frac{[\ln {(2\,Q)}]^{2}}{2\,\ln {\tau _{0}}}, \end{equation} \tag{ 16 }$$

where the second equality follows from the fact that (for the systems of interest) 2 ln (τ₀) ≫ ln (2 Q) is almost always true.

Combining this with Equation (14), we obtain

$$\begin{eqnarray} \zeta _{\rm min} &\approx & \frac{4}{3\,\gamma \,(\gamma -1)}\,\frac{\ln {(\tau _{0})}}{[\ln {(2\,Q)}]^{2}}\\ &=& \frac{4}{3\,\gamma \,(\gamma -1)}\,\frac{\ln {[t_{0}/(\mu ^{2}\,\Omega ^{-1})]}}{[\ln {(2\,Q)}]^{2}}. \end{eqnarray} \tag{ 17 }$$

We can immediately see some important consequences. Because of the stochastic nature of turbulent density fluctuations, ζ_min will never converge in time integration (assuming that 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 log (Q), so indeed higher-Q disks are "more stable," but there is no "hard" cutoff at a specific Q 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

$$\begin{equation} \langle t(N_{\rm frag}\sim 1) \rangle \approx \mu ^{2}\,\Omega ^{-1}\,\exp \left (\frac{3}{4}\,\zeta \,\gamma \,(\gamma -1)\,[\ln {(2\,Q)}]^{2} \right). \end{equation} \tag{ 19 }$$

As expected, this quickly becomes large for modest ζ and/or Q. 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 ζ_min, what do we expect in realistic systems such as protoplanetary disks? For a physical disk with ζ ≫ 1 we should expect γ ≈ 7/5–5/3 (depending on the gas phase), and in equilibrium Q ∼ 1; and we should integrate over the entire lifetime of the disk, τ₀ = τ_tot ∼ 10⁶–10¹⁰ (with values motivated in Section 5). We then expect

$$\begin{equation} \zeta _{\rm min}^{\rm int} \approx \frac{51}{\gamma \,(\gamma -1)}\,\frac{1+0.05\,\ln {(\tau _{0}^{\rm int}/10^{8})}}{(1+1.4\,\ln {Q})^{2}}. \end{equation} \tag{ 20 }$$

This is fairly sensitive to Q—note $\zeta _{\rm min}^{\rm int}\approx 13/[\gamma (\gamma -1)]$ if Q = 2 instead—and weakly sensitive to τ₀ for reasonable variations. But this implies that only disks with extremely slow cooling, t_cool ≳ 50 Ω⁻¹ (corresponding to steady-state Mach numbers $\mathcal {M}_{c}\lesssim 0.02$), are truly statistically stable with Q ≈ 1 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 γ = 2 and a steady-state Q ≈ 2.46, evolving their simulations for typical timescales t₀ ∼ 50 Ω⁻¹ (though they consider some longer-scale runs). Because these were two-dimensional shearing-sheet simulations, the appropriate μ is somewhat ambiguous, but recall (h/r_*) = Q μ by our definitions, and for the assumptions in Gammie (2001) their "standard" simulation corresponds to an h/r_* ≈ 0.01 Q (where we equate the "full disk size" to the area of the box simulated). Plugging in these values, then, we predict ζ_min = 3.4, in excellent agreement with the value ζ ≈ 3 found by trial of several values therein. Paardekooper (2012) considered very similar simulations but with Q ≈ 1 and all sheets run for t₀ ∼ 1000 Ω⁻¹; for this system we predict ζ_min = 20.4—again almost exactly their estimated "fragmentation boundary." Meru & Bate (2012) and Rice et al. (2012) consider three-dimensional global simulations; here μ = M_d/M_* = 0.1 is well defined, a more realistic γ = 5/3 is adopted, and the disks self-regulate at Q ≈ 1. The simulations are run for a shorter time ∼50–100 〈Ω〉⁻¹ (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 ζ_min ≈ 21–25, in very good agreement with where both simulations appear to converge (using either smoothed particle hydrodynamics or grid-based methods). This is also in good agreement with the earlier simulations in Rice et al. (2005), for γ = 5/3 (predicting ζ_min = 7.5, versus their estimated 6–7) and γ = 7/5 (predicting ζ_min = 14, versus their estimated 13). 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 Figure 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 ∼0.01–0.1 μ² M_d (depending on the exact parameters). This is equivalent to a spatial resolution of _r ∼ 0.02–0.2 Q^1/2 h, i.e., a small fraction of the disk scale height h. 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 Q ≫ 1. We therefore next consider the simple case where there is no gravo-turbulent instability (t_cool → ∞), 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 $\langle v_{t}^{2} \rangle ^{1/2} \approx v_{\rm A}$. In terms of the traditional β parameter (ratio of thermal pressure to magnetic energy density; β → 0 as magnetic field strengths increase), $\beta = 2\,c_{s}^{2}/v_{\rm A}^{2}$, so the rms driving-scale Mach number is $\mathcal {M}\approx \sqrt{2\,\beta ^{-1}}$. Magnetically driven turbulence is close to purely solenoidal, so b ≈ 1/3 and $\mathcal {M}_{c}=b\,\mathcal {M}\approx (\sqrt{2}/3)\,\beta ^{-1/2}$.

We stress, though, that what is important is the saturated local plasma β = β_sat, 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 β₀ ≲ 10⁴, β_sat ∼ 1/3–2/3 (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 transsonic or even super-sonic ($\mathcal {M}_{c}\sim \sqrt{2/3} \sim 0.8$). In the weak-field limit, however, with β₀ ≳ 10⁴, the saturation is much weaker, with β_sat ∼ 10–20, so $\mathcal {M}_{c}\sim 0.1$.

These same simulations allow us to directly check our simple scaling with v_A; the authors directly measure the rms standard deviation in the (linear) density δ = ρ/〈ρ〉, which for a lognormal density distribution is (by our definitions) identical to $\mathcal {M}_{c}$. For four simulations with β₀ = 10², 10³, 10³ (but higher resolution), 10⁴, they see (midplane) saturation β_sat = 0.4, 1.1, 0.7, 18 and 〈ρ²〉^1/2/〈ρ〉 = 0.60, 0.43, 0.53, 0.13 (compared to a predicted $\langle \rho ^{2}\rangle ^{1/2}/\langle \rho \rangle = \mathcal {M}_{c}= (\sqrt{2}/3)\,\beta _{\rm sat}^{-1/2} = 0.72,\,0.44,\,0.55,\,0.11$, 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 $\mathcal {M}_{c}$ and Q, 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 Equation (B5). But for a given β, this simply amounts (to lowest order) to the replacement $c_{s}^{2}\rightarrow c_{s}^{2} + v_{A}^{2} = c_{s}^{2}\,(1 + 2\,\beta ^{-1})$. Because Q∝c_s, near the Toomre scale, this is approximately equivalent to raising the stability parameter as $Q\rightarrow Q_{\rm eff} \equiv Q(v_{\rm A}=0)\,\sqrt{1 + 2\,\beta ^{-1}}$ (where Q(v_A = 0) is the Q including only thermal support). The energy and momentum of the bulk flows in the gas turbulence also provide support against collapse, so the "effective dispersion" in Equation (B5) includes all three effects; however, this is already explicitly accounted for in our previous calculations for any $\mathcal {M}$. But while the effective Q_eff increases in the strong-field limit with β^−1/2, so does $\mathcal {M}_{c}$, and the Q needed for statistical stability on long timescales (Figure 2) increases exponentially with $\mathcal {M}_{c}$—so the net effect of MRI is always to increase the probability of stochastic collapse.

Putting this into our general criterion Equation (9), we can write the statistical stability requirement

$$\begin{eqnarray} \nonumber Q(v_{\rm A}=0) &\equiv & \frac{c_{s}\,\kappa }{\pi \,G\,\Sigma _{\rm gas}} \approx \frac{c_{s}\,\Omega }{\pi \,G\,\Sigma _{\rm gas}} \\ \nonumber & \ge & \frac{1}{2\,\sqrt{1+2\,\beta _{\rm sat}^{-1}}}\,\exp \left [\sqrt{2\,\ln {(\tau _{0})}\,\ln {\left(1+2/9\,\beta _{\rm sat}^{-1} \right)}}\right] \\ &\gtrsim & \frac{1}{2\,\sqrt{1+2\,\beta _{\rm sat}^{-1}}}\,\exp \left[(2/3)\,\beta _{\rm sat}^{-1/2}\,\sqrt{\ln {\tau _{0}}} \right], \end{eqnarray} \tag{ 21 }$$

where τ₀ ≡ t₀ μ⁻² Ω 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

$$\begin{eqnarray} \frac{c_{s}\,\Omega }{\pi \,G\,\Sigma _{\rm gas}} & \gtrsim \frac{1}{2\,\sqrt{1+2\,\beta _{\rm sat}^{-1}}}\,\exp {\left(2.9\,\beta _{\rm sat}^{-1/2} \right)}. \end{eqnarray} \tag{ 22 }$$

This increases rapidly with increasing magnetic field strength: Q_min(v_A = 0) ≈ 15, 7, 3, 1.2 for β_sat = 1/3, 2/3, 2, 10.

So MRI with saturation β_sat ≲ 10 ("seed" β₀ ≲ 10⁴) will make even Q > 1 disks statistically unstable, without the need for any other source of turbulence. On the other hand, weak-field MRI with β_sat ≳ 10 produces only small corrections to statistical stability.

6.3. Convective Disks

A number of calculations have also shown that protoplanetary 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

$$\begin{equation} \frac{{d}\ln {T}}{{d}\ln {P}} > \frac{\gamma -1}{\gamma }. \end{equation} \tag{ 23 }$$

Following Lin & Papaloizou (1980), Bell & Lin (1994), and Rafikov (2007), using the fact that disk opacities can be approximated by κ ≈ κ₀ P^α T^β (P the midplane pressure), this can also be written as (1 + α)/(4 − β) > (γ − 1)/γ. For the appropriate physical values, this implies strong convective instability in disks with T ≲ 150 K (where κ is dominated by ice grains) and at higher temperatures ≳ 1.5 × 10³ 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 $\mathcal {M}^{2}\approx 0.25\hbox{--}1$ (depending on the driving gradients; recall also that this is the three-dimensional $\mathcal {M}$) at the scale height where the disk becomes optically thin.¹² Buoyancy-driven turbulence is primarily solenoidal forcing, so b ≈ 1/3 while $\mathcal {M}\sim 0.5\hbox{--}1$, leading to a "maximal" $\mathcal {M}_{c}\sim 0.2\hbox{--}0.3$ (assuming that the convection cannot become super-sonic; this is approximately what is measured in these simulations). If this saturation level is independent of Q (provided that the disk is convectively unstable at all), we then simply need to examine Figure 2 to determine Q_min for statistical stability; from Equation (9) this is approximately

$$\begin{equation} Q \ge 0.5\,\exp {\left [\mathcal {M}_{c}\,\sqrt{2\ln {\tau _{0}}}\right]} \sim 3. \end{equation} \tag{ 24 }$$

Thus, while this does not dramatically alter the behavior of the stability criterion Q, it does systematically increase the threshold Q 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 Q ≈ 1.4–1.8 > 1.

6.4. Additional Sources of Turbulence

We now apply the statistical stability criteria derived above to a specific model of a protoplanetary disk. This is highly simplified, but it allows us to estimate physically reasonable sound speeds, cooling times, and other parameters and 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

$$\begin{equation} \Sigma = \Sigma _{0,\,1000}\,{\rm 1000\,{\rm g\,cm^{-2}}}\,(r_{\ast }/{\rm AU})^{-\alpha }. \end{equation} \tag{ 25 }$$

The minimum-mass solar nebula (MMSN) corresponds to Σ_0, 1000 ∼ 1 and α ≈ 1.5, but we consider a range in these parameters below.

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

$$\begin{equation} T_{\rm eff,\,\ast } \approx \left (\frac{\alpha _{T}\,R_{\ast }^{2}}{4\,r_{\ast }^{2}} \right)^{1/4}\,T_{\ast }, \end{equation} \tag{ 26 }$$

where for a solar-type star α_T  ≈  0.005 (r_*/AU)⁻¹  + 0.05 (r_*/AU)^2/7, R_* = R_☉, and T_* = 6000 K. If, instead of irradiation, disk heating is dominated by energy from steady-state accretion with some $\dot{M}$, energy balance requires an effective temperature (i.e., disk flux)

$$\begin{equation} T_{\rm eff,\,acc} \approx \left [\frac{3}{8\pi }\,\frac{\dot{M}\,\Omega ^{2}}{\sigma _{\rm SB}} \right ]^{1/4}, \end{equation} \tag{ 27 }$$

where σ_SB is the Stefan–Boltzmann constant. The temperature of interest for our purposes, however, is the midplane temperature T_mid, since this is where the disk densities are largest and what provides the c_s resisting collapse; this is related to T_eff by a function of opacity and Σ, which we detail in Appendix A. But having determined T_eff and Σ, it is straightforward to calculate T_mid; the sound speed is $c_{s} = \sqrt{k_{B}\,T_{\rm mid}/\mu }$, where k_B 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 ζ ≡ t_coolΩ. Rafikov (2007) calculates 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

$$\begin{eqnarray} t_{\rm cool} &=& \frac{\Sigma \,c_{s}^{2}}{\sigma _{\rm SB}\,T^{4}}\,f(\tau) \nonumber\\ &\approx & 2\times 10^{4}\,{\rm yr}\,\left (\frac{\Sigma }{10^{3}\,{\rm g\,cm^{-2}}}\right)\,\left (\frac{T}{100\,{\rm K}} \right)^{-3}\,\frac{f(\tau)}{10^{3}},\quad \end{eqnarray} \tag{ 28 }$$

where f(τ) 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—the compressive Mach number $\mathcal {M}_{c}$ can be calculated following Section 6. We also technically need to assume the details of the turbulent spectral shape, for which we will assume a spectral index p = 2 and non-intermittent T = 0, as well as the gas equation of state, for which we take γ = 7/5, appropriate for molecular hydrogen. But these choices have small effects on our results, as shown in Figure 1.

7.2. The Characteristic Initial Fragment Mass

In Figure 3, we use this model to calculate the expected mass of a self-gravitating "fragment." Varying Σ_0, 1000, α, and M_*, we calculate the expected T and Q, assuming a constant accretion rate of $\dot{M} = 3\times 10^{-7}\,M_{\odot }\,{\rm yr^{-1}}$ at all radii,¹³ which dominates the disk temperature inside ∼10 AU. Given this, we calculate the mass spectrum as Figure 1 and define an average 〈M_collapse〉 (the mass-weighted, spectrum-integrated mass). Technically this depends on $\mathcal {M}_{c}$ and hence the turbulent driving mechanism, but the dependence is weak so we just assume $\mathcal {M}_{c}=0.1$ in all cases.

Zoom In Zoom Out Reset image size

In each case, 〈M_collapse〉 ∼ μ² Q² M_disk(< r_*) as expected. Since $\mu \equiv (\pi \,\Sigma \,r_{\ast }^{2})/M_{\ast }$, this increases with disk surface density or mass, and also with increasing disk-to-stellar mass ratio. Recall that $T_{\rm eff,\,acc}\propto (\dot{M}\,\Omega ^{2})^{1/4} \propto r_{\ast }^{-3/4}$ and $T_{\rm eff,\,\ast }\propto r_{\ast }^{-1/2}$. So modulo opacity corrections we expect $\langle M_{\rm collapse} \rangle \propto r_{\ast }^{3\,(1.5-\alpha)}$ at small radii ≲ a few AU, weakly increasing with r_*, and $\langle M_{\rm collapse} \rangle \propto r_{\ast }^{0.5+3\,(1.5-\alpha)}$ at large r_*, increasing more rapidly.

This is only the initial self-gravitating, bound mass—it may easily evolve in time, as discussed in Section 8. However, it is interesting that there is a broad range of masses possible, with Earth and super-Earth-like masses more common at ≲ 1 AU and giant planet masses more common at ≳ 10 AU.

7.3. Disk Temperatures at Which Direct Collapse Occurs

Given a mass profile, then for a source of turbulence in Section 6 we can translate the criteria for statistical stability—a probability $P_{\rm frag}^{\rm int}\ll 1$ of forming a fragment in a characteristic timescale ∼Myr—into a range of midplane temperatures T_mid.

Figure 4 shows this for a disk with Σ_0, 1000 = 1 and α = 1, around a solar-type star, as well as a disk with Σ_0, 1000 = 10. Choosing α = 1.5 gives a similar result but with the curve slopes systematically shifted. Together this spans a range in disk mass at ∼10 AU of ∼0.004–0.07 M_*. We compare the expected T_mid(r_*, Σ), for accretion rates $\dot{M}=0$ and $\dot{M}=3\times 10^{-7}\,M_{\odot }\,{\rm yr^{-1}}$.

Zoom In Zoom Out Reset image size

There is some T(r_*) below which Q < 1, so the disk will catastrophically fragment. But this is quite restrictive: Q < 1 only when Σ and r_* are large (for Σ_0, 1000 ∼ 10–100, r_* ≳ 10–100 AU for α ∼ 1–1.5, respectively; roughly where M_disk(< r_*) ≳ 0.1 M_*).

If the disk is convectively unstable, the resulting Mach numbers lead to a large temperature range over which Q > 1, so the disk is classically stable, but $P_{\rm frag}^{\rm int}\sim 1$; this is less likely to be relevant for an MMSN (r_* ≳ 50 AU for Σ_0, 1000 = 1) but can be sufficient at r_* ≳ 2–10 AU for Σ_0, 1000 ≳ 10. As discussed in Section 6.3, the convective Mach numbers are somewhat uncertain, so we show the calculation for the range therein.

Strong-field MRI—if/where it is active—produces even larger $\mathcal {M}$; this can greatly expand the range where $P_{\rm frag}^{\rm int}\sim 1$ even in an MMSN. Again there is a range of possible $\mathcal {M}_{c}$ in the saturated state, shown here. In the weak-field regime, $\mathcal {M}_{c}$ is smaller, giving results very similar to the convection prediction.

Gravitoturbulence (again, if/where it is active), interestingly, has the opposite dependence on temperature. Because cooling rates grow rapidly at higher T_mid, the expected $\mathcal {M}_{c}$ is larger and $P_{\rm frag}^{\rm int}$ is larger at higher T (despite higher Q). It is also much less sensitive to the disk surface density. If the process operates, this is sufficient to produce $P_{\rm frag}^{\rm int}\sim 1$ at r_* ≳ 2–5 AU, regardless of $\dot{M}$ for nearly all reasonable temperatures, and even at r_* ≲ 1 AU if there is a modest $\dot{M}$ to raise the temperature (hence cooling rate) to ≳ 300–1000 K. For the disk surface densities here, $P_{\rm frag}^{\rm int}\ll 1$ at r_* ≲ 1 AU requires $\dot{M} \lesssim 3\times 10^{-7}\,M_{\odot }\,{\rm yr^{-1}}$. However, we caution that at sufficiently high T and Q, the mechanism may not operate at all.

Beyond a certain radius (which depends on which combination of these mechanisms are active), perhaps the most important result is that there may be no temperature where $P_{\rm frag}^{\rm int}\ll 1$, for a given surface density.

7.4. Surface Densities at Which Direct Collapse Occurs

In Figure 5, we calculate the surface densities (as a function of radius r_* around a solar-type star) where disks are statistically unstable ($P_{\rm frag}^{\rm int}\sim 1$) if/when different turbulent driving mechanisms are active, as in Figure 4. Whereas in Figure 4 we allowed the temperature to be free, here we assume that it follows our best estimate T_mid(r_*, Σ) for either an accretion rate $\dot{M}=0$ or $\dot{M}=3\times 10^{-7}\,M_{\odot }\,{\rm yr^{-1}}$, but freely vary Σ(r_*).

Zoom In Zoom Out Reset image size

Roughly speaking, convection and MRI systematically lower the density at all radii where $P_{\rm frag}^{\rm int}\sim 1$. Classical instability (Q < 1) requires Σ ≳ 30 Σ_MMSN. If the disk is convective, it can have Q > 1 but be statistically unstable (vulnerable to direct fragmentation via turbulence density fluctuations) for Σ ≳ 10 Σ_MMSN at r_* ≳ 10 AU; at smaller radii the threshold is sensitive to accretion heating raising Q. If strong-field MRI is active, the threshold surface density for statistical instability is much lower: for small accretion rates and/or large radii ≳ 2–5 AU, even the MMSN can have $P_{\rm frag}^{\rm int}\sim 1$. Gravitoturbulence, if active, generates sufficient turbulence for statistical instability ($P_{\rm frag}^{\rm int}\sim 1$) even at Σ ∼ 0.1 Σ_MMSN, at radii r_* ≳ 1 AU regardless of accretion rate, and r_* ≲ 1 AU for $\dot{M}\gtrsim 3\times 10^{-7}\,M_{\odot }\,{\rm yr^{-1}}$.

Traditionally, a disk with Toomre Q > 1 is classically stable against gravitational collapse on all scales (modulo certain global gravitational instabilities). However, we show here that this is no longer true in a turbulent disk. Random turbulent density fluctuations can produce locally self-gravitating regions that will then collapse, even in Q > 1 disks.¹⁴ Formally, the probability of such an event is always non-zero, so strictly speaking turbulent disks are never "completely" stable, but can only be so statistically, if the probability of forming a self-gravitating fluctuation is small over the timescale of interest. Moreover, we can analytically predict the probability, as a function of total self-gravitating mass, of the formation of such a region per unit time in a disk (or disk annulus) with given properties.

We do this using the excursion-set formalism developed in Paper I and Paper II, which allows us to use the power spectra of turbulence to predict the statistical properties of turbulent density fluctuations. In previous papers, this has been applied to the structure of the ISM in galactic or molecular cloud disks; however, we show that it is straightforward to extend to a Keplerian, protoplanetary disk. The most important difference between the case here and a galactic disk is that in the galactic case, turbulence is highly super-sonic ($\mathcal {M}\sim 10\hbox{--}100$), as opposed to sub/transsonic. And in galactic disks, cooling is rapid and the disk is globally self-gravitating (non-Keplerian), so systems almost always converge rapidly to Q ≈ 1 (see Hopkins et al. 2012); here, we expect a wider range of Q. And finally, protoplanetary disks are very long-lived relative to their local-dynamical times, so even quite rare events (with a probability of, say, 10⁻⁶ per dynamical time) may be expected over the disk lifetime.

At Q ≈ 1, disks with $\mathcal {M}_{c}\gtrsim 0.1$ are classically stable but statistically unstable: they are likely, over the lifetime of the disk, to experience at least an order-unity number of "fragmentation" events (formation of self-gravitating, collapsing masses). As expected, higher-Q disks are "more stable" (although again we stress that this is only a probabilistic statement); for Q ∼ 3–5, values of $\mathcal {M}_{c}\gtrsim 0.3\hbox{--}0.5$ are required for statistical instability. If the turbulence is transsonic ($\mathcal {M}_{c}\sim 1\hbox{--}3$), values as large as Q ∼ 40–1000 can be statistically unstable! Generally speaking, we show that statistical stability (i.e., ensuring that the probability of a stochastic direct collapse event is ≪1) integrated over some timescale of interest t₀ requires

$$\begin{eqnarray*} Q_{\rm min} &\approx & 0.5\,\exp \left [\sqrt{2\,\ln {(t_{0}\,\mu ^{-2}\,\Omega)}\,\ln {(1+\mathcal {M}_{c}^{2})}} \right] \\ &\sim & 0.5\,\exp \left (6\,\sqrt{\ln {\left(1+\mathcal {M}_{c}^{2}\right)}}\right) \end{eqnarray*}$$

(Figure 2 and Equation (9)). At the Mach numbers of interest, this is exponentially increasing with $\mathcal {M}$!

This is a radical revision to traditional stability criteria. However, the traditional Toomre Q criterion is not irrelevant. It is a necessary, but not sufficient, criterion for statistical stability, which should not be surprising since its derivation assumes a homogenous, non-turbulent disk. If Q ≪ 1, then "catastrophic" fragmentation occurs even for $\mathcal {M}_{c}\rightarrow 0$; all mass in the disk is (classically) unstable to self-gravity, and collapse proceeds on a single free-fall time. If Q > 1, fragmentation transitions to the stochastic (and slower) statistical mode calculated here, dependent on random turbulent density fluctuations forming locally self-gravitating regions.

Likewise, the criterion in Gammie (2001), that the cooling time be longer than a couple times the dynamical time, is not sufficient for statistical stability. For a given turbulent Mach number and Q, we show that assuming a stiffer equation of state has quite weak effects on the "stochastic" mode of fragmentation; in fact, even pure adiabatic gas (no cooling) produces very similar statistics if it can sustain a similar $\mathcal {M}_{c}$. Again, the key is that the Gammie (2001) criterion is really about the prevention of catastrophic fragmentation; as noted therein, a sufficiently slow cooling time allows turbulence to maintain a steady-state Q ∼ 1, and dissipation of that turbulence (driven by gravitational density waves and inflow) can maintain the gas thermal energy (c_s). Thus, faster cooling leads to catastrophic fragmentation of most of the mass on a single dynamical time (Q < 1 and $\mathcal {M}_{c}\gg 1$). And indeed this is the case from Paper I in a galactic disk, where t_cool ≪ t_dyn and the mass is only "recycled" back into the diffuse medium by additional energy input (from stellar feedback).

However, provided that the Gammie (2001) criterion is met and Q > 1 everywhere, we still predict fragmentation in the "stochastic" mode if gravitoturbulence operates. The cooling time is then important insofar as it changes the equilibrium balance between turbulent and thermal energy, i.e., appears in Q and governs $\mathcal {M}_{c}\propto (t_{\rm cool}\,\Omega)^{-1/2}$. This, indeed, has now been seen in a growing number of simulations (see references in Section 1), involving either larger volumes and/or longer runtimes. We consider the application of our theory to these specific models in Section 6.1; this allows us to predict a revised cooling-time criterion in this "mode," required for statistical stability over any timescale of interest:

$$\begin{eqnarray} &&\zeta \equiv \frac{t_{\rm cool}}{\Omega ^{-1}} > \frac{4}{3\,\gamma \,(\gamma -1)}\,\frac{\ln {[t_{0}\,\mu ^{-2}\,\Omega)]}}{[\ln {(2\,Q)}]^{2}}. \end{eqnarray} \tag{ 29 }$$

This provides an excellent explanation for the results in these simulations and resolves the apparent discrepancies between them noted in Section 1.

If another process is able to drive turbulence, then stochastic direct collapse might occur with even longer cooling times. We show that if the MRI is active and saturates at strong-field β_sat ∼ 1, the required Q for complete suppression of fragmentation can be very large (≳ 10–15; scaling with β_sat as Equation (21)), independent of the cooling time. If the MRI is not active (in the dead zone, for example), or if it saturates at weak-field values β_sat ≳ 10, and cooling is slower than the limit above, then convection may be the dominant driver of turbulence. This is sufficient to produce stochastic fragmentation events in the range Q ∼ 1–3, though probably not much larger.

We apply these calculations to specific models of protoplanetary disks that attempt to self-consistently calculate their temperatures and cooling rates. Doing so, we show that the parameter space where statistical instability and stochastic fragmentation may occur is far larger than that of classical instability (where Q < 1) and can include most of the disk even in an MMSN. Gravitoturbulence appears to be the most important channel driving stochastic fragmentation when it is active, and is sufficient to produce an order-unity number of events in disks with Σ ≳ 0.1 Σ_MMSN at distances ≳ 1 AU (independent of accretion rate) and even at ≲ 1 AU (if the disk is heated by accretion rates ≳ 3 × 10⁻⁷ M_☉ yr⁻¹). At low accretion rates and/or large radii, strong-field MRI (if active) is also sufficient to drive statistical instability if Σ ≳ Σ_MMSN. And we show that beyond a few AU, the combination of gravitoturbulence, convection, and Toomre instability at low T means there may be no disk temperatures at which $P_{\rm frag}^{\rm int}\ll 1$ in a modest-density disk.

Ultimately, regardless of our (admittedly speculative) discussion of theoretical models for the sources of turbulence in protoplanetary disks, the key question is empirical. Are the actual Mach numbers in such disks sufficient, for their Q values, to be "interesting" here? This remains an open question. However, Hughes et al. (2011) present some early indications of detection of turbulent linewidths in two protoplanetary disks, with inferred Mach numbers of ∼0.1 and ∼0.4. Although uncertain, these essentially bracket the most interesting regime of our calculations here! Because of the exponential dependence of stochastic collapse on Mach numbers, future observations which include larger statistical samples and more/less massive disks, as well as constraints on whether the turbulence appears throughout the disk (since it is the midplane Mach numbers that matter most for the models here), will be critical to assess whether the processes described in this paper are expected to be relatively commonplace or extremely rare.

We also predict the characteristic mass spectrum of fragmentation events. Sub-sonic turbulence produces a narrow mass spectrum concentrated around the Toomre mass ∼μ² Q² M_disk; angular momentum and shear suppress larger-scale collapse, while thermal (and magnetic) pressure suppresses the formation of smaller-scale density fluctuations. For typical mass ratios μ ∼ 0.1 in the early stages of disk evolution, this corresponds reasonably well to the masses of giant planets. For smaller mass ratios μ ∼ 0.01, which should occur at somewhat later stages of evolution, this implies that direct collapse to Earth-like planet masses may be possible. Of course, such collapse will carry whatever material is mixed in the disk (i.e., light elements), so such a planet would presumably lose its hydrogen/helium atmosphere as it subsequently evolved (see references in Section 1). As the turbulence approaches transsonic, the mass spectrum becomes much more broad. As shown in Paper II, this owes to the greater dynamic range in which turbulence is important; in the limit $\mathcal {M}\rightarrow \infty$, the mass spectrum approaches a power law with equal mass at all logarithmic intervals in mass (up to the maximum disk Jeans mass). Thus, within a more turbulent disk, direct collapse to a wide range of masses—even at identical disk conditions and radii—is expected. This may explain observed systems with a range of planet masses within narrow radii (e.g., Carter et al. 2012), and it also predicts a general trend of increasing average collapse mass with distance from the central star. However, we are not accounting for subsequent orbital evolution here, and subsequent accretion (after collapse) will modify the mass spectra.

We discuss and show how these predictions change with different turbulent velocity power spectra, different gas equations of state, including or excluding magnetic fields, changing the disk mass profile, or allowing for (quite large) deviations from lognormal statistics in the density distributions. However, these are generally very small corrections and/or simply amount to order-unity re-normalizations of the predicted object masses, and they do not change any of our qualitative conclusions. Because of the very strong dependence of fragmentation on Mach number, the critical Mach numbers we predict as the threshold for statistical instability are insensitive to most changes in more subtle model assumptions.¹⁵

Throughout, we restrict our focus to the formation of self-gravitating regions.¹⁶ Following the subsequent evolution of those regions (collapse, fragmentation, migration, accretion, and possible formation of planets) requires numerical simulations to treat the full nonlinear evolution (see, e.g., Kratter & Murray-Clay 2011; Rice et al. 2011; Zhu et al. 2012; Galvagni et al. 2012, and references therein). Ideally, this would be within global models that can self-consistently follow the formation of these regions. However, this is computationally extremely demanding. Even ignoring the detailed physics involved, the turbulent cascade must be properly followed (always challenging), and most important, since the fluctuations of interest can be extremely rare, very large (but still high-resolution) boxes must be simulated for many dynamical times. As discussed in Section 1, this has led to debate about whether or not different simulations have (or even can be) converged. Most of the longest-duration simulations to date have been run for times ∼1000 Ω⁻¹, which for plausible $\mathcal {M}_{c}$ and Q may be shorter by a factor of ∼10³–10⁵ than the timescale on which of order a single event is expected to occur in the entire disk. But certainly in the example of planet formation, a couple of rare events are "all that is needed," so this is an extremely interesting case for future study.

We thank Jim Stone, Eugene Chiang, and Eliot Quataert for insightful discussions that helped inspire this paper. Support for P.F.H. was provided by NASA through Einstein Postdoctoral Fellowship Award Number PF1-120083 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060.

Recall, from the text, in the case of disk irradiation by a central solar-type star, the effective temperature is

$$\begin{equation} T_{\rm eff,\,\ast } \approx \left (\frac{\alpha _{T}\,R_{\ast }^{2}}{4\,r_{\ast }^{2}} \right )^{1/4}\,T_{\ast } \end{equation} \tag{ A1 }$$

(Chiang & Goldreich 1997), with α_T ≈ 0.005 (r_*/AU)⁻¹ + 0.05 (r_*/AU)^2/7, R_* = R_☉, and T_* = 6000 K. In this regime the external radiation produces a hot surface dust layer that re-radiates ∼1/2 the absorbed light back into the disk, maintaining T_mid; if the disk is optically thick to the incident and re-radiated emission, this gives an approximate $T_{\rm mid,\,\ast }\approx T_{\rm eff,\ast }^{1/4}/2^{1/4}$. Meanwhile, accretion produces an effective temperature

$$\begin{equation} T_{\rm eff,\,acc} \approx \left [\frac{3}{8\pi }\,\frac{\dot{M}\,\Omega ^{2}}{\sigma _{\rm SB}} \right ]^{1/4}. \end{equation} \tag{ A2 }$$

In the optically thick limit, this is just related to T_mid by $T_{\rm mid,\,acc}^{4}\approx (3/4)\,(\tau _{R} + 2/3)\,T_{\rm eff,\,acc}^{4}$, where τ_R is the Rosseland-mean optical depth, τ_R = κ_R(T_mid) Σ (so T_mid, acc is determined implicitly).

In the text, Figures 4 and 5 simply use the optically thick relations for T_mid and, when accretion is included, interpolate with $T_{\rm mid}^{4} \approx T_{\rm mid,\,\ast }^{4}+T_{\rm mid,\,acc}^{4}$. We approximate the opacities with the simple κ_R ∼ 5 cm² g⁻¹ at T > 160 K, and κ_R ∼ 2.4 × 10⁻⁴ T² cm² g⁻¹ K⁻² at lower temperatures (approximately what is obtained with a galactic gas-to-dust ratio; see Adams & Shu 1988; Bell & Lin 1994).

We then follow Rafikov (2007) and estimate the cooling time as

$$\begin{equation} t_{\rm cool} = \frac{\Sigma \,c_{s}^{2}}{\sigma _{\rm SB}\,T_{\rm mid}^{4}}\,f(\tau), \end{equation} \tag{ A3 }$$

where $c_{s} = \sqrt{k_{B}\,T_{\rm mid}/\mu }$ with μ ≈ 2.3 appropriate for molecular, dusty gas. The function f(τ) is given by the interpolation between the convective and radiative terms at the photosphere

$$\begin{eqnarray} f(\tau) &= \chi \,\tau ^{\eta } + \phi \,\tau ^{-1}, \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} \eta &= \frac{4\,\gamma ^{-1}\,(\gamma -1)}{1+\alpha +\beta \,\gamma ^{-1}\,(\gamma -1)}, \end{eqnarray} \tag{ A5 }$$

where χ and ϕ are constants; from the detailed estimates therein ϕ ≈ 1 and χ ≈ 0.19–0.31, depending on the disk parameters, so we assume χ = 0.31 to be conservative (since this gives larger cooling times). This interpolates between the optically thick limits (dominated by χ τ^η) and optically thin cases (ϕ τ⁻¹, where the cooling flux becomes ${\propto} \tau \,\sigma _{\rm SB}\,T_{\rm mid}^{4}$).

The scaling index η is determined from the temperature gradient to the photosphere for a disk in vertical hydrostatic equilibrium, under the assumption that over some (limited) local temperature range the opacity κ can be approximated by κ ≈ κ₀ P^α T^β, which is valid for our assumptions. We assume γ = 7/5 in both this and the turbulent density fluctuation calculation in Figures 3–5. Based on the scaling of opacities in Semenov et al. (2003), at T < 160 K, α = 0 and β ≈ 2, so η ≈ 7/11; at 160 < T < 1500 K, α = 0 and β ≈ 1/2–1, so η ≈ 7/8–7/9 (we adopt 7/9, but this makes no significant difference); and at T > 1500 K (when all grains sublimate and molecular opacity dominates), α ≈ 2/3 and β ≈ 7/3, so η ≈ 3/7.

In the text, we restrict to these simple estimates because the quantities of interest are fairly uncertain. We can, however, examine a more detailed approximation here. First, we take the opacities κ_R(T_mid) from the full tabulated values in Semenov et al. (2003). Second, a more accurate estimate of the midplane temperature is given by solving the implicit equation

$$\begin{equation} T_{\rm mid}^{4} = \frac{3}{4}\,\left [\tau _{V} + \frac{4}{3}+\frac{2}{3\,\tau _{V}} \right ]\,T_{\rm eff,\,acc}^{4} + \left [1 + \frac{1}{\tau _{V}} \right ]\,T_{\rm eff,\,\ast }^{4}, \end{equation} \tag{ A6 }$$

where τ_V = τ_V(T_mid) is the vertical optical depth from the midplane, τ_V = κ_R(T_mid) Σ/2. This allows for an appropriate interpolation between the optically thick case and the case where the disk is optically thin to its own re-radiation. Third, we can switch between the f(τ) above, appropriate for a convective disk, and f(τ) ≈ τ + τ⁻¹ appropriate for a purely radiative, convectively stable disk, when the disk falls below the (temperature-dependent) criteria for convective instability ∇₀ ⩾ ∇_ad with ∇₀ ≡ (1 + α)/(4 − β) and ∇_ad ≡ (γ − 1)/γ (see Lin & Papaloizou 1980; Rafikov 2007), where α and β depend on T (using the full explicit derivatives from the opacity tables).

Figures 6–8 repeat our calculations from the main text with this more detailed temperature calculation. We find that the quantitative results in Figures 3–5 are all changed at the factor ≲ 2 level and the qualitative conclusions are completely unchanged. The sense of the quantitative change tends to slightly expand the regions of parameter space where statistical instability and turbulence-promoted fragmentation can occur. The more detailed opacity calculation imprints some small features on the parameter space, the most significant of which is the elimination of predicted temperatures much larger than ∼1500 K at small radii (because grains sublimate and cooling becomes optically thin), but the "gravito-turbulent" thresholds shift with the predicted temperatures so the statistical stability is essentially identical.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here, we review the basic framework of the model developed in Paper I and Paper II, specifically some key equations needed to reproduce the results in this paper. Readers interested in a full derivation and explanation of these equations should see Paper II.

Consider, for simplicity, the isothermal (lognormal) case: if density fluctuations are lognormal, then the variable δ(x) ≡ ln [ρ(x)/ρ₀] + S/2, where ρ(x) is the density at a point x, ρ₀ is the global mean density, and S is the variance in ln ρ, is normally distributed according to the PDF:¹⁷

$$\begin{equation} P_{0}(\delta \,|\,S) = \frac{1}{\sqrt{2\pi \,S}}\,\exp {\left(-\frac{\delta ^{2}}{2\,S} \right)}. \end{equation} \tag{ B1 }$$

More generally, we can evaluate the field δ(x | R), which is the δ(x) field averaged around the point x with some window function of characteristic radius R; this is also normally distributed (see Paper II; Appendix F), with a variance at each scale S(R) that is directly related to the density power spectrum.

Paper I derives the "first-crossing" distribution for the general form of these fields. This corresponds to the mass and initial size spectrum of regions that are sufficiently dense so as to be self-gravitating averaged on the scale R (specifically defined as the largest scale on which the region is self-gravitating, i.e., excluding bound sub-units already counted "within" the parent, although these can be counted separately if desired). This corresponds to the statistics of regions where δ(x | R) > B(R), where B(R) (the "barrier") is some (scale-dependent) critical value. In Paper I we derive S(R) and B(R) from simple theoretical considerations for all scales in a galactic disk and/or molecular cloud. However, the derivation proceeds almost identically for a protoplanetary disk, following the most general form presented in Paper II, which we outline here.

It is well established that the contribution to density variance from the velocity variance on a given scale goes as $S\approx \,\ln {(1+b^{2}\mathcal {M}^{2})}$, where $\mathcal {M}$ is the Mach number. For a given turbulent power spectrum, then, S(R) is determined by summing the contribution from the velocity variance on all scales R' > R:

$$\begin{eqnarray} S(R,\,\rho) &= \int _{0}^{\infty } |\tilde{W}(k,\,R)|^{2} \ln \left [1 + \frac{b^{2}\,v_{t}^{2}(k)}{c_{s}^{2}(\rho,\,k) + \kappa ^{2}\,k^{-2}} \right ] {d}\ln {k}, \nonumber\\ \end{eqnarray} \tag{ B2 }$$

where W is the window function for the smoothing,¹⁸ v_t(k) is the turbulent velocity dispersion averaged on a scale k (trivially related to the turbulent power spectrum), c_s is the thermal sound speed (both c_s and S can depend locally on ρ if the gas is not isothermal), and b is the fraction of the turbulent velocity in compressive (longitudinal) motions (discussed below). Here κ is the epicyclic frequency; since we are interested in Keplerian disks, we take κ ≈ Ω, the disk orbital frequency. Note that on large scales, angular momentum (κ² k⁻²) enters in a similar way to $c_{s}^{2}$ and suppresses fluctuations, which follows directly from the form of the dispersion relation for density perturbations (e.g., Lin et al. 1969; Toomre 1977; Lau & Bertin 1978); accounting for this is necessary to ensure mass conservation.

Since we are interested in the formation of self-gravitating gas objects, we define B(R) corresponding to the critical density averaged on a given scale, ρ_crit(R), at which an overdensity will collapse. Given δ(R) ≡ ln [ρ(R)/ρ₀] + S/2 defined above, then B(R) follows from the dispersion relation for a density perturbation in a disk with self-gravity, turbulence, thermal and magnetic pressure, and angular momentum/shear (Vandervoort 1970; Aoki et al. 1979; Elmegreen 1987; Romeo 1992):

$$\begin{equation} B(R) = \ln {\left(\frac{\rho _{\rm crit}(R)}{\rho _{0}} \right)} + \frac{S(R)}{2}, \end{equation} \tag{ B3 }$$

where ρ_crit is the critical density above which a region is self-gravitating. This is the (implicit) solution to

$$\begin{eqnarray} &&\frac{\rho _{\rm crit}(R)}{\rho _{0}} \equiv \frac{Q}{2\,\tilde{\kappa }}\,\left(1+\frac{h}{R} \right) \left [ \frac{\sigma _{g}^{2}(R,\,\rho _{\rm crit})}{\sigma _{g}^{2}(h,\,\rho _{0})}\,\frac{h}{R} + \tilde{\kappa }^{2}\,\frac{R}{h}\right ], \quad \end{eqnarray} \tag{ B4 }$$

where ρ₀ is the mean midplane density of the disk, h is the disk scale height, $\tilde{\kappa }\equiv \kappa /\Omega =1$ for a Keplerian disk, and Q ≡ (σ_g[h, ρ₀] κ)/(π G Σ_gas) is the Toomre Q parameter. The "total" dispersion σ_g is

$$\begin{equation} \sigma _{g}^{2}(R,\,\rho) = c_{s}^{2}(\rho) + \left\langle v_{t}^{2}(R) \right\rangle + v_{\rm A}^{2}(\rho,\,R). \end{equation} \tag{ B5 }$$

The map between scale R and the total mass in the collapsing region is

$$\begin{eqnarray} &&M(R) \equiv 4\,\pi \,\rho _{\rm crit}\,h^{3}\, \left [\frac{R^{2}}{2\,h^{2}} + \left (1+\frac{R}{h}\right)\,\exp \left (-\frac{R}{h}\right)-1 \right ].\nonumber\\ \end{eqnarray} \tag{ B6 }$$

It is easy to see that on small scales, these scalings reduce to the Jeans+Hill criteria for a combination of thermal, turbulent, and magnetic support, with M = (4π/3) ρ_crit R³; on large scales it becomes a Toomre-like criterion with M = πΣ_crit R².

For any B(R) and S(R), Paper I shows that the instantaneous "first-crossing" mass function (i.e., instantaneous mass function of collapsing objects, uniquely defined to resolve the "cloud-in-cloud" problem) is

$$\begin{equation} \frac{{d}n}{{d}M} = \frac{\rho _{\rm crit}(M)}{M}\,f_{f}(M)\,\left |\frac{{d}S}{{d}M} \right |, \end{equation} \tag{ B7 }$$

where f_f(S) is a function shown by Zhang & Hui (2006) to be the solution of the Volterra integral equation:

$$\begin{eqnarray} &&f_{f}(S) = \tilde{g}_{1}(S) + \int _{0}^{S}\,{d}S^{\prime }\,f_{f}(S^{\prime })\,\tilde{g}_{2}(S,\,S^{\prime }), \end{eqnarray} \tag{ B8 }$$

with

$$\begin{eqnarray} \tilde{g}_{1}(S) &=& \left [2\,\left |\frac{dB}{dS}\right | +\frac{B(S)}{S} \right ]\,P_{0}(B(S)\,|\,S), \end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray} \tilde{g}_{2}(S,\,S^{\prime }) &=& \left [-\frac{B(S)-B(S^{\prime })}{S-S^{\prime }} -2\,\left |\frac{dB}{dS}\right |\right ]\\ \nonumber &&\times \, P_{0}[B(S)-B(S^{\prime })\,|\,S^{\prime }-S]. \end{eqnarray} \tag{ B10 }$$