{'fraction_non_alphanumeric': 0.02963282785115738, 'fraction_numerical': 0.01541118522499978, 'mean_word_length': 4.768820210440706, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gravitational Instabilities in Circumstellar Disks

[Abridged] Star and planet formation are the complex outcomes of gravitational collapse and angular momentum transport mediated by protostellar and protoplanetary disks. In this review we focus on the role of gravitational instability in this process. We begin with a brief overview of the observational evidence for massive disks that might be subject to gravitational instability, and then highlight the diverse ways in which the instability manifests itself in protostellar and protoplanetary disks: the generation of spiral arms, small scale turbulence-like density fluctuations, and fragmentation of the disk itself. We present the analytic theory that describes the linear growth phase of the instability, supplemented with a survey of numerical simulations that aim to capture the non-linear evolution. We emphasize the role of thermodynamics and large scale infall in controlling the outcome of the instability. Despite apparent controversies in the literature, we show a remarkable level of agreement between analytic predictions and numerical results. We highlight open questions related to (1) the development of a turbulent cascade in thin disks, and (2) the role of mode-mode coupling in setting the maximum angular momentum transport rate in thick disks.


INTRODUCTION
Protostellar and protoplanetary disks stand at the center of our study of origins. They regulate the formation of stars and the planetary systems that they host. Over the last three decades, the observational study of protostellar disks has evolved from quantifying hints of infrared excess (Kenyon & Hartmann 1987) to resolving ring-like structures at 1AU scales (ALMA Partnership et al. 2015). Long wavelength studies with facilities such as the SMA, Plateau de Bure, CARMA, and now ALMA, combined with sensitive infrared measurements from IRAS and Spitzer enable detailed studies of the sizes, temperatures, and masses of a variety of disks in young clusters. While many uncertainties remain in translating high resolution interferometric measurements into reliable disk properties, we can begin to see the various stages of star and planet formation occurring in these disks.
Theoretical studies have made similar progress, from simple, 1-D diffusive models to threedimensional (magneto)hydrodynamics simulations that incorporate radiative transfer. These advances have been enabled by the burgeoning computational resources at our disposal. While the most complex of these models provide many insights, which we review in detail, simple models prove valuable to this day.
In this review we discuss the role that self-gravity plays in regulating the growth and evolution of protostellar and protoplanetary disks. Self-gravity in general refers to the mutual gravitational influence of relatively diffuse gas on itself, in contrast to the gravity exerted on the gas by an external body. In this context, self-gravity refers to the effect of the gravitational field of the disk itself, rather than the protostar, on the evolution of the system. Although the dividing line between protostellar and protoplanetary disks is fuzzy, we use the former to refer to younger disks surrounding protostars that are not yet near their final mass, while the latter refers to the remnant disk in which the planet formation process completes (of course it may begin much earlier). The physics we discuss is relevant for all stellar masses, but much of our scalings are for sun like stars.
We focus our attention on hydrodynamics, and thus we do not discuss the important role that magnetic fields and MHD effects have on the evolution of protostellar disks. It is well known that the hotter parts of protostellar disks are subject to the magneto-rotational instability (Balbus & Hawley 1998), that can provide the bulk of angular momentum transport. On the other hand, the colder, outer disks, which are the regions most likely to become gravitationally unstable, are probably less affected by MHD turbulence, given their low ionization state. Very few numerical studies have considered the combined effects of MHD and gravitational instabilities in disks (Fromang et al. 2004a,b), whose interplay might be important in protostellar disk evolution.
In Section 2 we begin with a brief observational overview with an eye towards the applicability of the physics of self-gravitating fluids to observed systems. In Section 3 we provide an in-depth review of the nature of self-gravity and gravitational instabilities, including the non-linear outcome of the instability: fragmentation of a disk. The physics discussed here is applicable to accretion disks in many different astrophysical contexts. In Section 4 we hone in on disks around young stars, and consider how they fit into the various regimes outlined in Section 3. In Section 5 we discuss the relevance of gravitational instabilities to the formation of planets both directly (gas-collapse) and indirectly (concentration of dust particles). We conclude with a brief summary and prospects for future.

Observational Evidence for Self-Gravitating Disks
Before reviewing the theoretical underpinnings of the role of self-gravity in protostellar and protoplanetary disk evolution, we address the observational evidence that this physics is relevant. As we detail in later sections, understanding the behavior of self-gravitating disks is complex and depends on a range of disk properties including size, surface density, temperature, and thermal physics. The standard reference for quantifying the degree to which a disk is self-gravitating is the parameter Q, defined as (Toomre 1964): where cs is the sounds speed, κ is the epicyclic frequency (which is equal to Ω in a Keplerian disk), and Σ is the disk surface density. As Q → 1 self-gravity becomes increasingly important. We discuss the origin of Q in more detail in section 3. To translate this constraint into one motivated by observations, we note that Q can be written as: where M * is the stellar mass, M d is the disk mass, H is the disk scaleheight, and r the radius.
We subsume numerical factors of order unity, which depend on the surface density profile, into a pre-factor f . While we cannot escape the dependence on the disk aspect ratio, H/r, which depends on the temperature, we can use this relation to obtain an order of magnitude estimate for the minimum mass required for a disk to be self-gravitating. Colder, thinner disks are more prone to self gravity. Star forming regions typically do not have temperatures below 10K, so we take that as a lower limit. What about disk radii? Many lines of argument, from measures of cloud angular momentum (Goodman et al. 1993), to stellar interactions (Adams et al. 2006), to direct observations in both absorption (Eisner et al. 2008) and millimeter wavelengths (Andrews et al. 2011) suggest that 10 2 − 10 3 AU are reasonable limits for disk radii. At constant temperature, the aspect ratio increases with disk radius as r 1/2 , so we use 100AU to set H for our scaling. The requirement that Q < 1 (implying gravitational instability, see below) then implies: Thus for self-gravity to be important, i.e. for Q to be near unity given optimistic constraints on temperature, we expect disk-star mass ratios greater than ≈ 10 −2 .
With this benchmark in mind, we now consider the evidence for disks with mass ratios in this regime or higher.
limitations are related to the fact that most disk masses are measured by proxy through detection of millimeter wavelength fluxes associated with dust grains of up to a millimeter in size. Based on studies of the ISM, the ratio of gas to dust is thought to be of order 10 2 , but this is not well confirmed in disks. Secondly, because larger grains and rocks do not emit at millimeter wavelengths, these data do not reveal the potentially sizeable fraction of dust which has already grown, something for which there is already strong evidence (e.g., Pérez et al. 2012). Moreover, the conversion of a flux to a <mm-size dust mass requires adopting an opacity, which is itself dependent on unknowns like the grain-size distribution, and grain properties such as chemical composition and porosity (see, e. g. Natta et al. 2004;Ricci et al. 2010). Finally, mass measurements from dust are only possible if the dust is optically thin. If the dust is optically thick, which is possible in the inner regions of disks (soon to be resolved by ALMA) only a lower limit on the mass is possible.
In addition to the uncertainties generated by the dust properties,  have highlighted the estimation errors introduced by contamination of disk flux measurements due to the presence of a core at early phases, leading to overestimated masses at early times. Similarly, as disks grow to larger radii, high resolution interferometric studies can also resolve out large scale structures where much of the mass resides. This error leads to a more severe underestimate of disk masses.
With these caveats in mind, we review the best estimates of disk masses.
2.2.1. Observational Results. The most comprehensive survey for disk masses at this time is that of Andrews et al. (2013), who survey all known members of the nearby star forming region Taurus. These disks are mostly Class II sources, however some class I sources are included. Andrews et al. (2013) compiles photometry at 1.3mm using the SMA to derive both disk masses, and disk-star-mass ratios. The stellar masses are estimated using the isochrones from D' Antona & Mazzitelli (1997); Baraffe et al. (2002) and Siess et al. (2000). Disk masses are estimated from the submillimeter flux assuming a constant distance of 140pc, dust-gas ratio ζ = 0.01, dust opacity: κν = 2.3cm 2 /g and a disk temperature constant with radius, which scales with stellar luminosity as T d = 25(L * /L ) 1/4 K. The masses are derived from the following equation for optically thin, isothermal dust: log M d = log Fν + 2 log d − log(ζ · κν ) − log Bν ( T d ), where Fν is the observed flux, Bν is the Planck function and d is the distance to the source. It is apparent from Andrews et al. (2013) figure 11 that most disk masses are far too low in mass by 2-3 Myr (predominantly class II stage) for GI to be relevant. Note additionally that the assumed outer disk temperatures for this sample are higher than our optimistic (for triggering GI) estimate of 10K. Depending on order unity factors, at most 20%, and probably more like 10% of Class II disks are plausibly massive enough for self-gravity to be important. A self-consistent treatment of the disk temperatures in our instability threshold (Equation (3)) and mass estimates would yield critical mass ratios closer to 0.1, higher than nearly every disk in the sample.
As this sample is dominated by Class II disks, we may conclude that both angular momentum transport and planet formation reliant on strong self-gravity are rare at best for systems which have evolved to the Class II stage.
What about evidence for younger, massive Class I and Class 0 disks? Unfortunately at the time of this writing the evidence remains sparse. Class 0 and I sources are deeply embedded and thus intrinsically harder to both identify and measure. Eisner (2012) and Sheehan & Eisner (2014) have isolated and modelled in detail the Class I sources in Taurus. They find that the masses are typically higher than in the Class II phase by a factor of two: the median mass of Class I disks is 0.008M compared to 0.005M for older sources. Among Class I sources, roughly 10% have masses higher than 0.01M (Eisner 2012), at the margins of susceptibility to strong self-gravitating effects. Some individual sources have been measured to have relatively massive and compact disks in the Class 0/I phase, for example IRAS 16293-2422B (Rodríguez et al. 2005, or WL12 (Miotello et al. 2014).
A recent survey by Mann et al. (2015) also hints that younger sources may indeed have more massive disks. These authors find that in NGC2024 (thought to be < 1Myr old as compared to 2-3 Myr for Taurus) as many as 10% of disks have masses greater than 0.1M , with 20% above 0.01M . Despite the 2σ significance of higher mass measurements, selection biases alone (related to stellar mass and multiplicity) could account for the increase if the cluster properties of Taurus and NGC2024 are similar.
Finally, we turn to observations of disks in the Class 0 stage. The most conclusive evidence for a relatively massive disk comes from Tobin et al. (2013). Using the SMA and CARMA, they detect an edge on ≈ 150 − 200AU disk surrounding the Class 0 protostar L1527. Using a similar optically thin dust model they derive a mass comparable to Class I sources of 0.075M , however the protostellar mass is thought to be < 0.5M , and thus the disk-star mass ratio is easily high enough for self-gravity to be important. There are several other tentative detections of Class 0 disks (see e.g. Looney et al. 2000;Enoch et al. 2009). A larger survey in Perseus by Tobin et al. (2015) resolved several more likely candidates and found that more than 50% of 9 sources were more consistent with envelope + disk models rather than pure envelope models. Additionally several sources show evidence for rotation (though not necessarily Keplerian motion). The masses measured in this survey cannot reliably separate out all envelope contamination, thus the increase in masses (up to 0.1M ) is not necessarily indicative of strong self-gravity. Nevertheless it is quite plausible that disk masses are higher in the Class 0 phase. Maury et al. (2010) have made strong claims against the presence of disks in the Class 0 phase, but the well-resolved L1527 was included in their sample, suggesting that perhaps the data was insufficient to rule out some extended disks. Tobin et al. (2013) had slightly better resolution, and also a different beam alignment relative to the disk direction, which might explain the Maury et al. (2010) non-detection.
While dust mass measurements are most common, recent efforts have been made to directly detect molecular gas tracers to constrain the total mass and dust-to-gas ratio. Bergin et al. (2013) use Herschel data to detect HD in TW Hya, finding a relatively large disk mass for a nominally 10 Myr system. Williams & Best (2014) used optically thick and thin CO isotopologues and conclude that most Class II disks have dust-to-gas ratios up to a factor of 10 below ISM values, implying that either planet formation is well underway, or total disk masses are even smaller.
The above discussion focused on solar type and Herbig AeBe stars, for which there is more data, but there is some evidence that high mass stars host massive, self-gravitating disks (Cesaroni et al. 2007 and references therein). Precise measurements remain elusive because high mass stars are typically further away, and thus harder to resolve. Additionally, the only evidence for disks coincides with the equivalent of the Class 0-I phases, when the stars are deeply embedded in a gas envelope. For example, rotationally supported disks are difficult to distinguish from rotationally flattened envelopes (Johnston et al. 2011). One of the more convincing examples (Shepherd et al. 2001) reveals a rotating massive disk around a B-star that is likely gravitationally unstable based on the estimated parameters. Krumholz et al. (2007) has shown that a fully functional ALMA may be able to validate these observations using molecular lines.

Direct Evidence for Disk Self-Gravity
While the population at large may not show strong evidence for disk self-gravity, a handful of individual sources do show indications of the spiral arm structure typically associated with spiral density waves. These wave may be triggered either by a young planet, or by disk self-gravity. Muto et al. (2012) and Wagner et al. (2015) find evidence for spiral structure in scattered light images of transition disks. Self-gravity might be at play, but the low sub-mm masses make this questionable. Benisty et al. (2015) also find evidence for spiral structures in scattered light images, but claim a planetary origin. ALMA data might settle these two theories by providing disk mass constraints, but at present neither theory can be ruled out (Dong et al. 2015).

Indirect Evidence for Disk Self-Gravity
While observational examples of disks with strong self-gravity remain scarce, there are other indirect sources of evidence from comprehensive studies of star formation. There are a variety of reasons to expect that disks may be more massive at earlier evolutionary states. First, since the discovery of the first statistically significant sample of protostars, observers and theorists alike have pointed out the "luminosity problem" (Hartmann & Kenyon 1996). This problem arose due to the mismatch between the observed protostellar luminosities, of the order of a few L , and the expected luminosities derived from basic energetic arguments that predict that the accretion luminosity: where Rin is the inner disk radius andṀ is the infall rate, estimated based on a star formation timescale of roughly 1Myr. This over prediction of the typical luminosity has been used as evidence for large, short-lived accretion events, either bursts, or more continuous but rapid events early in a protostars history. A variety of theories have been proposed for producing these bursts. Before considering specific theories, we note the following simple scaling argument. Consider a star forming core or turbulent filament. Observational evidence suggest that these bodies eventually undergo a modified form of free fall collapse when self-gravity overcomes magnetic or turbulent pressures. Thus the rate at which material will fall onto any disk like structure formed in the collapse is dictated by free fall collapse. The oversimplified isothermal inside out collapse mode of  provides an order of magnitude estimate for this infall rate: Somewhat more realistic Bonnor Ebert sphere infall rates are the same order of magnitude, and are consistent with the timescale on which star formation occurs. In contrast, the observed accretion rates from disks on to stars in the Class I -Class II phase range from 10 −9 − 10 −7 M /yr, with the majority of measurements at the lower end (Gullbring et al. 1998). There are two possible means to resolve the mismatch between infall and stellar accretion rates. Either disks at earlier times are more massive, and thus can process more material at higher rate. This option is somewhat suggestive of self-gravitational effects, but not required. Alternatively if initially disk-star accretion rates are not higher, the mass of the disk will rise, bringing it ever closer to the regime where self-gravity is important. It seems unlikely that the infall and disk-star accretion rates, which rely on entirely different physics (free fall and cloud turbulence vs either self-gravity or magnetic fields), could conspire to match up at all times in order to maintain a constant, tiny fraction of the total mass in the disk. Most of the stellar mass is thought to pass through the disk at some point, making such a balance improbable. Armitage et al. (2001) first suggested that stellar accretion and FU Orionis type outbursting systems might be generated by an interplay between gravitational instability driven accretion and magnetically dominated accretion driven by MRI turbulence. Vorobyov & Basu (2006) suggest that these same outbursts could be triggered by clumps formed due to gravitational instability in the outer disk that are tidally disrupted as they migrate inwards, and subsequently accrete onto the star. See the chapter by Hartmann in this Annual Reviews (and recently, the Protostars and Planet VI chapter, Audard et al. 2014) for more details on variable accretion onto pre-main-sequence stars.
One final piece of evidence for higher disk masses relies on observations of exoplanets. Najita & Kenyon (2014) have shown that typical protoplanetary disks do not have enough observed dust mass to account for the mass of planets implied by occurence rates derived from Kepler. They infer that planet formation, and thus grain growth likely commenced at even earlier times. This implies that the observed dust masses are missing a substantial amount of solids. Accounting for this missing mass might well push a larger fraction of disks into the self-gravitating range.
In conclusion we see that while inferred disk masses are often lower than the nominal requirement for strong self-gravity, it is by no means ruled out. Given the many ways to underestimate disk masses, the observational evidence demands the study of disk self-gravity .
3. Physics of self-gravity and gravitational instability GI: Gravitational Instability

SPH: Smoothed Particle Hydrodynamics
In the previous section we have shown that disks may be subject to the effects of self-gravity when the ratio M d /M * > 10 −2 . We first address how self-gravity affects the basic hydrostatic structure of the disk, and subsequently explore how gravitational instability, hereafter GI, may affect disk structure and evolution. In this section we focus on the physical nature of self-gravity and GI, which applies to disks in a wide range of astrophysical contexts. We begin by discussing isolated disks with a fixed surface density profile for simplicity, and subsequently consider the changes that arise in the presence of infall and external heating. In Section 4, we discuss the extent to which such physical processes apply to actual protostellar disks.

Effects of self-gravity on the structure of accretion disks
Consider a thin axisymmetric disk surrounding a protostar of mass M . Let ρ(r, z), Σ(r) and cs(r) be the volume density, the surface density and sound speed of the gas as a function of cylindrical radius r, and height above the disk midplane z, respectively. Standard, non self-gravitating disk models are based on the assumptions of centrifugal balance in the radial direction and hydrostatic balance in the vertical direction. Both equations must be modified to account for self-gravity.
A standard non self-gravitating disk in hydrostatic balance is characterised by a Gaussian profile for ρ, where the disk scale height is given by where Ω k = GM * /r 3 is the Keplerian angular velocity and we have assumed the disk to be vertically isothermal. Self-gravity modifies vertical hydrostatic balance. In the limit where selfgravity dominates, the vertical density profile is given by (Spitzer 1942) where Hsg = c 2 s πGΣ .
In order to neglect self-gravity in the vertical direction, we require that Hnsg/Hsg 1. This immediately translates into the condition csΩ k πGΣ 1.
From this requirement alone we thus recover the usual stability parameter for GI (Toomre 1964). This means that whenever self-gravity is important in the vertical direction, the disk is near the instability threshold for GI. This constraint is also equivalent to the constraint on disk mass presented in Equation (2) in Section 2. In the general case where both contributions (of the central star and of the disk) are important, there is no analytical solution to the hydrostatic balance equation, but a simple and accurate interpolation formula between the two results above has been found by Bertin & Lodato (1999) (see also Lodato 2007;Muñoz et al. 2015): We now consider centrifugal balance. To lowest order, centrifugal balance prescribes the angular velocity of the disk Ω to be Keplerian, Ω(r) = GM /r 3 . Even in the absence of self-gravity, a small correction is required to account for pressure gradients, which cause the gas to orbit at a slightly sub-Keplerian speed. This sub-Keplerian motion has a dramatic effect on the dynamics of small solids in the disk, whose orbits are Keplerian due to the absence of pressure forces. Taking these corrections into account, we have which gives sub-Keplerian rotation because in general P is a decreasing function of radius. Deviations from Keplerian motion due to the pressure gradient are of order (H/r) 2 . Self-gravity also modifies centrifugal balance by a force term which is of the order of GM disk /r 2 (a more detailed calculation, obtained by solving Poisson's equation, can be found in Bertin & Lodato 1999). Clearly, in order to dominate the rotation curve, the disk mass needs to be implausibly large for protostellar disks. However, it is worth noting that when the disk is marginally gravitationally stable (Q ≈ 1), the disk mass is of the order of (H/r)M (Equation (2)) and thus the deviation from Keplerian motion due to self-gravity is of the order of H/r, stronger than the effect of pressure gradients, though both solids and gas feel the change in Ω(r) due to disk self-gravity.

Dispersion relations
The most important effect of disk self-gravity on the dynamics of protostellar disks is connected with the development of GI. The linear stability of self-gravitating disks has been investigated for decades in the context of galactic dynamics. In the tightly wound approximation, the WKB dispersion relation for an infinitesimally thin disk is the famous Lin & Shu (1964) relation: where ω is the wave frequency, k is the radial wavenumber, m is the azimuthal wavenumber (that is, the number of spiral arms produced by the instability), and κ is the epicyclic frequency, which, for a Keplerian disk is simply Ω. We now see that the aforementioned Q from Equation (1) arises as the threshold for exponential growth of axisymmetric modes, where stability requires The above dispersion relation and stability criterion have been derived under several important limiting approximations. The first is that the equilibrium structure of the disk is axisymmetric, the second (and most important) is that the disk is infinitesimally thin. It can be shown that finite disk thickness (Vandervoort 1970;Bertin 2000;Binney & Tremaine 1987) has the effect of stabilizing the disk by diluting the self-gravity term in Eq. (13). A simple way to account for this effect is to modifiy eq. (13) as follows: that can be expanded to first order to give where now the stabilizing effect of disc thickness becomes apparent, as an additional 'pressure-like' term in the dispersion relation. By inserting Eq. (11) in Eq. (16) we finally get From Equation (17) we see that the stability of the disk depends only on the single parameter Q, and one can show that the marginal stability value decreases to Q ≈ 0.6, for a Keplerian disk.
In the context of galactic dynamics, it was recognized many years ago (Ostriker & Peebles 1973;Hohl 1971Hohl , 1973) that disks that are locally stable according to the Q criterion, might still generate large scale, spiral waves. This is due to the fact that such global modes are not captured by the WKB tightly wound approximation. However, a WKB description of such modes can still be obtained under less restrictive conditions than the tightly wound approximation (Lau & Bertin 1978). The resulting dispersion relation is more complicated (it is a cubic rather than quadratic expression in k) and depends on a new dimensionless parameter J: where the last approximation holds for a Keplerian disk. While the Q parameter balances shear, thermal pressure and disk mass, J is a measure of the disk-star mass ratio. The cubic dispersion relation reduces to the standard Lin and Shu expression for the case of light disks (J or M d /M much smaller than unity), while for massive disks, where J ≈ 1, or M d /M ≈ 1/m) lower m, less tightly wound, spiral arms arise. In both regimes, the most unstable mode scales with H. When H is small, the behavior is well captured by a local, even 2D analysis. Otherwise, any manifestation of the instability arises over a significant fraction of the radial extent of the disk.

The onset of linear instability
When the disk becomes linearly unstable, the onset of GI is similar in both regimes of J. A spectral analysis of GI in numerical simulations confirms that the dominant modes excited have a wavenumber k ≈ 1/H, as predicted from the Lin & Shu dispersion relation in Equation (13) (Boley et al. 2007;Cossins et al. 2009;Michael et al. 2012). Similarly, the azimuthal wavenumber m scales inversely with H so that thick / massive disks are characterized by fewer spiral arms and by a more open spiral structure. We provide a more detailed description of the onset of linear instability in each regime below. Figure 1 illustrates the morphological differences between the instability with small and large H/r.
Tightly Wound, small H/r. The development of small scale instabilities can be captured by both local shearing box / shearing sheet simulations and by thin disk global simulations. This regime is commonly referred to as gravito-turbulence, although some caution should be used when talking about 'turbulence' in this context.  used a Fourier analysis of GI in the  shearing sheet approximation to show that the power spectrum of perturbations generated by the instability is contained in modes with wavelengths of a few H.  finds essentially no power at scales much smaller than H, which are resolved in the calculations. A similar result has recently been obtained by Shi & Chiang (2014) using three dimensional shearing box simulations. Young & Clarke (2015) also show -in the context of fragmentation (see below) -that collapse of unstable modes is never initiated on scales H. These results imply that 'gravito-turbulence' extends over a very small range of lengthscales. Additionally, there is no obvious indication of an energy cascade with dissipation occurring on the smallest scales. On the contrary, in global, thin disk simulations, Cossins et al. (2009) find that mode dissipation occurs through large scale, almost sonic shocks. In this respect, the dynamics introduced by local GI is not akin to what one usually calls turbulence. Morphologically, structures in shearing box simulations look more like turbulence than the high m, tightly wound spirals seen in global simulations. This might be an artifact of the shearing sheet/box approximation, but merits further investigation. Whether gravitoturbulence shares properties with isotropic, isothermal turbulence has important implications for fragmentation, as we discuss in Section 3.8.
Large H/r. As the mass of the disk is increased, the disk becomes unstable at larger values of H/r and higher values of Q, as predicted from the Lau & Bertin (1978) dispersion relation. Once the disk-star mass ratio increases above ≈ 0.1, the evolution of the instability changes character somewhat, transitioning from quasi-stationary spiral structures, to recurrent strong episodes of spiral activity, followed by brief quiescent phases. This behavior has been observed in a range of simulations including isolated disks with adiabatic equations of state with and without optically thin cooling Mejia et al. 2005;Laughlin et al. 1997), isothermal and irradiated embedded disks (Krumholz et al. 2007;Kratter et al. 2010a), and even two-dimensional, embedded disks with radiative heating and cooling (Zhu et al. 2012). A similar, recurrent behaviour has also been identified in N -body galaxy simulations, where dynamical cooling occurs via the injection of low-velocity dispersion, or 'cold' particles (Sellwood & Carlberg 1984).
The episodic behavior observed in global simulations is likely due to one of two phenomena. The first is the growth of the m = 1 mode, as the disk mass becomes comparable to the stellar mass (M d ≥ 0.3M * , in analytic estimates Adams et al. (1989)). Non-axisymmetric perturbations grow due to interaction with the"indirect potential" (Adams et al. 1989). The indirect potential, caused by the displacement of the central body from the system center of mass, manifests as an extra term in the gravitational potential when working in a reference frame centered on the central mass. These modes arise at higher values of Q than the axisymmetric modes. Simulations of massive disks that neglect this term may underestimate the growth of instabilities. Under some circumstances, SLING amplification may occur (Shu et al. 1990), where the m = 1 mode can amplify by interacting with the disk edge and the outer and inner Lindblad resonances. Because the m = 1 mode leads to very large mass transport rates, the disk mass drains on such short timescales so that it can (at least temporarily) become stable and quiescent. The second source of the quasi-periodic behavior is likely mode coupling, which we discuss in Section 3.5.3.

Accretion and transport
We now turn to the most important consequence of GI in disks: its ability to drive angular momentum transport, and therefore accretion of matter onto the central object. The fact that a spiral structure in the disk can transport angular momentum was recognized initially by Lynden- Bell & Kalnajs (1972) in the context of galactic dynamics. The Rφ component of the associated stress is given by where gr and g φ are the radial and azimuthal component of the perturbed self-gravity field and the brackets indicate azimuthal averaging. To this stress one should also add the induced Reynolds stress: where δvr and δv φ are the perturbed fluid velocities (see e.g. Balbus & Papaloizou 1999). For a magnetized disk (not considered here), one should also add the relevant Maxwell stress to the two terms described above. The Shakura & Sunyaev (1973) α prescription relates the stress T rφ to the local disk pressure by The dimensionless coefficient, α, often called the anomalous viscosity parameter, is expected to be less than unity, To facilitate comparison with an eddy viscosity, this is often expressed as implying that sub-sonic eddies on scales smaller than H facilitate angular momentum exchange. The introduction of this parameterization is especially useful for broad parameter studies in a range disks at moderate computational cost by solving a 1D diffusion equation, often coupled with an energy equation, as we discuss in Section 3.7. One can also compute numerically the stress induced by GI (or another process like the MRI) by running a (magneto-)hydrodynamical simulation of an unstable disk, and compare the resulting total stress to Equation (21) to obtain an effective α.
Let us now consider the torque and power associated with a spiral wave induced by GI. Standard wave mechanics (Toomre 1969;Shu 1970;Fan & Lou 1999) link the energy and angular momentum of a wave, Ew and Lw respectively, to the wave action A: where Ωp = ω/m is the pattern speed of the perturbation. Thus the ratio of energy and angular momentum (and the ratio of power and torque exerted on a disk when launching the wave) is the pattern speed Ωp, rather than the local Ω as for a viscous stress. More explicitly, the wave energy is given by where δΣ/Σ is the fractional amplitude of the perturbation. Equation (25) can be rewritten as the sum of two components: The second term in the above expression is equal to the wave angular momentum times the local angular velocity and thus represents a local, viscous-like term. The first term is an 'anomalous' energy transport term, which is inherently non-local (Balbus & Papaloizou 1999). The ratio of these two terms, |Ωp − Ω|/Ω, is therefore a dimensionless parameter that measures the degree of non-locality of the disturbance. In practice, the transport induced by spiral waves will be global in character if the waves are able to travel within the disk significantly away from their corotation radius.
We have seen already that the local dispersion relation involving the parameter J (see Equation (18) above) predicts that the transition between local and global behavior is related to the mass ratio M d /M . Is the mass ratio condition equivalent to a division of regimes based on the distance that waves propagate before dissipating? Clearly, to answer this question we have to look at the spectrum of modes excited by the instability. Using a spectral analysis of the dominant modes excited by GI, one finds a clear correspondence between morphology and the dominant mode of angular momentum transport. Cossins et al. (2009) computed the pattern speed of the dominant modes, finding that independent of the disk-star mass ratio, the Doppler shifted Mach number M of the dominant perturbation is remarkably close to unity: This implies that the shocks that develop as soon as the waves travel far enough from corotation to become sonic quickly dissipate the waves. Based on this sonic condition it is now possible to translate the 'locality' parameter defined after Equation (26), into the disk-star mass ratio: For a Q ≈ 1 disk, H/r ≈ M d /M , thus the 'corotation condition' and the 'disk mass' condition are equivalent, so long as waves dissipate once they become sonic. We immediately see that the fraction of transport induced by non-local processes is proportional to the mass ratio.

Saturation of the Instability vs Non-linear growth
Thus far we have only discussed the linear growth of GI. However, to properly assess the outcome of the instability we need to consider how it evolves in the non-linear regime. We have introduced two dimensionless parameters, Q and M d /M * that capture the linear phase. In order to understand the non-linear evolution, we require a third, independent dimensionless parameter, which captures the radiative timescale of the disk.
3.5.1. Thermal Saturation. The non-linear evolution of the instability and its saturation is best understood using numerical hydrodynamical simulations of self-gravitating disks, which have been the focus of most of the attention in this field over the last twenty years. It is still useful to begin with a few simple analytic arguments. Consider first an isolated disk with an initial surface density and temperature profile such that Q ∼ 1 over a range of radii. In this scenario, Q will evolve either due to mass redistribution or heating and cooling. We begin by examining changes in temperature rather than surface density because the viscous timescale is often slow compared to the thermal timescale. Since the excitation of waves will inevitably lead to shocks, it is natural to imagine a self-regulation mechanism (originally proposed by Paczyński 1978) whereby the parameter Q plays the role of an effective 'thermostat' for the disk temperature. In the absence of other significant heating terms, the disk cools down until it reaches marginal stability. At this point, GI sets in, causing shocks that heat the disk. If sufficient heat is retained, the disk will heat back up, quenching the instability, which in turn decreases the heating rate, allowing the disk to cool back towards instability. Thus in this self-regulated state, the disk should remain marginally unstable, with Q close to unity. In such a state the disk permits growing modes, but the spiral perturbations do not grow exponentially, instead saturating at some finite value (see Equation (31) below).
One can estimate analytically the conditions under which the disk should self-regulate. Since the instability provides heat on a timescale that is of order the dynamical timescale in the disk, t dyn = Ω −1 , the 'thermostat' can only work if cooling is slower than this timescale, that is if β = Ωt cool 1, where t cool is the cooling timescale 1 . Thus we see that β is a third, independent dimensionless number that controls not the onset of the instability like Q, or its morphology, like M d /M * , but its evolution and saturation.
The importance of the parameter β, and the relation to self-regulation, was first illustrated numerically in a paper by . This work introduced an explicit cooling term in a self-gravitating shearing sheet calculation where u is the internal energy of the fluid and t cool is the cooling timescale, which is set as a free parameter in the simulations to be some constant times the dynamical timescale Ω −1 . It is important to stress that the above description is not meant to reproduce any specific cooling law, but is a toy model for exploring the role of the cooling timescale in the outcome of GI. Different investigators use different prescriptions for the specific form of the cooling time, which in some cases is taken to be a constant (Mejia et al. 2005) while in other cases is taken to be proportional to the dynamical timescale as a function of radius, so that β = Ωt cool is kept constant (e.g., . In particular, the so-called β-cooling prescription is now commonplace when simulating self-gravitating disks in a simplified manner. One noteworthy limitation of this model is it can only capture the behavior of optically thin gas in numerical simulations, because it does not account for optical depth effects, which can be important. If one aims to capture more realistic thermal physics, radiative transfer must be properly accounted for (Boley et al. 2006;Stamatellos & Whitworth 2009b;Krumholz et al. 2009). Nevertheless, the evolution produced by simplified cooling models is remarkably similar to that in more complex simulations.
In general, there is now a coherent picture emerging from all such simulations of isolated, nonirradiated discs, despite the large variety of hydrodynamic techniques (SPH or grid based codes, in 2D or in 3D) and cooling prescriptions employed. If cooling is relatively inefficient, so that the cooling timescale is long (β 10 − 20, though see Section 3.8), the disk settles into a selfregulated state, where the spiral structure steadily transports angular momentum providing the heating required to balance cooling and keep the disk in thermal equilibrium. We call this thermal saturation of GI. We discuss in Section 3.8 uncertainties in the critical value of β and 2D versus 3D effects.
As discussed above, if the disk is low mass (with H/R ∼ M disk /M 0.1), the transport induced by GI is local, and analogous to a viscosity. If dissipation is the only heating source, then a disk undergoing GI in thermal equilibrium (Pringle 1981;, will have the following relation between α and the cooling rate: Since the energy density and the angular momentum density associated with density waves scales with the square of the perturbed gas density (e.g. Eq. 25), the condition for thermal regulation (Eq. 30) requires that the instability saturates at an amplitude that scales with the square root of the cooling time: a behavior confirmed by the numerical simulations of Cossins et al. (2009). This relation is consistent with the critical β being of order a few, since it is the point at which density perturbations have nonlinear amplitude. In practice, energy dissipation occurs through the development of roughly sonic shocks. This simple relationship suggests that the strength of the instability, and its efficiency as an angular momentum transport process, are correlated with the thermodynamics in steady state.
We review GI models that use α as a proxy for transport in Section 3.7.
3.5.2. Thermal saturation in irradiated disks. The discussion and the simulations described above assume that the only source of heating in the disk is that provided by the instability itself as it dissipates energy through shocks. However, protostellar disks are generally expected to be irradiated by either the central protostar, their natal envelope, or nearby stars. In the presence of irradiation, clearly the thermal saturation described above must change, since the disk alone is not responsible for setting the equilibrium temperature. The thermostat cannot function according to Equation (30). When comparing GI in irradiated versus self-heated disks, one must consider the behavior at fixed Q, so that the disk has the same propensity to the instability. Even at fixed Q, the impact of irradiation is not intuitive. On one hand, in the presence of external heating a much lower level of dissipation due to the instability is required to prevent runaway growth, which would imply a stabilization of the disk. On the other hand, a very strongly irradiated disk will behave more like an isothermal disk, where the dissipation of energy due to even a large perturbation has little effect on thermal balance (Kratter & Murray-Clay 2011). In this case it is hard to envision a mechanism that would stop an unstable perturbation from growing non-linearly, leading to fragmentation. Indeed irradiated disk simulations have more similarities to isothermal simulations than barotropic or cooling simulations Offner et al. 2009). For disks in the moderately irradiated regime, both arguments apply.  conducted local 2D simulations of disks including a 'standard' β-cooling term and a constant heating term, representing the effects of irradiation. On the basis of cooling times, irradiation stabilizes the disk in that a marginally faster cooling time is required for fragmentation in the irradiated case as compared to the non-irradiated case. However, Eq. 31 does not hold anymore, and for a given value of β the average perturbation amplitude is much smaller. Thus, in the irradiated case, even perturbations with a very small average amplitude result in a runaway collapse into fragmentation. In this sense, irradiation has the effect of destabilizing the disk, because smaller amplitude perturbations can collapse. It is unclear whether the use of 2D simulations in this study might affect the measured critical β (Young & Clarke 2015).
In irradiated disks, there is another route to saturation that does not relate to thermal balance, but rather is dynamical, in the form of non-linear coupling between modes.
3.5.3. Mode-mode coupling. The numerical investigation of self-gravitating disks began well before sophisticated radiative transfer cooling algorithms were computationally feasible. Early numerical simulations found that both isothermal (Laughlin & Bodenheimer 1994), and barotropic (Laughlin & Ròżyczka 1996;Laughlin et al. 1998) disks could develop strong gravitational instabilities that lead to fast redistribution of angular momentum within the disk. In these regimes, GI cannot self-regulate in the manner describe above as they represent β → 0 and β → ∞ respectively. Instead, non-linear mode coupling appears responsible for saturating the amplitude of GI-induced perturbations. Laughlin et al. (1997) demonstrate that coupling between m = 0 and m = 2 modes prevents runaway growth and thus fragmentation, while allowing rapid angular momentum transport. This same behavior is evident in the isothermal simulations of Kratter et al. (2010a), and is also likely responsible for the transient spiral patterns seen in the massive disk simulations of .
The exact conditions under which mode-mode coupling dominates remain unclear, as little detailed work has been done since Laughlin et al. (1997). Based on the simulations in which a similar effect has been observed, we expect this saturation mechanism to dominate in massive disks, as well as disks that are irradiation dominated. These two cases represent (1) the regime in which local dissipation breaks down and (2) the regime in which the energy balance due to GI induced heating and radiative cooling breaks down. Ironically, this saturation mode may be the dominant mechanism preventing fragmentation in protostellar disks, despite the dearth of attention in the literature.

When saturation fails: Fragmentation
If the saturation mechanisms listed above fail to prevent runaway growth of perturbations, the unstable region of the disk will break up into bound objects, or fragments. As the perturbations grow to non-linear amplitudes, their own self-gravity causes them to collapse and separate from the background disk. Fragmentation is in some ways an alternative means for shutting off the instability: the separation of large amplitude perturbations form the disk means a reduction in local surface density, and thus an increase in Q. In the context of protostellar disks this is of great importance, as fragmentation can in principle lead to the formation of either stellar, brown dwarf, of perhaps planetary mass companions. We provide a brief overview of the two common pathways to fragmentation, deferring a discussion of the astrophysical implication to Sections 4 and 5 Fragmentation driven by cooling. A great majority of disk simulations have considered disk fragmentation as consequence of rapid cooling in isolated disks. In this case, Q decreases due to a progressive drop in temperature, until GI sets in. If the cooling is rapid, GI cannot provide sufficient heating to counteract cooling, and the disk will fragment. Exactly how fast the disk must cool (the so-called critical value of β) has been the subject of many investigations. Studies have examined the role of the temperature dependence of cooling (Johnson & Gammie 2003;Cossins et al. 2010a), the role of the evolution of β in time at fixed radius , and the role of irradiation (Boley et al. 2007;Stamatellos & Whitworth 2008;. The critical value is a function of the equation of state, and is of order β = 10 − 20. The deviations from the critical value of β due to more complex physical effects are of order a factor of two. We review the numerical complications of measuring this value in Section 3.8. Cooling as a path to fragmentation is most likely relevant for disks that are thermally selfregulated because, as discussed above, irradiation and global transport can hamper the operation of the disk thermostat. AGN disks, or protostellar / protoplanetary disk which are shadowed by inner disk wall or warps may be the best candidates for this type of fragmentation. Alternatively disks which undergo rapid changes in accretion rate, which in turn can effect the stellar luminosity and thus the irradiated disk temperature, might also be subject to fragmentation driven by cooling (Kratter & Murray-Clay 2011).
Fragmentation driven by accretion. The second path to fragmentation is due to evolution of the disk surface density either due to infall from the external environment, or due to changes in the accretion rate within the disk that lead to mass pile-ups. For disks undergoing infall at some rateṀin, there exists a regulation mechanism akin to thermal saturation. Consider a disk with Q > 1 that accretes material such that Q declines towards unity. As the disk approaches Q ≈ 1, the linear instability transports angular momentum. If the accretion induced by GI givesṀ <Ṁin, the disk will become more massive, and thus susceptible to a wider variety of modes, enabling more transport. IfṀin ≤Ṁmax,GI, the maximum possible rate that GI can provide (as yet unspecified, see below), the disk will transport material at the rate it is being fed, with the disk mass acting as the regulator of transport (Kratter et al. 2010a;Zhu et al. 2012). IfṀin >Ṁmax,GI the disk will ultimately fragment (Kratter et al. 2008). The ability of the disk to regulate transport up to a critical rate is likely related to the spiral-mode coupling described above.
What sets the maximum transport rate induced by GI, and how does one know when disks will fragment versus increase in mass? Kratter et al. (2010a) argued that two dimensionless parameters, benchmarked to the infall rate, are strong predictors of fragmentation driven by accretion. These two are where M * d , the total system mass. The first parameter, ξ references the infall to the isothermal sphere collapse rate, c 3 s /G  but uses the disk, not core temperature. For a disk in steady state, ξ and α are interchangeable, using the steady-state accretion relationṀ = 3αc 3 s /(GQ) (Frank et al. 2002). However in many cases, the disk accretion rate does not come into equilibrium with the infall rate. The second parameter, Γ is the ratio of the orbital timescale to the disk mass doubling timescale. When infall is significant, Ω is controlled by the circularization radius of accreting material (not viscous spreading) and is thus a proxy for disk angular momentum.
Using a suite of 3D numerical simulations of isothermal disks undergoing self-similar accretion, Kratter et al. (2010a) find that disks undergoing infall with Γ < ξ 2.5 /850 are unstable to fragmentation. Subsequent work by Offner et al. (2010) and Zhu et al. (2012) confirm these relations hold when more complicated thermal physics and more realistic, turbulent infall conditions are included. Such a parameterization is applicable to AGN disks as well, where it might be relevant for seed black hole growth at high redshift (Lodato & Natarajan 2006).
The ξ, Γ parameterization predicts that the maximum transport rate is a function of mass ratio (in steady state with Q ∼ 1, we have that H/R ≈ M d /M * ≈ (Γ/ξ) 1/3 ). For low mass disks, where GI-induced transport is local,  have shown that the existence of a critical cooling rate for fragmentation translates, through Equation (30), into a more general upper limit of α ≈ 0.1. The boundary for fragmentation in ξ − Γ space is consistent with this value of α for equivalently low mass disks, despite the absence of thermal regulation. If the disk mass grows above 0.1−0.2M * , however, transport becomes inherently global, and thus an increase in mass ratio corresponds to an increase in the maximum rate of transport. This is consistent with the Lau & Bertin (1978) dispersion relation. According to the ξ − Γ boundary derived above, a more massive disk corresponds to a larger value of Γ for fixed ξ, which will thus be stable at the same infall rate. These simulations suggest that global modes can provide close to α ≈ 1 when M d → M * .
Other fragmentation triggers. Nelson (2000); Mayer et al. (2005); Boss (2006) and Lodato et al. (2007) have studied the influence of a massive companion on disk fragmentation, reaching apparently contradictory conclusions. While Boss (2006) argues that the presence of a companion increases the tendency for fragmentation, all of the other studies obtain the opposite result. They find that the tidal heating due to the companion has the effect of inhibiting fragmentation in a marginally stable disk.

Capturing GI with α models
As mentioned in the previous sections, it is often convenient to parametrize angular momentum transport measured in numerical simulations in terms of the Shakura & Sunyaev (1973) α-parameter. We have relied on this parameterization heavily in discussing thermal saturation (see Equation (30)). Despite the inherent assumption of local transport in the α-model (Balbus & Papaloizou 1999), the approximation is also often used for global transport, where it approximates the accretion rates measured remarkably well. Kratter et al. (2008) compiled transport measurements from numerous simulations, and found that a two component α model -one to account for local transport, one for global, could accommodate a wide range of disks in terms of both mass and thermal physics. As noted above, the biggest advantage of invoking an effective α to capture GI is that one can conduct parameter studies over a wide swath of physical variables using semi-analytic or 1D models. We provide a brief overview of the most common models used.
The first proposed αsg models of Lin & Pringle (1987 posit that since Q controls the onset of instability, the strength of transport should scale inversely with Q. They propose: where ζ is an additional parameter measuring the strength of the gravitationally induced torque, and Q is the maximum value (≈ 2) at which GI sets in. This prescription, at face value, is orthogonal to α models based on Equation (30). Although the dependencies appear distinct between Equation (30) and Equation (34), they can be tuned to capture similar behavior, but are better suited to different regimes. Equation (30) is appropriate for disks that are in steady-state and thermally self-regulated. It cannot be used to capture the behavior of a disk with Q much different from unity, or one that is out of equilibrium because the relation is predicated on energy balance. On the other hand, Equation (34) can capture GI in a disk that is non-self regulated and may be evolving (either in surface density or temperature) in time. Because Equation (34) has a free normalization parameter, it can also describe a disk in steady-state that is marginally unstable by tuning ζ so that the transport rate matches up with the R.H.S. of Equation (30).
These α models have proved indispensable for studying multiple transport processes simultaneously. Armitage et al. (2001) explored 1D viscous models including both MRI and GI-induced transport, showing time variable accretion capable of reproducing FU-Ori like episodic accretion. Bertin & Lodato (2001) discuss how GI would make the outer disk in FU Ori objects hotter compared to non self-gravitating disks and hence produce a flatter SED (see also Adams et al. 1989). Martin & Lubow (2011) similarly identify regions of disk parameter space lacking steady-state solutions as responsible for outbursting behavior, in a fashion akin to the limit-cycle instability for dwarf novae. Matzner & Levin (2005) first identified, through simple α models, the existence of a critical radius at ≈ 50 − 100 AU beyond which the disk is expected to fragment. Clarke (2009) use an α model for GI to estimate the relevant cooling time as a function of radius, and show that fragmentation occurs beyond ∼ 70 AU, while self-regulation is expected between ∼ 20 − 70 AU. Crucially, Clarke (2009) neglected the effect of irradiation from the central star, that in general makes the inner disk gravitationally stable rather than self-regulated.
Recently, Rafikov (2015), following Rafikov (2009), developed a more comprehensive viscous α prescription, which encompasses a wider range of disk states, although the caveats about locality of energy deposition remain. Rafikov (2015) finds that the predicted disk properties are quite similar to those derived under α(Q) models (e.g. Kratter et al. (2008); Zhu et al. (2010)). The important feature of these viscous models is the smooth growth of transport rates as Q declines, with a cutoff at some maximum α for a minimum Q.

Convergence of numerical results
In our analysis of the saturation of GI, we have relied heavily on numerical simulations. The two values often pulled from these simulations are the critical value of β, or the maximum transport rate α, at which saturation fails, leading to fragmentation. These two quantities (directly correlated in some regimes) are only measurable in simulations. Frustratingly, these quantities can be difficult to disentangle from numerical artifacts. No astrophysical simulations operate at realistic Reynolds numbers, and all suffer from numerical errors that manifest as an effective viscosity. In the case of SPH simulations, an artificial viscosity must be included explicitly. For grid codes, Krumholz et al. (2007) for example, has measured the effective α from diffusion as a function of resolution, showing that it can be competitive with physical mechanisms at moderate resolution.
At the time of this writing, the convergence problem in self-gravitating disks simulations is an open issue. We review them below for completeness, but emphasize that for application to protostellar and protoplanetary disks, uncertainty in the critical value of β has little impact on disk evolution. Indeed, in the region where protostellar disks are gravitationally unstable, the cooling time is a steep function of disk radius, so that changes in the critical value of β only lead to minor changes to the radius at which fragmentation occurs (Clarke 2009;).
Until recently, most investigations of resolution focused on resolving fragmentation. Early studies found that the Jeans length (Bate & Burkert 1997;Truelove et al. 1997) must be resolved by ≈ 10 grid cells to accurately capture fragmentation. More recent work has argued that the Jeans length must be resolved by as many as 64 grid cells to achieve convergence (Turk et al. 2012). Since in SPH the resolution is set by the smoothing length h, and the Jeans length is of the order of the disk thickness H, the Truelove and Bate & Burkert criteria imply that h/H 1 to resolve fragmentation. This condition is amply satisfied by most recent simulations. Nelson (2006) lists three different conditions: one is a variant of the Jeans length criterion, the second is that h/H 1, and the third is that adaptive smoothing lengths be used in SPH (which is done in all modern SPH calculations). That convergence remains a challenge even when the above resolution criteria are met is not so surprising, since none explicitly capture one of the most important aspects of the physics: dissipation. Meru & Bate (2011a,b) found that in simulations with several million SPH particles (factors of several above most previous work) the fragmentation criteria as measured by a critical β did not appear to converge. The disk fragmented at larger β with higher resolution. Similar results have come from grid-based simulations. Paardekooper et al. (2011) found that fragmentation could be artificially triggered at sharp boundaries between the stable and unstable portions of the disk, while Paardekooper (2012) claimed, based on 2D simulations, that fragmentation might be a stochastic process, and therefore allowed (rarely) at values of β associated with much smaller amplitude perturbations. This interpretation would be consistent with gravito-turbulence fully sampling a turbulent power spectrum of density perturbations (Hopkins & Christiansen 2013). These results cast doubt on the numerical value of the fragmentation boundary obtained in previous investigations (that pointed to a critical value βcrit ≈ 10).  address the importance of the magnitude of the artificial viscosity as a function of resolution, which inherently alters the disk thermostat. In SPH codes, artificial viscosity provides an equivalent αart ∝ h/H.  find that in order for artificial viscosity to be negligible as a source of heating, one requires h H 40 β .
Note that as β is increased (as one does when probing the fragmentation boundary), the resolution requirement becomes progressively harder to satisfy, because a smaller amount of numerical dissipation represents a larger fractional contribution to heating. Such a trend of increasing resolution requirements with increasing β is indeed consistent with the Meru & Bate (2011a,b) results, but do not fully explain them, as it turns out that, based on these simple arguments, fragmentation would be artificially quenched by numerical viscosity when the latter produces only 5% of the heat required by thermal equilibrium. Meru & Bate (2012), on the other hand, point out that most SPH simulations for self-gravitating disks have employed an artificial viscosity formulation that strongly reduces viscosity for high Mach number shocks. Paradoxically, this reduction of the viscosity coefficient has the counter-intuitive effect of increasing artificial dissipation, since it generates artificially large velocity gradients (Lodato & Price 2010). Even more extreme challenges to the fragmentation boundary have come from Paardekooper (2012) and Hopkins & Christiansen (2013) who suggest that fragmentation may be inherently stochastic, and thus allowed at much longer cooling times. This relies on GI producing a wide power spectrum of turbulent fluctuations Investigations by Young & Clarke (2015) question this for two reasons. First, they show that independent of numerical method, there are inherent limitations in the treatment of the small scale gravitational field in 2D. As a result, measurement of the fragmentation boundary (or stochastic fragmentation) in a 2D simulation is questionable, and better left to 3D simulations. They further argue that because of the quasi-regular nature of selfgravitating structures, stochastic fragmentation should be inhibited at very large β, because the growth and development of spiral structure occurs over a modest number of orbital periods, and thus may not be akin to a true turbulent cascade (see also Shi & Chiang (2014)).
Despite the lingering questions about resolution and critical cooling times, we emphasize that the critical β for fragmentation has very little impact on whether disk fragmentation, and especially planet formation via GI, is viable because most disks that have Q ∼ 1 also attain relatively small values of β for realistic protostellar disk conditions.

Application to Protostellar Disks
In Section 3 we have given an overview of the basic physics that governs the behavior of disks when self-gravity becomes important. These principles apply to disks in any astrophysical context. Here we narrow our focus to the case of protostellar and protoplanetary disks in order to understand which of the regimes and corresponding behaviors are most relevant. The details of how star formation proceeds inevitably effects the characteristics of protostellar disks, and thus the importance of self-gravity and GI. It is beyond the scope of this review to describe how different modes of star formation might produce different disk properties; instead we endevaor to provide a broad overview of all plausible disk states.
Narrowing the Parameter Space. We have shown that GI can be described by three independent dimensionless numbers that describe the state of the disk, Q, M d /M * and β. However, a protostellar or protoplanetary disk does not attain arbitrary combinations of these three quantities because the disk's internal energy and cooling rate are both temperature dependent. Consider β in terms of these physical variables: where τ = κΣ/2, and κ is the opacity. The cooling time, and thus β, depend on the disk surface density, temperature, and dust opacity. Once one selects a disk model with a given Q profile and mass ratio, the profile β(r) is also determined (albeit with uncertainty due to dust opacity models). We present typical disk scalings for temperature, mass, size, and accretion rate in order to survey GI in protostellar and protoplanetary disks, highlighting three case studies.

Setting the Thermal Physics
We begin with a simple model for the temperature profile of the disk. Ideally disk temperatures should be found from accurate radiative transfer calculations, which account for wavelength dependent opacities and heating sources. Nevertheless, the agreement between the simple analytic formulae below and radiative transfer calculations is quite good, making these useful for estimation purposes (Zhu et al. 2012;Boley et al. 2010). Disks are heated primarily by three sources: accretion energy, host star irradiation, and external radiation. The disk midplane temperature can be found via the following relation by balancing heating and cooling terms: where F * and Fext are the energy fluxes associated with stellar irradiation and external heating respectively, while Facc is the energy flux associated with accretion and (Rafikov 2005). Here, τ is the Rosseland mean opacity, and f (τ ) reasonably captures how the accretion energy diffuses outward from the midplane in both the optically thick and thin regimes. We now examine each term in detail. The first term on the R.H.S, Facc, describes the release of gravitational energy as material moves through the disk, falling deeper into the potential well. This is commonly referred to as viscous dissipation. Independent of whether the disk is truly viscous, the potential energy of accreted material must be released, some portion of which will go into the thermal energy of the disk. The flux (far from the stellar surface) associated with this term if all of the accretion energy goes into heating is: Because the heat must propogate from the midplane to the surface, this term is modified by the optical depth. As discussed in Section 2, disk opacities are uncertain, and likely evolve in time due to grain growth and changes in the dust-to-gas ratio. At very early times, when grains are ISM-like (µm − 10µm), the opacity is temperature dependent, typically κR ≈ κ2T 2 (where T is in Kelvin), (Bell & Lin 1994), in the cold, icy limit T < 155K relevant for GI. Semenov et al. (2003) find that the Rosseland mean opacity for small grains is well fit by κ2 ≈ 5 × 10 −4 cm 2 /g. Once grain growth commences, grains quickly become large compared to the blackbody wavelengths of the disk. In this regime, the opacity is both lower and temperature independent. Based on Pollack et al. (1985) we use a temperature independent κ0 = 0.24cm 2 /g. We benchmark our results to this large grain case for several reasons. First, after no more than 10 5 yrs, there is good evidence that grain growth has commenced. Although GI may be present earlier, at this epoch the disk is optically thick for either opacity law. Secondly, lower opacities associated with grain growth produce a wider range of disks susceptible to GI; since our intent is to review all the possible phase space we consider this choice "conservative." The second term on the R.H.S. of Equation (37) captures irradiation from the host star. The stellar luminosity comprises two components: accretion luminosity and the intrinsic stellar luminosity that comes either from gravitational radiation or nuclear burning depending on age and mass. For low mass stars the accretion luminosity typically dominates. High mass stars, which reach the ZAMS while accreting, are dominating by H-burning (Tout et al. 1996). Unlike the dissipation of accretion energy, the impact of stellar irradiation on disk temperatures is highly dependent on the envelope and disk structure, because the absorption of radiation depends on both the intervening material reprocessing the radiation and on the flaring angle of the disk absorbing radiation. We consider two models for stellar irradiation. The first, from ray tracing calculations of Matzner & Levin (2005), finds that for deeply embedded disks, the heating term is well described by: with f ≈ 0.1 derived from ray tracing calculations where much of the stellar flux impacts the infalling envelope (ala Terebey et al. (1984)) and is reprocessed back down onto the disk. Note that L * can include both internal luminosity and stellar accretion luminosity. The second model refers to the case of a non-embedded disk, for which the analytic model of Chiang & Goldreich (1997) accounts self-consistently for the flaring due to incident irradiation, and heating of dust grains in the outer disk layers. They find: where αF measures the grazing angle at which starlight hits the disk; this in turn is a function of the height of the disk photosphere, which depends on grain properties. When disks reach the threshold for instability, they are typically optically thick, making the use of this approximation reasonable. We can rewrite this explicitly as a function of stellar luminosity and disk radius as: We note that when using this relation for temperature due to stellar irradiation in conjunction with other heating terms, there is a small inconsistency in that the actual disk scale height will differ slightly from that solved for in the irradiation model.
Including both heating terms, we see that an individual disk may span both the self-luminous and irradiated regime at different radii. It is instructive to examine the radial scaling of the first and second terms in Equation (37): where the last two lines cover the embedded and isolated cases. In either case, the viscous term will fall off more steeply with radius than either irradiation term. Thus the inner disk is more likely to be self-luminous, while outer disk may be dominated by irradiation from either the central star or even the external environment. Finally, external radiation from other stars can contribute to the disk temperature (see Thompson 2013 for an extreme example). We consider the interstellar radiation field to create a temperature bath of Text (such that Fext = σT 4 ext ), which sets a lower threshold for the disk temperature as described above. In more active, higher mass star forming regions this can be up to Text = 20−30K (Caselli & Myers 1995). For deeply embedded cores, the envelope may insulate the disk from external sources, providing a somewhat lower temperature thermal bath. Neither irradiation term is accompanied by an optical depth because they set the temperature at the surface, which controls the cooling (L.H.S. of Equation (37)).

Disk Masses and Radii
In Section 2 we showed evidence for disks ranging from M d /M * ≈ 0.001 − 0.1, but argued that in many cases these may be lower limits based on uncertain grain growth, dust-gas conversions, and biases against observing the youngest, most massive disks. Here we consider disks with 0.05 < M d /M * < 0.5 to encompass both the least massive disks where self-gravity can be relevant, and a generous upper limit on those observed. Recalling Equation (2), one can immediately see that a Q = 1 disk at the upper boundary has an aspect ratio so large as to call into question the disk geometry.
Direct observations show that disks extend out to 50 − 500au, based on submm and optical absorption of proplyds (Andrews et al. 2009;Eisner et al. 2008). Disk radii at early times (Class 0, I phases) are controlled by the infalling angular momentum from the envelope or feeding filament (Terebey et al. 1984). Most of the disk mass will be contained within the circularization radius: where j is the average specific angular momentum of accreting material, which can be obtained from consideringJ (t) =Ṁ (t)j(t).
In a simple core-collapse model one can specify j(t), however it is more realistic to consider the impact of accretion from turbulent filamentary structures that characterize star forming regions (Offner et al. 2010). Turbulence in the interstellar medium follows a line-width size relation, such that increasingly larger scales will typically carry more angular momentum (Larson 1981). Kratter & Matzner (2006) used this to generate models for disk radii as a function of time. It is not guaranteed, however, that the direction of the angular momentum vector will remain constant throughout infall, and thus disks may not grow monotonically in time. Recent simulations by Bate et al. (2010) and Fielding et al. (2015) show that there can be enough change in the infalling angular momentum vector to tilt the disk with respect to the spin-axis of the star, especially if variable accretion causes substantial oscillations in the disk mass, and if stellar interactions truncate the disk. When disks transition from the protostellar to protoplanetary phases, infall of mass and angular momentum becomes less important. In this case one expects the disk to evolve towards the selfsimilar viscous spreading solution of Lynden- Bell & Pringle (1974). While this will lead to larger outer disk radii than that given by the circularization radius, the outer regions typically hold very little mass and thus do not have a significant impact on disk self-gravity.
In order to determine Σ(r), we require a model for the radial profile. Here again we rely on observations, which suggest a power law, or a power law plus exponential tail (Andrews et al. 2009;Pérez et al. 2012). Since the exponential part contains an irrelevant amount of the mass, we consider only the power law component. For our calculations, we choose Σ ∝ r −1 , which puts more mass at large radii where instabilities might arise, lest we underestimate the role of self-gravity. Note that this surface density profile, coupled with an irradiation dominated temperature profile is consistent with steady-state disk models withṀ (r) =const. To give a quantitative idea of the scaling for our parameter ranges:

Accretion Rates
The only remaining parameter to specify is the accretion rate through the disk, which has two effects. First, as discussed in Section 3.6, infall of material on to the disk impacts the disk mass, angular momentum transport rate, and susceptibility to fragmentation. Secondly, the disk temperature is a function of the accretion rate. We consider disks with steady-state accretion at rates ranging fromṀ = 10 −8 − 10 −3 M yr −1 , where the lower end is consistent with accretion rates observed on actively accreting T Tauri stars towards the end of their lives, and the highest rates represent a generous upper limit set by infall from a supersonic core for more massive stars (McKee & Tan 2003;Klaassen et al. 2011) .
In reality, a disk with a given surface density and temperature profile cannot support an arbitrary accretion rate in steady state. A physical mechanism for providing angular momentum transport is required. Indeed there is no guarantee that a disk will typically be in steady-state witḣ M (r) = const (Armitage et al. 2001;Martin & Lubow 2011).
As discussed in Section 3, saturated GIs provide transport at rates such that 10 −2 < α < 1 depending on the disk parameters. Obviously magnetic effects, both small scale turbulence from the MRI (Balbus & Hawley 1994;Simon et al. 2013), or large scale magnetic winds (Bai 2015) may also be important for various disk properties at different times.

The instantaneous disk state
We demonstrate how GI manifests in protostellar and protoplanetary disks with three case studies shown in Figure 2 that cover a range of parameters. Note we have not yet specified how disks might evolve into and out of these parameters, and some areas of the parameters space are unrealistic. For each case study we provide two contour plots of fundamental disk parameters as a function of disk radius and accretion rate. We fix the disk mass, outermost radius and host star properties, and then calculate the implied properties of the disk at all locations for a range of accretion rates.
The left hand panels show contours of the disk temperature with all three heating terms from Equation (37). For comparison we also show in white contours the temperatures that would be derived neglecting stellar irradiation, as this has been omitted in many numerical simulations. The accretion rate not only contributes to the viscous heating term, but also the stellar heating term, where for low mass stars we use the models of Baraffe et al. (2015) to estimate the stellar radius at 1Myr. The right hand panels show contours of Q and β (calculated from Equation (36)). When relevant, we demarcate the boundary where the accretion rate assumed would require α > 1. This region is inconsistent with GI as a transport mechanism. 2 We describe each of the three cases in detail below: This disk is most consistent with observations of Class II protoplanetary disks. It is immediately clear that these disks will not be subject to GI at any radius. Q is lowest at the outer edge, where cooling times are also very rapid, and the temperature is controlled by stellar irradation. In order to make the outer regions of this disk unstable to GI, the disk mass must be increased by roughly a factor of 4. In this case the cooling time remains short, suggesting that such disks would fragment, but never enter the self-regulated regime, as first recognized by Rafikov (2005); Matzner & Levin (2005).
Disk 2:. Figure 2b shows a disk modelled after observations of IRAS 16293-2422b (Rodríguez et al. 2005). This disk is most likely associated with an isolated Class 0 source. The disk appears to be massive and compact, consistent with a very young age, with M d /M * = 0.35, Rout = 30au around a 0.8M protostar. The Rodríguez et al. (2005) model suggests an accretion rate of order 10 −6 M /yr, but we show the full range as for the other examples. Because the disk has low Q and is slowly cooling, we expect one of two outcomes depending on the rate at which mass is accreted onto the disk. Up to some critical rate, the disk can process the material via GI-induced transport and remain quasi-stable. This would require a combination of thermal regulation and mode-mode coupling. At larger infall rates, the disk will begin to grow in mass, and either revert to the first state of self-regulation, or increase in mass until it fragments.
Disk 3. Figure 2c shows a disk model appropriate for a massive star, such as that identified by Shepherd et al. (2001), consistent with theoretical models of high mass star formation (Kratter & Matzner 2006;Krumholz et al. 2007). Here M d /M * = 0.75, Rout = 150au and M * = 8M . Unlike the low mass stars in the previous two examples, high mass stars are expected and observed (Klaassen et al. 2011) to have higher infall rates. Therefore the upper half of the plot is more representative of likely parameters. Here we see that the entire outer logarithmic decade in radius will have Q ∼ 1. While at distances of the order of tens of au the cooling time is relatively long and the disk might survive in a self-regulated state, beyond ∼ 50au the cooling time becomes short enough to induce fragmentation. The disk temperatures are higher in part because massive stars arrive on the ZAMS and start fusion while still accreting, but they still have low Q values. The high mass ratio represented in this example may be typical because the high infall rates tend to drive up disk masses. In this case, simulations suggest that disk undergo strong m = 1 dominated transport, and in many cases fragmentation into stellar companions (Adams et al. 1989;Kratter et al. 2010a). The susceptibility of these disks to fragmentation is consistent with the observation that nearly all massive stars are found in binaries (Sana & Evans 2010;Sana et al. 2014).

Global Models for Disk Evolution
We have laid out the instantaneous properties of protostellar and protoplanetary disk in the previous sections. However, protostellar disks do not exist in isolation, but are constantly fed material from the background core or filament. We argued in Section 3 that infall can control the linear and non-linear development of GI. There, we simply imagined that a given disk, and considered what happened as we changed the infall rate. In order to understand what properties disks have during their evolution, star formation models are required that dictateṀin andJin. The state of the disk is controlled by competition between the infall rate and the disk-star accretion rate. Several authors have considered this global evolution: Kratter et al. (2008), who used analytic models tied to numerical benchmarks for non-linear transport terms, Vorobyov & Basu (2007, who consider 2D models for the disk considering only transport driven by GI, Zhu et al. (2010Zhu et al. ( , 2012, who consider long timescale semi-analytic models, and shorter timescale 2D simulations. These calculations suggest that low mass stars, as demonstrated by disk case studies (1) and (2), are susceptible to GI at early times, but stable against fragmentation. In order to push these disks into the fragmenting regime, one must invoke highly variable accretion or lower stellar illumination. Highly variable accretion can force disks like case study (2) up to the fragmentation threshold, where binary formation is the most likely outcome. Reduced irradiation is more relevant for older disks, where fragmentation into low mass objects is most likely as we discuss in the following section. These studies also suggest that burst-like behavior is probable when GI is coupled with other modes of angular momentum transport like the MRI. Finally, Kratter et al. (2008), who explored the evolution of higher mass stars, found that more massive stars often evolve from the self-luminous global mode regime at early times into the irradiated, massive, and fragmenting regime at late times. The primary difference is that more massive stars are thought to experience higher mass infall rates, which more easily exceed the maximum GI accretion rate discussed in Section 3.6.
To summarize, we have shown that protostellar and protoplanetary disks sample several of the regimes discussed in the previous section. Compact, Class 0 disks (see, for example Rodríguez et al. 2005), can be GI-unstable, self-luminous, and regulated via thermal saturation and modemode coupling. For most other cases, a disk that reaches the threshold for gravitational instability will have a temperature controlled primarily by stellar irradiation, with short cooling times. Under these circumstances, the disk will either fragment, or the instability will saturate through modemode coupling, depending on whether or not the infall rate exceeds that which can be processed by strong global modes. As a result, it is imperative to consider the impact of both stellar irradiation and infall in characterizing the behavior of these disks.

GI and planet formation
The non-linear phase of GI can lead to fragmentation of the disk into bound objects in orbit about the primary star. What are these objects? Beginning with Adams et al. (1989), there is a growing consensus that GI-borne objects are most likely to evolve into either brown-dwarfs or stellar binary companions around sufficiently massive stars (Kratter & Matzner 2006;Kratter et al. 2010b;Stamatellos & Whitworth 2009b;Zhu et al. 2012;Forgan & Rice 2013). In contrast, Boss (1997) argued that the direct collapse of clumps from gravitational instabilities can create massive gas giant planets (see also  and references therein). More recently, a series of papers starting with Nayakshin (2010), have proposed that a combination of GI and tidal stripping can even produce rocky planets. Considering which disks are susceptible to GI and fragmentation reveals parameters unfavorable to the formation of planets, but favorable to the formation of more massive wide-orbit companions. In Section 4 we have seen that unstable disks are massive, which leads to large initial masses for fragments, often undergoing rapid infall, which can enhance the subsequent growth rates, and experiencing rapid angular momentum transport, which can cause rapid migration and destruction of marginally bound objects (which are the most likely to be low in mass in the first place). All of these features favor the formation of objects above the deuterium burning limit.
The disks best suited to forming lower mass objects are those that are unstable in the absence of continuing infall. These disks tend be less massive, and are thus colder when Q = 1, which translates to smaller initial masses, less continued growth, and the possibility of slower migration in a more laminar disk. While these conditions are much more favorable for forming low mass objects, the formation of close-in companions < 20au, with low mass < 1MJup, is the exception.

Initial clump conditions
The mass at which a GI perturbation becomes a bound fragment is the outcome of the non-linear phase of the instability, and is thus challenging to predict from simple analytic models. Although numerical simulations are now in principle well enough resolved to answer these questions, the masses evolve from this initial state quickly, and simulations become poorly resolved on a similar timescale. It is thus valuable to have analytic estimates with which to compare numerical models.
The simplest approximation is that: is the most unstable wavelength for the axisymmetric instability when Q = 1. The factor of 4 comes from assuming that only 1/2 of the unstable wavelength generates the over density. As noted by Kratter & Murray-Clay (2011), such a fragment would by definition be unbound at the Σ associated with Q = 1, because the size implied by λ is larger than the Hill radius for the mass M f in Equation (49). This is consistent with the instability producing non-linear perturbations. Many authors have argued for different, order unity factors in front of the dimensionally motivated Σλ 2 estimate (e.g. Boley et al. 2010;Forgan & Rice 2011). In Figure 5.1 we illustrate several different fragment mass estimates, which span an order of magnitude in mass. To compare these mass estimates we assume a temperature profile:

What is a planet?
To determine whether gravitational instability makes planets, one requires a working definition of a planet. At present, there is no consensus in the community. In this chapter, we use nuclear burning as the dividing line, although fragmentation is ignorant of this boundary. We also suggest that system architecture be considered in separating planetary systems from stellar multiples. While planets typical orbit a dominant central mass (or occassionally a tight binary), higher order stellar systems typically have hierarchical orbits.
which is (conservatively) a ∼ 20% reduction from the fiducial irradiation model of Kratter et al. (2010b). We set the disk surface density such that Q = 1 at every radius shown. We choose a critical value of β = 13, (relevant only for the prediction from Forgan & Rice 2011). In addition to the initial fragment mass estimates, we also show several limiting cases. We plot the canonical gap opening mass (using α set by the cooling time) and the isolation mass (Lissauer 1987). These scales are described in more detail in Section 5.1.3. Even the lowest mass estimate predicts initial fragments above 1MJup for typical disk properties. The estimates of the mass at which fragments might cease to grow are all well above the deuterium burning limit. Clearly fragments must be prevented from growing if they are to become planets.
5.1.1. Clump Survival. Regardless of the mass at which clumps are born, they face a variety of threats to long term survival. These can be categorized as: 1. Dispersal due to interaction with another over density (clump or spiral arm) 2. Tidal shredding due to slow cooling 3. Tidal shredding due to rapid migration All three of these rely on the balance between contraction due to cooling, which allows the fragment's gravity to become ever more dominant, and disruption due to tidal forces from either the star or another overdensity in the disk.
First consider the cooling rate of a newborn fragment. In Section 3 we discussed the critical role of the cooling time in determining whether an unstable disk can thermally self-regulate before the instability becomes non-linear. In modelling clump survival, we must consider the cooling properties of the clump itself as it separates from the disk and increases in density. Although the numerical values of the critical cooling time are not wildly different between disk and fragment, they can in principle be different if there are sharp changes in the equation of state with density. Moreover, the appropriate radiative cooling timescale for a clump is, in fact, different from that which governs thermal self-regulation. The former is a comparison of the disk's thermal energy to the radiative cooling rate, assuming a vertically isothermal disk. The cooling timescale relevant for clump survival describes the evolution of a midplane perturbation. The appendix of Kratter et al. (2010b) works out this timescale as: where T0 is the temperature set by irradiation alone and κ ∝ T B is the dependence of the opacity on temperature, typically 0 < B < 2. Comparison with Equation (36) shows that the perturbation cooling timescale is only different by a factor of a few at most. The cooling rate sets the size of Here we illustrate the range of initial fragment masses with lower and upper bounds set by Boley et al. (2010) and Forgan & Rice (2011), respectively. We also show estimates for the masses at which accretion onto the fragment from the background disk might cease: the gap opening mass and the isolation mass.
See text for details about the disk model used for this calculation.
a fragment as a function of time. For typical fragmenting disks, which have β 10, cooling is sufficiently fast that fragments initially contract on their free fall timescale (neglecting tidal forces).
The contracting fragment's radius must be compared with the scale on which tidal forces can shred it. Naively, one requires that: where a is the current semi-major axis. Consider the origin of this limit: the Hill radius RH determines where an object's self-gravity will dominate over tidal gravity. In the case of a collapsing fragment, pressure forces may not be negligible, weakening the inward pull of self-gravity. Kratter & Murray-Clay (2011) find that tidal shredding can occur at radii of order R Hill /3 if pressure support is important. However, even if fragments suddenly slow their cooling, and undergo Kelvin-Helmholtz contraction rather than free fall collapse, most will contract to < R Hill /3 within a dynamical time, almost guaranteeing survival at a fixed radius.
Clump interaction. Fragments may also be disrupted due to tidal interaction with other nearby fragments or spiral arms. Shlosman & Begelman (1989) define a critical cooling time for clumps to collapse prior to interacting. This is relevant if the disk breaks up into multiple fragments with separations of order their own size. Since the encounter velocity of two fragments is roughly the Hill velocity vH = GM clump /RH , fragments cross each others Hill radius in one dynamical time, implying that fragments must collapse to well within the Hill radius on this very rapid timescale. This stringent requirement of collapse in one orbital period is unlikely to be necessary in fragmenting protostellar disks, which have H/r sufficiently large that fragments form relatively far apart (Boley et al. 2010;Kratter et al. 2010a, though see Stamatellos & Whitworth 2009a).
In addition to interacting with other clumps, clumps may also interact with a passing spiral wave, since we know these should be present when fragmentation begins. Young & Clarke (2016) argued that the quasi-regular nature of the spiral patterns in self-gravitating discs (as first characterized by Cossins et al. 2009) should encourage proto-fragment disruption via collisions with spiral features over many orbital timescales. Their quantification of the intermittency of spiral features in simulations suggested that 'stochastic' fragmentation should not be possible for β > a few ×10.
5.1.2. The physics of clump cooling. We have outlined above some disruption mechanisms and the critical dimensionless cooling rates. But what determines the cooling rate physically? The cooling physics is nearly identical to that at the onset of star formation (Larson 1969), though self-gravitating clumps will generally be optically thick from the start (Rafikov 2005). At low temperatures, as in the background disk, cooling rates are set primarily by dust opacity. If grain growth has already commenced in such disks the opacity might be somewhat reduced compared to the ISM case. As the clump contracts under its own self-gravity, compressional heating causes the central temperatures and densities to rise. At temperatures above roughly 1200K (Bell & Lin 1994;Semenov et al. 2003), dust sublimation removes the dominant source of opacity, accelerating cooling somewhat. Once temperatures rise all the way to 2000K, triggering molecular hydrogen dissociation, rapid collapse ensues toward planetary sizes. Galvagni et al. (2012) has carried out a careful study of the cooling timescale of clumps including a complex equation of state. They extract a clump from Boley et al. (2010) (neglecting subsequent growth, which can affect both the mass and cooling on similar timescales). They include radiative cooling due to a mixture of atomic and molecular hydrogen, with an opacity appropriate for ISM grains. They find that clumps that do not have their rotation artificially suppressed collapse to roughly 40% of their initial radius within 10 dynamical times. The same clump with rotation suppressed collapses even faster to roughly 15% of the initial radius over 10 dynamical times. This suggests that the clumps described by these models would satisfy the constraints given above, predicting survival at fixed radius.

Continued Growth.
A detailed treatment of the evolution of clumps based on their initial mass provides a useful benchmark, but omits the influence of continued accretion from the background disk. The absolute upper limit is set by the isolation mass, which is the mass that a secondary potential embedded in the disk can attain by accreting disk material within its own Hill sphere (Lissauer 1987). The importance of continued accretion following fragment formation is likely a strong function of initial mass. If the initial mass is very large, it will clear a deep gap in the surrounding disk, possibly reducing the accretion rate. In contrast, low mass fragments, or massive fragments in rapidly accreting disks, would remain semi-embedded in the natal disk.
Gaps are opened in a disk when a relatively massive companion produces gravitational torques which repel material faster than viscosity can replenish it. Using the standard gap opening require-ments of Lin & Papaloizou (1986) and Bryden et al. (1999) where fg = 0.23 is a geometric factor derived by Lin & Papaloizou (1993). By comparing Equation (54), the red line, with the initial mass estimates shown in Figure 5.1, we see that typical fragments will be near, but below the gap opening mass, and likely able to grow at least this much. We benchmark our gap-opening estimate with a relatively large value of α as this is typical for strongly self-gravitating disks.
Gap-opening, which occurs at masses well below the isolation mass, may not entirely starve growth. Numerous simulations of Type II planet migration and gap opening show that at least some material does flow across the gap, and enter into circumplanetary disks. Indeed Duffell et al. (2014) find that the standard Type II migration rates that assume no mass flow are incorrect for the same reason. Lissauer et al. (2009) found in numerical simulations that gaps roughly 5RH in width are capable of entirely shutting off accretion. Whether this limit, obtained in much lower viscosity disks, applies, is uncertain.

Migration of clumps.
Even if clumps survive at fixed radius, as discussed above, they are subject to tidal torqueing by the parent disk. Fragments may initially be just below the gap opening mass, suggesting that they begin migrating similarly to planets undergoing type I migration. Baruteau et al. (2011) have confirmed in 2D simulations that the migration of such objects in gravitationally unstable disks is quite rapid, and that the gap opening criteria is accurate even in such disks. They find that despite the presence of stochastic outward forcing due to GI, the overall migration direction is inwards. Similarly, Zhu et al. (2012) find that planets suffer divergent fates depending on whether or not they open a gap; gap-opening fragments first move in, then stall their migration, while smaller planets migrate inwards on < 10 dynamical times. This is consistent with the work of Vorobyov (2013) in terms of migration timescale. The scenarios in which migration is most likely to stall or even reverse direction are when fragments become massive compared to the disk out of which they formed. While very interesting for understanding the distribution of brown dwarf and stellar binary companions (Kratter et al. 2010a), these are not relevant to the question of giant planet formation. 5.1.5. Long timescale disruption: Tidal downsizing. As noted at the start of this section, a second GI planet formation paradigm has arisen in the last few years wherein an initially massive fragment becomes enriched in solids and migrates inwards until its Hill Radius shrinks below the fragment radius. The envelope of the fragment is completely or partially stripped, leaving behind a low mass, rocky core with or without an atmosphere (Nayakshin 2010;Boley et al. 2011). Successive refinements of this model include the effect of radiative feedback (Nayakshin & Cha 2013), heavy element sedimentation (Nayakshin et al. 2014), and late stage metal enrichment (Nayakshin 2015) 3 . One significant concern with these models is the attempt to disentangle the initial disk conditions from the clump evolution. While very sensible from a modeling perspective, this can have an outsize impact on the conclusions. For example Nayakshin & Cha (2013) use artificially inserted fragments in a non-fragmenting disk as the initial conditions, which does not adequately capture the resultant disk-fragment interactions. Similarly, Nayakshin et al. (2014) neglect accretion from the disk onto the fragment, which might fundamentally change the conclusions. These issues aside, we show in the following section that populations synthesis models of this process produce objects which are inconsistent with the rocky planet population known thus far.

Comparison with Observed Population
What, if any, insight can be gained from looking at the known population of planets and low mass companions? The GI model makes some predictions for the populations. The two simplest predictions come directly from the required initial conditions for clump formation: more massive stars seem more prone to developing gravitational instabilities and when clumps survive, they are prone to growing beyond the deuterium burning limit. Forgan & Rice (2013) have carried out the most extensive population synthesis calculations to date, including a range of initial masses, accretion of solids, migration, and tidal disruption. Notably, this work excludes any gas accretion after the initial clump is formed, meaning that their final masses are a lower limit (and much of the evolution might change if mass growth occurred early). They find that roughly half of the initially formed clumps are totally destroyed, and 90% of the remaining objects are above the Deuterium burning limit for a wide variety of initial models. A strong prediction of these models is that GI mostly makes BDs or more massive objects when it occurs, and even if tidal truncation / downsizing can in principal produce some rocky cores, it is a finely tuned part of the parameter space and unlikely to be responsible for the bulk of observed rocky planets.
In Figure 4, we illustrate the population of low mass companions as compiled by the exoplanets.eu catalog. This catalog includes companions up to 30MJup. We show companions as a function of mass ratio, rather than measured mass, as this has more relevance to the question of GI. Note that most of the Kepler planets are excluded from this diagram as their masses are unknown. Based on their periods and radii, they primarily occupy the region within 1AU of the host star and below 10 −5 in mass ratio. GI is most likely to produce the objects at 10s of au above mass ratios of 10 −2.5 . In order for GI to account for more of the "normal" planet population this outer region would likely be more populated than allowed for by current observational limits (Forgan & Rice 2013).
On the opposite end of the mass spectrum, is their observational evidence for GI as a mechanism for BD and binary formation? Yes and no. On the one hand, the propensity of massive stars to binary formation via disk fragmentation is generally consistent with the abundance of binaries around more massive stars (Raghavan et al. 2010;Sana & Evans 2010). On the contrary, Kraus et al. (2011) find that for young stars in Taurus, the lower mass objects M * < 0.7M typically have companions with smaller separations consistent with disk fragmentation (< 100AU ,) whereas stars with 0.7M < M * < 2.5M have more companions at larger separations. Recent work on Class 0 and I binaries now finds evidence for a bimodal distribution of separations, characteristic of a disk and core mode of fragmentation (Tobin et al. 2016). Dupuy & Liu (2011) found that the orbital distribution of BD binaries was inconsistent with ejection in pairs from fragmenting disks as predicted by Stamatellos & Whitworth (2009a). As of this writing there is not a complete model for binary formation from cores and disks. Bate (2012) provides the largest statistical sample of simulated binary formation, which matches well with observations of the field. Unfortunately the disk fragmentation process is not sufficiently well resolved in these calculations to attribute parts of the binary distribution to one mechanism. 5.3. Role of GI in the core accretion model for planet formation GI may have an entirely separate means of kickstarting the planet formation process: by concentrating grains in spiral arms. Protoplanetary disks can be well-modelled as two fluids: gas and 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 Semi-Major Axis (AU) The distribution of planetary and other low mass companions from Exoplanets.eu as a function of mass ratio versus semi-major axis. All objects are shown in green, those above the deuterium burning limit are demarcated by red dots, those with more massive primaries are shown in yellow, while those with M-star hosts are green solids. As noted in Section 3, while the gas feels pressure forces, which counteract gravity and slow orbital velocities, solids act as a pressureless fluid. The relative velocities between the two fluids generates drag, which is a function of particle size. This drag force results in slow radial drift of particles towards pressure maxima (Whipple 1972;Adachi et al. 1976;Weidenschilling 1977). We neglect in this review the separate physics of self-gravity of solid layer in protoplanetary disks (Goldreich & Ward 1973;Weidenschilling 1980), or self-gravity of planetesimals concentrated by the streaming instability (Youdin & Goodman 2005). Spiral arms in a self-gravitating disk are an excellent candidate for producing pressure maxima that can enhance particle growth, much like MRI turbulence for the Streaming Instability (Johansen et al. 2009). Rice et al. (2004) using 3D, global SPH simulations of gas and dust coupled via drag with a simple algorithm first demonstrated numerically that these traps work effectively in selfgravitating disks. Later, Rice et al. (2006), including the dust self-gravity, showed that the dust sub-structure induced by the traps can become self-gravitating itself, and lead to a rapid formation of large mass rocky cores. The actual mass of these cores is related to the Jeans mass of the dust substructure, which in turn depends on the velocity dispersion induced by the gas flow in the particle distribution. Whether the velocity dispersion measured in such simulations (typically of the order of the sound speed) is physical or affected by numerics is still to be fully understood (see Booth et al. 2015).  have raised the concern that this effect might require the cooling timescale in the disc to be relatively small based on the magnitude of the radial velocity induced by an overdensity of a given size. This would confine the effect to the outer regions of the disk. However, Gibbons et al. (2012), using 2D grid based simulations conclude that significant clumping occurs even for t cool ≈ 40Ω −1 , much above the fragmentation boundary. The particles with stopping times τs of order unity will typically be 1-10cm in size for the densities associated with self-gravitating disks at radii of 10-100AU. Gibbons et al. (2014) also consider the self-gravity of the particles and their back reaction on the gas. Again, similar to the Streaming Instability (Youdin & Goodman 2005), they find that particle concentrations in spiral arms can become high enough to become self-gravitating and collapse directly into large planetesimals.
A global view of particle clumping due to GI with an appropriate size distribution for the dust would be required to fully characterize the effect that GI might have in the core accretion model. This requires the use of an accurate and versatile algorithm for simulating the coupled dust and gas dynamics (Laibe & Price 2012, 2014Price & Laibe 2015).
Another issue to be fully understood is the fate of the rocky cores or planetesimals formed in a self-gravitating disk as described above. Generally, the mass of the rocky objects formed by particle concentrations within spiral arms is not expected to be high enough to directly trigger core accretion and significant further growth is required. However, such growth is significantly hampered during the self-gravitating disk phase, where the gas structure (even in the absence of drag) induces gravitational potential fluctuations that pump the eccentricity of the planetesimal swarm up to very high values, thus inhibiting significant growth (Britsch et al. 2008). Such an early population of planetesimals can thus produce large rocky cores only if they are retained in the disk long enough that the disk has become GI stable (Walmswell et al. 2013).

Conclusions and Outlook
In this review, we have provided a broad overview of the physical nature of GI, illustrating the variety of ways in which the instability manifests itself. It is first and foremost an angular momentum transport mechanism, and under some circumstances it may cause disk fragmentation, which can lead to the formation of secondary companions. We demonstrate that protostellar disks around low mass stars may be prone to a self-regulated form of the instability at early times, if their mass ratios are high. At later times and larger radii, disk are irradiation dominated and prone to fragmentation if unstable. The resultant fragments are unlikely to be below the deuterium burning limit. Protoplanetary disks are not expected to be unstable based on current observations. Disks around more massive stars are more likely prone to both instability and fragmentation because rapid infall drives the disks to higher masses.
There are a handful of theoretical issues that remain unresolved. The first set of open questions relate to saturation via mode coupling. Can mode coupling operate in tandem with thermal regulation? Under what circumstance does mode-mode coupling dominate over thermal regulation? Is there a critical disk-star mass ratio at which a transition occurs? How is the maximum mode amplitude set? Does mode coupling change under rapid infall? Since this form of saturation may dominate in protostellar disks, a detailed understanding of the mechanism is warranted. The second set of open questions relates to the critical cooling time, the nature of so-called gravito-turbulence, and stochastic fragmentation. Does GI produce a turbulent energy cascade and power spectrum akin to isotropic turbulence in the ISM? Or is gravito-turbulence observed in simulations the high resolution manifestation of small-m spiral mode growth? If GI does lead to turbulence, then the stochastic fragmentation at a range of cooling times is viable. If, on the other hand, the instability does not produce rare, high amplitude density perturbations, there should be a concrete value of the cooling time at which self-regulated disks transition to fragmentation On the observational front, much progress has been made in recent years, and ALMA is poised to make significant breakthroughs as all dishes come online. Protostellar disks are now being imaged with high enough resolution to reveal their internal structure and morphology with unprecedented details. As mentioned above, in several cases a prominent spiral structure can be observed (for example, in the NIR scattered light images of MWC 758, Benisty et al. 2015). Several papers have anticipated the exceptional ability of ALMA to detect spiral structures induced by GI both through molecular line emission (Krumholz et al. 2007) and through dust emission (Cossins et al. 2010b;Dipierro et al. 2014). It will be possible to distinguish the large scale, open spirals expected in massive disks (Krumholz et al. 2007) from the more tightly wound spirals expected in the low mass case (Dipierro et al. 2014). ALMA might even be able to detect the spatial changes in opacity due to dust trapping in spiral arms (Dipierro et al. 2015).

SUMMARY POINTS
1. Gravitational Instability provides efficient angular momentum transport in relatively massive protostellar disks. Evidence of this phase observationally remains sparse, but is likely to improve in the coming years with more ALMA data. 2. Based on the models presented in this review, we find that protostellar disks which achieve Q ∼ 1 often have large disk-star mass ratios and temperatures controlled by stellar irradiation. In this regime, Q evolves due to changes in Σ rather than T , and therefore infall is the dominant driver of instability. 3. When fragmentation does occur, BD or stellar companion formation is more likely than planet formation. Evidence that the instability ever forms planets is currently lacking. Uncertainties in the fragmentation boundary observed in numerical simulations have minimal impact on the likelihood of planet formation via this mechanism.

FUTURE ISSUES
1. Demonstration of the relevance of the instability to realistic disks requires high resolution, high sensitivity observations with ALMA coupled with sophisticated molecular line and radiative transfer models that can more reliably ascertain disk masses. 2. In protostellar disks, perhaps the main channel by which gravitational instability saturates -mode-mode coupling -remains poorly understood and should be the subject of future investigation. 3. Gravito-turbulence has not convincingly been shown to behave like turbulence with a multiscale power spectrum and energy cascade. The extent to which gravito-turbulence exhibits these characteristics is intimately tied to the issue of a critical cooling time and stochastic fragmentation. Stochastic fragmentation requires a very well sampled PDF, which may or may not occur in realistic disks.