1 INTRODUCTION
It is likely that self-gravity will, at times, play an important role in the evolution of accretion discs. In particular, we expect it to be important in discs around active galactic nuclei (Nayakshin, Cuadra, & Springel Reference Nayakshin, Cuadra and Springel2007) and in discs around very young protostars (Lin & Pringle Reference Lin and Pringle1987). Although we will discuss some general properties of self-gravitating discs, we will focus more on protostellar discs, than on other types of potentially self-gravitating accretion discs.
An accretion disc will be self-gravitating, or susceptible to the growth of the gravitational instability, if the Toomre parameter, Q, is close to unity (Safronov Reference Safronov1960; Toomre Reference Toomre1964)
$$\begin{equation}
Q = \frac{c_{\rm s} \Omega }{\pi G \Sigma } \sim 1,
\end{equation}$$
where c s is the sound speed, Ω is the epicyclic frequency (equal to the orbital angular frequency in Keplerian discs), Σ is the disc’s surface density, and G is the gravitational constant. This parameter, however, only tells us if the disc’s self-gravity is important and if it will likely develop a gravitational instability. It doesn’t, however, tell is how such discs will actually evolve.
There are two likely pathways for a self-gravitating disc. It could settle into a quasi-steady, self-regulated state (Paczynski Reference Paczynski1978), or it could become so unstable that it fragments into bound objects (Boss Reference Boss1998). The latter could be a potential star-formation process in discs around active galactic nuclei (Levin & Beloborodov Reference Levin and Beloborodov2003), or a mechanism for forming planets, or brown dwarfs, in discs around young stars (Boss Reference Boss2000). Even if a disc does not fragment, a self-gravitating phase could play a very important role in driving angular momentum transport and, hence, allowing mass to accrete onto the central protostar (Lin & Pringle Reference Lin and Pringle1987; Rice, Mayo, & Armitage Reference Rice, Mayo and Armitage2010).
In this paper, we will summarise our current understanding of the evolution of self-gravitating accretion discs. There are a number of factors that likely determine the evolution of such discs. In particular, the balance between heating and cooling, and the mass of the disc relative to the mass of the central object. We’ll discuss if we can establish the conditions which determine whether a self-gravitating disc will settle into a quasi-steady state, or fragment, and—if so—when such conditions are valid, and when they might be violated. We’ll look at the role of external factors, such as external irradiation and perturbations from stellar, or sub-stellar, companions. We’ll also discuss, in the context of protostellar discs, whether or not we expect fragmentation to actually occur and—if so—where it is likely to occur and what we might expect the outcome to be.
The paper is structured in the following way. In Section 2, we focus on the quasi-steady evolution of self-gravitating discs, how such discs evolve under the influence of external irradiation, what happens if the disc mass is high relative to the mass of the central object, and what might happen if we consider more realistic scenarios. In Section 3, we focus on the fragmentation of self-gravitating discs. In particular, we consider what conditions determine if a disc will fragment, when such conditions may apply and when they may be violated, could it be stochastic rather than precisely determined, and could fragmentation be triggered by external perturbations. In Section 4, we then discuss these results and draw conclusions.
2 QUASI-STEADY EVOLUTION
As already mentioned, a self-gravitating disc is one that has sufficient mass, or that is cold enough, so that the Toomre parameter, Q, is of order unity. Such a self-gravitating disc will generate an instability that will dissipate and heat the disc. If the cooling rate can match the rate at which the instability is heating the disc, the system may then sustain a state of marginal stability (Paczynski Reference Paczynski1978). If not, the system may fragment to form bound objects, potentially protoplanets or brown dwarfs in discs around young stars (Boss Reference Boss2000), or stars in discs around active galactic nuclei (Levin & Beloborodov Reference Levin and Beloborodov2003; Bonnell & Rice Reference Bonnell and Rice2008).
If such a system does maintain a state of marginal stability, one would expect the instability to transport angular momentum. A standard way to describe angular momentum transport in discs is to assume that it is viscous and can be described using the α-prescription (Shakura & Sunyaev Reference Shakura and Sunyaev1973). Energy transport in a self-gravitating disc, however, is not diffusive and so doesn’t, strictly speaking, satisfy this condition (Balbus & Papaloizou Reference Balbus and Papaloizoum1999). Numerical studies (Laughlin & Rozyczka Reference Laughlin and Rozcyczka1996; Nelson, Benz, & Ruzmaikina Reference Nelson, Benz and Ruzmaikina2000; Pickett et al. Reference Pickett, Cassen, Durisen and Link2000) have, however, shown that the α-prescription can be a reasonable way to parameterise angular momentum transport in self-gravitating discs.
One of the first studies to quantify this (Gammie Reference Gammie2001) used a local shearing sheet in which the cooling time, τ, was assumed to be a constant, β, divided by the local angular frequency, Ω:
$$\begin{equation}
\tau = \beta \Omega ^{-1}.
\end{equation}$$
Such a cooling formalism has now become known as β-cooling. What this work showed was that, as expected, if the disc settled into a quasi-steady state, the dissipation rate would match the imposed cooling, and the effective viscous α would then be given by
$$\begin{equation}
\alpha = \frac{4}{9 \gamma (\gamma - 1) \beta },
\end{equation}$$
where γ is the specific heat ratio. Figure 1 shows an example of the quasi-steady, surface density structure in a local shearing sheet simulation in which the imposed cooling is assumed to have β = 10 and in which γ = 2. Figure 2 shows the resulting α evolution. The simulation is initially run with no cooling, which is then turned on at t = 50. There is an initial burst and then the system settles into a quasi-steady state with an approximately constant α. The dashed line shows the α value to which we’d expect it to settle, based on Equation (3).
A snapshot showing the surface density structure in a shearing sheet simulation with β = 10 and α = 2. The disc has settled into a quasi-steady. self-regulated state in which heating via dissipation of the gravitational instability is balancing the imposed cooling.
The effective α plotted against time for the simulation shown in Figure 1. Cooling is turned on at t = 50, after which there is a sudden burst. However, by t = 100, the simulation has settled into a quasi-steady state with a roughly constant α, consistent with what is expected from energy balance (dashed line).
Other studies employing different types of numerical methods, such as Smoothed particle hydrodynamics (SPH, Lodato & Rice Reference Lodato and Rice2004; Rice, Lodato, & Armitage Reference Rice, Lodato and Armitage2005) and grid-based methods (Mejía et al. Reference Mejía, Durisen, Pickett and Cai2005), have similarly shown that a self-gravitating disc will tend to settle to a quasi-steady state in which Q is close to unity and in which angular momentum transport can be described using an α-formalism in which α is given by Equation (3).
There are, however, some caveats to the above that we will discuss in more detail later. It applies in situations where the disc mass is low enough that the local approximation is valid, and also requires that the cooling time is not so short that the instability becomes non-linear and the disc breaks up into fragments.
2.1. Irradiation
A fundamental assumption in some of the early work on the evolution of self-gravitating discs (e.g., Gammie Reference Gammie2001) was that there was only a single heating source; dissipation of the gravitational instabilty. In such a scenario, the cooling balances only this heating source and the effective α is given simply by Equation (3). In the presence of an additional heating source (such as external irradiation), it is more complex, since—in a state of marginal stability— the cooling balances both heating due to the instability and the additional heating source.
For a given cooling time, the effective α is reduced in the presence of an additional heating source, compared to a situation where the only heating source is the gravitational instability itself. It can be shown Rice et al. (Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011) that, in the presence of external irradiation, one can approximate the effective α using
$$\begin{equation}
\alpha = \frac{4}{9 \gamma (\gamma - 1) \beta } \left(1 - \frac{\langle \Sigma \rangle c^2_{{\rm so}}}{\langle \Sigma c_{\rm s}^2 \rangle } \right),
\end{equation}$$
where Σ is the local surface density, c so is the sound speed to which the disc would settle in the presence of external irradiation only, and c s is the sound speed to which the disc settles if it does tend to a state of marginal stability. This is similar to the form suggested by Kratter, Murray-Clay, & Youdin (Reference Kratter, Murray-Clay and Youdin2010). This can then be recast as
$$\begin{equation}
\alpha = \frac{4}{9 \gamma (\gamma - 1) \beta } \left(1 - \frac{Q^2_{\rm irr}}{Q^2_{\rm sat}} \right),
\end{equation}$$
where Q irr is the Q value to which it would settle in the presence of external irradiation only, while Q sat is the saturated Q value to which the disc settles in the presence of both the external irradiation and heating through dissipating the gravitational instability. Using shearing sheet simulations, Rice et al. (Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011) showed that this relationship is indeed reasonable, and the results are illustrated in Figure 3 (from Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011) which shows how the effective α varies with increasing levels of external irradiation, when the underlying cooling timescale is β = 9.
Figure showing how the quasi-steady α value varies with the level external irradiation (taken from Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011). The triangles are an initial semi-analytic estimate, while the diamonds with error bars are from the shearing sheet simulations. Although there is some discrepancy, it’s clear that the effective α decreases with increasing levels of external irradiation.
The important point here is that if there are other sources of heating, the strength of the instability is not set by the cooling timescale, β, alone. In the presence of other heating sources, the strength of the instability, and the magnitude of α, will be smaller than if the only heating source is the gravitational instability.
2.2. Local vs non-local
Another assumption underlying the idea that angular momentum transport via the gravitational instability can be treated as essentially viscous, is that the shear stress (or, equivalently, the viscous α) depends primarily on the local conditions in the disc. Gravity is, however, inherently global, and so such an assumption isn’t strictly correct. Balbus & Papaloizou (Reference Balbus and Papaloizoum1999) have shown that the energy equation for a self-gravitating disc has extra terms associated with global wave transport. They do, however, suggest that self-gravitating discs that hover near marginal gravitational stability may behave in a manner that is well described by a local α model.
Early global simulations (Laughlin & Bodenheimer Reference Laughlin and Bodenheimer1994) seemed to indicate that simple α models could reproduce the density evolution of self-gravitating accretion discs. Such simulations, however, initially ignored heating and cooling processes in the disc, which play a fundamental role in setting the quasi-steady nature of such systems (Gammie Reference Gammie2001). Imposing a β-cooling formalism in global disc simulations, Lodato & Rice (Reference Lodato and Rice2004) illustrated that quasi-steady, self-gravitating discs do indeed behave like α discs and that the local model is a reasonable approximation.
The local nature of self-gravitating angular momentum transport does, however, depend on the mass of the disc relative to the mass of the central object. As the disc-to-star mass ratio increases, the disc becomes thicker and global effects are likely to become more important. Lodato & Rice (Reference Lodato and Rice2004) only considered discs that had masses less than one quarter the mass of the central star. Later work (Lodato & Rice Reference Lodato and Rice2005) considered much higher mass ratios. What was found was that as the disc-to-star mass ratio increases, lower order spiral modes start to dominate. For mass ratios of ~ 0.5, a transient m = 2 spiral develops that rapidly redistributes mass and then allows the disc settle into a quasi-steady, self-regulated state. For mass ratios approaching unity, however, the disc never seems to settle into some kind of quasi-steady, self-regulated state. However, it does appear to satisfy the basic criteria in a time-averaged sense.
The basic picture therefore appears to be that for mass ratios below ~ 0.5 a self-gravitating disc will rapidly settle into a quasi-steady, self-regulated state in which Q ~ 1, heating and cooling are in balance, and in which angular momentum transport is pseudo-viscous and can be described by an effective α. As the mass ratio increases, lower m-modes begin to dominate, producing global spiral density waves that either rapidly redistribute mass, allowing the disc to then settle into a quasi-steady, self-regulated state, or—for sufficiently high mass ratios—that have time varying amplitudes leading to the disc only being quasi-steady in a time averaged sense.
2.3. Realistic cooling
A great deal of work on the evolution of self-gravitating accretion discs has used the basic β-formalism when implementing cooling. This, however, is clearly not a realistic representation of how such discs will actually cool, and was never intended to be so. It is simply a basic way in which to represent cooling so as to investigate how such discs evolve under a range of different cooling scenarios.
It’s clear that in reality, a self-gravitating disc will not have a cooling function that can be described via a single β parameter (Pickett et al. Reference Pickett, Mejía, Durisen, Cassen, Berry and Link2003) and that the actual evolution may depend on the form of the cooling function (Mejía et al. Reference Mejía, Durisen, Pickett and Cai2005). A number of studies have indeed used more realistic cooling formalisms. For example, Boss (Reference Boss2001, Reference Boss2003) and Mayer et al. (Reference Mayer, Lufkin, Quinn and Wadsley2007) used simulations with radiative transfer to understand the fragmentation of self-gravitating discs, Cai et al. (Reference Cai, Durisen, Michael, Boley, Mejía, Pickett and D’Alessio2006) considered how the disc metallicity might influence its evolution, Boley et al. (Reference Boley, Mejía, Durisen, Cai, Pickett and D’Alessio2006) considered the quasi-steady evolution of discs with realistic cooling and realistic opacities, and Meru & Bate (Reference Meru and Bate2010) probed the parameter space where disc fragmentation may occur, using three-dimensional hydrodynamical simulations with radiative transfer.
More recently, Forgan et al. (Reference Forgan, Rice, Cossins and Lodato2011) also considered the evolution of self-gravitating discs using a more realistic cooling formalism, but with the goal of trying to understand its quasi-steady evolution and, in particular, if the local approximation still applies when the cooling is more realistic. They performed three-dimensional SPH simulations of discs around young stars in which a radiative transfer formalism (Stamatellos et al. Reference Stamatellos, Whitworth, Bisbas and Goodwin2007; Forgan et al. Reference Forgan, Rice, Stamatellos and Whitworth2009) was used to estimate how the disc cooled.
Their results suggest that for discs with masses below about half that of the central star, the general behaviour is well described by the local approximation. An estimate of the local cooling time, when Q ~ 1, can be used to determine the local value of α. However, as the mass ratio increases, the behaviour becomes more global, and the local cooling time is not a good indicator of the effective α, except in a time-averaged sense.
Given the above, it is possible to use semi-analytic calculations to determine the realistic structure of quasi-steady, self-gravitating accretion discs (Rafikov Reference Rafikov2009; Clarke Reference Clarke2009; Rice & Armitage Reference Rice and Armitage2009), at least for those with masses less than about half that of the central star. For example, for a given opacity one can determine the cooling rate for a given surface density and, hence, the local effective viscosity (α). Using standard one-dimensional accretion disc models (Pringle Reference Pringle1981), one can then determine the surface density profile for a disc accreting at a specified accretion rate.
For example, Figure 4 shows the surface density profile of a quasi-steady, self-gravitating disc around a 1 M⊙ star accreting at ~ 10−7 M⊙ yr−1 (solid line) and accreting at ~ 10−8 M⊙ yr−1 (dashed line), for standard opacities (Bell & Lin Reference Bell and Lin1994) and in the absence of external irradiation. Given the opacities, the profile is essentially unique for a specific mass accretion rate. What Figure 4 also shows is that there is quite a strong mass dependence. Assuming an outer disc radius of 100 AU, the disc mass of the disc accreting at ~ 10−7 M⊙ yr−1 is ~ 0.35 M⊙, while that of the disc accreting at ~ 10−8 M⊙ yr−1 is ~ 0.25 M⊙. Less than a factor of two change in disc mass, produces about an order of magnitude change in accretion rate.
Examples of two surface density profiles for self-gravitating discs that are in a quasi-steady state. In these examples, there is no external irradiation and standard opacities are assumed. For given opacities, the profile is essentially unique for a specific mass accretion rate. The mass accretion rates here are ~ 10−7 M⊙ yr−1 (solid line) and ~ 10−8 M⊙ yr−1 (dashed line) and the disc masses are ~ 0.35 M⊙ (solid line) and ~ 0.25 M⊙ (dashed line). This illustrates that a relatively small change in disc mass can produce a substantial change in accretion rate, indicating that disc self-gravity is only likely to play an important role in angular momentum transport when these systems are very young, and the disc mass is relatively high.
This has a number of consequences. To explain TTauri-like, and higher, accretion rates, self-gravitating discs need to be quite massive, relative to the mass of the central star. Given that disc lifetimes are typically 5 Myr, or less (Haisch, Lada, & Lada Reference Haisch, Lada and Lada2001), self-gravity is unlikely to explain accretion rates during the later stages of a discs lifetime, since it would imply much longer lifetimes than those observed. Similarly, observed spiral structures in protostellar discs are unlikely to be due to disc self-gravity unless the disc mass is still relatively high (Dong et al. Reference Dong, Hall, Rice and Chiang2015). For discs with masses less than ~ 0.1 M⊙, the instability is likely to be weak, is unlikely to generate the stresses that can dominate the angular momentum transport process, and is unlikely to have spiral density waves that are sufficiently strong so as to be observable.
3 FRAGMENTATION
The previous section focussed on the quasi-steady evolution of self-gravitating accretion discs. For sufficiently unstable discs, however, one possible outcome is that they fragment to form bound objects. It has been suggested that these could be protoplanets in discs around young stars (Boss Reference Boss1998), or protostars in discs around supermassive black holes (Levin & Beloborodov Reference Levin and Beloborodov2003).
Again using shearing sheet simulations, Gammie (Reference Gammie2001) showed that a disc will fragment if the cooling time is sufficiently fast; in this case, the fragmentation boundary being at β = 3. This is illustrated in Figure 5 which shows the surface density structure in a shearing sheet simulation with β = 2 and γ = 2. It is clear that there are a number of dense clumps forming in the disc.
A snapshot showing the surface density structure in a shearing sheet simulation with β = 2 and γ = 2. There are a number of dense clumps, indicating that the rapid cooling in this simulation has led to the disc fragmenting, rather than settling into a quasi-steady, self-regulated state.
A similar fragmentation boundary was obtained using three-dimensional, global simulations (Rice et al. Reference Rice, Armitage, Bate and Bonnell2003). However, as Equation (3) shows, the effective α depends on the specific heat capacity γ. Consequently, it was shown that the fragmentation boundary is actually set by a maximum α that can be sustained in a quasi-steady, self-regulated disc, not by the cooling time specifically (Rice et al. Reference Rice, Lodato and Armitage2005). For example, the cooling time for fragmentation if γ = 2 is smaller than for γ = 1.4.
Additionally, one has to be careful because all of this early work was done under the assumption that dissipation of the gravitational instability was the only heating source. As discussed above, the presence of external irradiation weakens the instability and reduces the effective α (Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011). Therefore, it is more reasonable to regard the fragmentation boundary as being associated with the maximum self-gravitating α that a disc can sustain in a self-regulated, quasi-steady state, than as the minimum cooling time that can be imposed without the disc fragmenting.
Furthermore, the above analyses have all used simulations in which the local approximation is reasonably valid. As discussed above, as the disc becomes more massive (relative to the mass of the central object) global modes start to dominate, and the local approximations are only valid in a time averaged sense. When global modes start to dominate, it appears that a self-gravitating disc can sustain larger α values without fragmenting, than would be expected based on the local approximation only (Lodato & Rice Reference Lodato and Rice2005). Simulations of collapsing cloud cores with radiative transfer also indicate that the behaviour is quite different to what would be expected based on the local approximation (Forgan et al. Reference Forgan, Rice, Cossins and Lodato2011; Tsukamoto et al. Reference Tsukamoto, Takahashi, Machida and Inutsuka2015). Discs that don’t fragment can have α values considerably higher than the local approximation would suggest is possible. As Forgan et al. (Reference Forgan, Rice, Cossins and Lodato2011) illustrates, however, the Reynolds stress, driven by velocity shear from material from the infalling envelope contacting the disc, dominates, and the gravitational component of the stress still typically produces α values less than 0.1.
3.1. Perturbation amplitudes
If we focus again on the case where the local approximation is valid, one way to understand fragmentation is to consider the relationship between α and the perturbation amplitudes. Cossins, Lodato, & Clarke (Reference Clarke2009) show that the transport of energy and angular momentum, by the spiral waves driven by the gravitational instability, depends on the surface density perturbation δΣ. Since the angular momentum transport can be described via an effective α, which itself depends on the cooling timescale β, Cossins et al. (Reference Cossins, Lodato and Clarke2009) show that the perturbation amplitudes depend on β through
$$\begin{equation}
\frac{\left< \delta \Sigma \right>}{\left< \Sigma \right>} \simeq \frac{1.0}{\sqrt{\beta }}.
\end{equation}$$
Similarly, you can cast this as (Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011)
$$\begin{equation}
\frac{\left< \delta \Sigma \right>}{\left< \Sigma \right>} \simeq \sqrt{\alpha }.
\end{equation}$$
Figure 6 (taken from Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011) shows a plot of time-averaged α against < δΣ > / < Σ > illustrating that as α increases the perturbation amplitudes increase. Figure 6 also shows results from simulations with no external irradiation (Q irr = 0) and with external irradiation (Q irr ≠ 0) illustrating that it is α that determines the perturbation amplitudes, not simply the cooling time β.
A figure showing time-averaged α plotted against the mean perturbation amplitude (< δΣ > / < Σ >), for a variety of different shearing sheet simulations. It’s clear that there is a relationship between α and the perturbation amplitudes, with the perturbation amplitudes increasing with increasing α. In the absence of external irradiation (Q irr = 0), this would also indicate a relationship between the perturbation amplitudes and β (Cossins et al. Reference Cossins, Lodato and Clarke2009), but in the presence of external irradiation (Q irr ≠ 0) it is α that sets the perturbation amplitude, not simply β. (Figure from Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011.)
Essentially, the basic picture is that as the self-gravitating α increases, the perturbations become larger and, there will be an α value beyond which the perturbations collapse and the disc fragments. Most simulations suggest that this occurs at an α value of about 0.06 (Rice et al. Reference Rice, Lodato and Armitage2005), but this should probably not be seen as some kind of hard boundary, as we’ll discuss in more detail later.
3.2. Convergence
Recently, it has been suggested that the fragmentation boundary does not converge to a fixed α value as the numerical resolution of the simulation increases (Meru & Bate Reference Meru and Bate2011, Reference Meru and Bate2012; Baehr & Klahr Reference Baehr and Klahr2015). There have been a number of attempts to explain this. Lodato & Clarke (Reference Lodato and Clarke2011) and Meru & Bate (Reference Meru and Bate2012) suggest that numerical viscosity may play a larger role than previously thought (e.g., Murray Reference Murray1996) and that, consequently, simulations will need to consider higher resolutions than have been used in the past. Paardekooper, Baruteau, & Meru (Reference Paardekooper, Baruteau and Meru2011) show that the lack of convergence could be related to edge effects in simulations with very smooth initial conditions. This apparent non-convergence has, however, lead to a suggestion that fragmentation may occur for very long cooling times. Meru & Bate (Reference Meru and Bate2012) do, however, extrapolate their results to suggest that they may be tending towards convergence and that the critical value may be around β ~ 30.
One issue with the idea that fragmentation could occur for very long cooling times is, however, based on our basic understanding of the underlying processes. As discussed above, in the absence of other heating sources, there is a relationship between cooling time and α, such that as the cooling time gets longer, the effective α value gets smaller. However, this also implies that the perturbation amplitudes get smaller. So, one issue would be how do discs fragment if the perturbations are tiny?
Also, if fragmentation at longer cooling times requires higher resolution, then it implies that this fragmentation is occurring on scales that are not resolved by the lower resolution simulations. However, the Jeans mass/radius at Q ~ 1 is typically resolved in at least part of the simulation domains for these lower resolution simulations (Rice et al. Reference Rice, Paardekooper, Forgan and Armitage2014; Tsukamoto et al. Reference Tsukamoto, Takahashi, Machida and Inutsuka2015). This would suggest that at these longer cooling times, large amplitude, small-scale perturbations are driving Q ≪ 1 and hence producing fragments that are unresolved in the lower resolution simulations. In shearing sheet simulations, however, the power typically decreasing with increasing wavenumber, k, as shown in Figure 7. The issue then becomes why at long cooling times fragmentation happens at scales where there is much less power than at the scales at which fragmentation occurs when the cooling time is short?
Figure showing the power spectrum of the perturbations in a shearing sheet simulation with β = 10. The length of each side of the sheet is L, and so this figure shows that most of the power is at scales a few times smaller than the sheet size and that there is very little power at very small scales.
Consequently, Rice, Forgan, & Armitage (Reference Rice, Forgan and Armitage2012), suggest that the effect is largely numerical and is related to the manner in which cooling has typically been implemented in SPH simulations. More recently, an SPH simulation with 10 million particles and with a modified form of the standard β cooling formalism, was shown not to fragment for β = 10 (Rice et al. Reference Rice, Paardekooper, Forgan and Armitage2014), well below the values suggested by Meru & Bate (Reference Meru and Bate2011) and Meru & Bate (Reference Meru and Bate2012). Michael et al. (Reference Michael, Steiman-Cameron, Durisen and Boley2012) also considered convergence in a three-dimensional grid-based code and, although they didn’t test for the fragmentation boundary specifically, they did find that simulations with α < 0.06 typically did not fragment. Hence, it seems that it’s important to be careful of how the cooling is implemented and to avoid implementations that may promote fragmentation.
There has been some attempt to try and see if there is a physical explanation for this lack of convergence. In particular, is there evidence that simulations are suppressing fragmentation at small scales? One possibility, suggested by Young & Clarke (Reference Young and Clarke2015), is that there are two modes of fragmentation. One mode operates at the free-fall time and requires rapid cooling. The other is a mode in which fragments form via quasi-static collapse and requires that these fragments then contract to a size where they can survive disruption by spiral shock waves. Young & Clarke (Reference Young and Clarke2015) suggest that many two-dimensional simulations smooth gravity on the Jeans scale, inhibiting contraction to smaller scales, and hence preventing this second mode from operating. However, they also show that even if this second mode does operate, it probably still requires relatively fast cooling (β ≲ 12), but with quite a large uncertainty.
Part of the motivation behind Young & Clarke (Reference Young and Clarke2015) was a suggestion that turbulent discs are never stable against fragmentation (Hopkins & Christiansen Reference Hopkins and Christiansen2013) and, in particular, that rare high-amplitude, small-scale fluctuations could lead to gravitational collapse. This suggestion, however, is not entirely consistent with what we’re discussing here. The Hopkins & Christiansen (Reference Hopkins and Christiansen2013) analysis considered discs that were typically assumed to be isothermal, essentially implying that β = 0.
Also, the turbulence they were considering was not the gravito-turbulence driven by the instability itself. The standard analysis of self-gravitating discs considers a scenario in which the quasi-turbulent structures are a self-consistent response to the growth of the instability and depend—as discussed already—on the thermal balance in the disc. Hopkins & Christiansen (Reference Hopkins and Christiansen2013) were considering a situation in which isothermal turbulence is present in a disc in which self-gravity might also operate. In fact, the presence of such turbulence—if the disc were non-isothermal—should act as an extra heating source, which would then be expected to reduce the magnitude of the gravitational instability, rather than enhance it (Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011).
3.3. Triggered fragmentation
Even though the analysis in Hopkins & Christiansen (Reference Hopkins and Christiansen2013) is not quite consistent with the basic picture presented here, it is still interesting from the perspective of whether or not fragmentation could be triggered. This could be via some other form of turbulence (as suggested by Hopkins & Christiansen Reference Hopkins and Christiansen2013) or via some kind of external perturbation. Some early work (e.g., Boffin et al. Reference Boffin, Watkins, Bhattal, Francis and Whitworth1998; Watkins et al. Reference Watkins, Bhattal, Boffin, Francis and Whitworth1998) suggested that stellar encounters could promote fragmentation. Such early simulations typically, however, assumed an isothermal equation of state, which is known to promote fragmentation (Pickett et al. Reference Pickett, Benz, Adams and Arnett1998, Reference Pickett, Cassen, Durisen and Link2000). Later work suggested that the inclusion of compressive and shock heating (Nelson et al. Reference Nelson, Benz and Ruzmaikina2000; Lodato et al. Reference Lodato, Meru, Clarke and Rice2007) would tend to prohibit encounter-driven fragmentation. This was further confirmed using simulations with a realistic radiative transfer formalism (Forgan et al. Reference Forgan, Rice, Stamatellos and Whitworth2009).
The basic argument is that, typically, the encounter will compress the disc on a dynamical timescale, which is much shorter than the thermal timescale. Consequently, the encounter would be expected to heat—and stablise—the disc, rather than produce fragmentation. Fragmentation would only be expected to occur if the cooling time can be reduced by the compression (Whitworth et al. Reference Whitworth, Bate, Nordlund, Reipurth, Zinnecker, Reipurth and Keil2007). This may be possible if the temperature change pushes the gas into a regime where the opacity is reduced (the opacity-gap). However, the corresponding increase in gas density tends to act in the opposite sense, and so this may still make triggered fragmentation unlikely (Johnson & Gammie Reference Johnson and Gammie2003). There has, however, been a recent suggestion (Meru Reference Meru2015) that fragmentation in the inner parts of discs may be triggered by fragments that form, initially, in the outer parts of discs, and Thies et al. (Reference Thies, Kroupa, Goodwin, Stamatellos and Whitworth2010) suggest that tidal encounters might trigger fragmentation in massive, extended protostellar discs. Hence, there may be scenarios under which fragmentation could be triggered.
3.4. Stochasticity
Paardekooper (Reference Paardekooper2012) have suggested that fragmentation may be stochastic. Their basic argument is that a quasi-steady gravito-turbulent state consists of weak shocks and transient clumps that will contract on the cooling timescale (β). If such a clump can avoid being disrupted by one of the spiral shocks, then it could survive to become a bound fragment. For example, Figure 8—taken from Paardekooper (Reference Paardekooper2012)—shows the evolution of the maximum surface density (Σmax/Σ o ) and α in a shearing sheet simulation in which the disc appears to initially settle into a quasi-steady, self-regulated state. This is illustrated by both the maximum density and α settling to reasonably constant values, as indicated by the dashed line in the lower panel. However, at Ωt ~ 300, α starts to decrease and the maximum density increases to a value a few hundred times greater than the initial value, Σ o . This indicates that this system has undergone fragmentation after first appearing to settle to a quasi-steady state, and may be indicative of stochasticity.
Figure—from Paardekooper (Reference Paardekooper2012)—showing the maximum surface density relative to the initial surface density (Σmax/Σ o —top panel) and α (bottom panel) against time, from a shearing sheet simulation that initially appears to have settled into a quasi-steady, self-regulated state. However, at Ωt ~ 300, α starts to decrease and the maximum surface density starts increasing, eventually reaching values many hundreds of times greater than it was initially. The system, which appeared to have settled into a quasi-steady state, has now undergone fragmentation and this may be indicative of stochasticity.
Paardekooper (Reference Paardekooper2012) find that transient clumps are present in all of their simulations (up to β = 50) but that fragmentation only actually happens up to β = 20. Hopkins & Christiansen (Reference Hopkins and Christiansen2013) suggest that, given sufficient time, stochastic fragmentation may actually occur both if the cooling time is very long and the Q value is very high. However, as we suggest above, the scenario they’re considering isn’t quite the same as that being discussed here.
Young & Clarke (Reference Young and Clarke2016) argue that this stochasticity is consistent with the idea behind their suggestion that there are two modes of fragmentation (Young & Clarke Reference Young and Clarke2015). In discs with rapid cooling, fragmentation can occur on the free-fall time. In discs that are cooling slowly (for example, β > 10–20), fragmentation can still occur if a clump can cool and contract prior to encountering a spiral shock. Young & Clarke (Reference Young and Clarke2016) find that the wait time between shocks is not strongly dependent on the cooling time. Given that clumps contract on the cooling timescale, this means that the likelihood of fragmentation decreases with increasing cooling time. They also find an exponential decay in the likelihood of a patch remaining unshocked, and so the probability of a transient clump contracting sufficiently so as to survive a shock encounter is effectively zero for very long cooling times. Given this, even if fragmentation is stochastic, this is unlikely to substantially change the conditions under which fragmentation actually occurs. Similarly, Baehr & Klahr (Reference Baehr and Klahr2015) impose a cooling implementation that attempts to incorporate optical depth effects, and find that this also appears to inhibit stochastic fragmentation.
3.5. Does fragmentation actually occur?
If we focus specifically on discs around young stars (protostellar discs), then we can determine where fragmentation could happen, and even if it does actually happen. It is likely that the inner regions of protostellar discs are optically thick and, consequently, have long cooling times. This is illustrated in Figure 9 which shows α plotted against radius for a quasi-steady disc in which the cooling time at each radius is calculated using the actual optical depth, rather than by simply assuming some β-value a priori. The α values are extremely small in the inner disc but increase with increasing radius, reaching values where fragmentation may be possible in the outer parts of the disc. Consequently, we expect fragmentation to only occur in the outer parts of such systems (Matzner & Levin Reference Matzner and Levin2005; Rafikov Reference Rafikov2005; Stamatellos & Whitworth Reference Stamatellos and Whitworth2008), beyond ~ 30–50 AU.
Figure showing α against radius in a quasi-steady, self-gravitating disc with realistic cooling. The inner disc is optically thick, leading to long cooling times and very small α values. The α value increases with radius, reaching value where fragmentation might be possible in the outer parts of the disc. Hence, we expect fragmentation to only be possible in the outer parts of protostellar accretion discs.
Even if fragmentation is stochastic, or can occur for longer cooling times than we currently think, this has little impact when it comes to fragmentation in protostellar discs. The radial dependence of the cooling time is very strong in the inner parts of these discs (Clarke Reference Clarke2009; Rice & Armitage Reference Rice and Armitage2009). Even a factor of a few in the cooling time for fragmentation has little effect on where fragmentation likely happens (Young & Clarke Reference Young and Clarke2016). If the opacity is significantly lower than interstellar, then the fragmentation radius could to move towards the inner parts of the disc (Meru & Bate Reference Meru and Bate2010). However, it is still unlikely that fragmentation can take place inside ~ 10–20 AU in such discs (Meru & Bate Reference Meru and Bate2010).
In addition to where we might expect fragmentation to occur, we can also estimate typical fragment masses (Boley et al. Reference Boley, Hayfield, Mayer and Durisen2010; Forgan & Rice Reference Forgan and Rice2011). Typically, the initial fragment masses are a few Jupiter masses. Given that we expect fragmentation to only occur in the outer parts of protostellar discs, and to have initial masses of at least a few Jupiter masses, it’s been suggested this could explain some of the directly imaged planets on wide orbits (Kratter et al. Reference Kratter, Murray-Clay and Youdin2010; Nero & Bjorkman Reference Nero and Bjorkman2009). It’s also been suggested that in some cases such objects could rapidly spiral into the inner disc (Baruteau, Meru, & Paardekooper Reference Baruteau, Meru and Paardekooper2011; Michael, Durisen, & Boley Reference Michael, Durisen and Boley2011; Malik et al. Reference Malik, Meru, Mayer and Meyer2015) and could lose mass via tidal stripping (Nayakshin Reference Nayakshin2010; Boley et al. Reference Boley, Hayfield, Mayer and Durisen2010). This could potentially then be the origin of some of the known exoplanets, even terrestrial/rocky exoplanets if the cores can survive disruption (Nayakshin Reference Nayakshin2011). If not, the disruption of an embryo could lead to the formation of planetesimal belts (Nayakshin & Cha Reference Nayakshin and Cha2012) or change the overall disc chemistry (Boley et al. Reference Boley, Hayfield, Mayer and Durisen2010). In some cases, however, it may be possible for the protoplanet to open a gap and survive (Stamatellos Reference Stamatellos2015).
Forgan & Rice (Reference Forgan and Rice2013), however, argue that typically fragments are either destroyed, or remain at large radii with relatively high masses (many Jupiter masses). According to their analysis, the formation of planetesimal belts or terrestrial planets via this mechanism is also rare. The high destruction rate, however, is consistent with other simulations (Boley & Durisen Reference Boley and Durisen2010; Zhu et al. Reference Zhu, Hartmann, Nelson and Gammie2012) and could lead to outbursts, such as the FU Orionis phenomenon (Vorobyov & Basu Reference Vorobyov and Basu2005; Dunham & Vorobyov Reference Dunham and Vorobyov2012). Galvagni & Mayer (Reference Galvagni and Mayer2014), however, suggest that a large fraction of the clumps could survive inward migration and that such a process could explain a reasonable fraction of the known population of ‘hot’ Jupiters and other gas giant planets.
In addition to the rapid inspiral of fragments, it’s also possible that objects forming via fragmentation at large radii could be scattered by stellar, or other, companions (Forgan, Parker, & Rice Reference Forgan, Parker and Rice2015). Some of these could then be circularised onto a short-period orbit, via tidal interactions with the parent star (Lin, Bodenheimer, & Richardson Reference Lin, Bodenheimer and Richardson1996), producing either a ‘hot’ Jupiter or a proto-hot Jupiter (Dawson & Johnson Reference Dawson and Johnson2012). Rice et al. (Reference Rice, Lopez, Forgan and Biller2015), however, argue that the population of known ‘hot’ and proto-hot Jupiters is inconsistent with this being common and suggests that disc fragmentation rarely forms planetary-mass objects. This is also consistent with the relatively low frequency of known wide orbit planetary and brown dwarf companions (Biller et al. Reference Biller2013; Brandt et al. Reference Brandt2014).
4 DISCUSSION AND CONCLUSION
We’ve tried here to summarise our current understanding of the evolution of self-gravitating accretion discs, focusing—in particular—on those around young’s protostars. When self-gravity is important, we expect such discs to either settle into a quasi-steady, self-regulated state in which heating and cooling are in balance and Q ~ 1 (Paczynski Reference Paczynski1978), or—if the instability is particularly strong—to fragment to form bound objects.
Typically, the boundary between a disc fragmenting, or settling into a quasi-steady state, is determined by the rate at which it is losing energy, with rapid cooling potentially leading to fragmentation (Gammie Reference Gammie2001). The exact boundary does, however, depend on whether or not there is an additional heating source (Rice et al. Reference Rice, Armitage, Mamatsashvili, Lodato and Clarke2011), and on the mass of the disc relative to the mass of the central object; due to the global nature of the gravitational instability (Balbus & Papaloizou Reference Balbus and Papaloizoum1999), very massive discs may avoid fragmentation under conditions that may lead to fragmentation in less massive discs (Lodato & Rice Reference Lodato and Rice2005; Forgan et al. Reference Forgan, Rice, Cossins and Lodato2011). Similarly, very massive discs may not settle into a true quasi-steady, self-regulated state, but instead show variability, and are only quasi-steady in a time averaged sense. It is also possible that fragmentation itself is stochastic and could, at times, occur for longer cooling times than this basic analysis suggests (Paardekooper Reference Paardekooper2012).
In protostellar discs, however, the exact fragmentation boundary is not that relevant when it comes to determining where fragmentation might happen. The inner parts of such discs are so optically thick, and cool so slowly, that fragmentation is only actually likely in the outer parts of these discs (Rafikov Reference Rafikov2005). Even if fragmentation could occur for slightly longer cooling times than originally thought, or if fragmentation is stochastic, this is unlikely to significantly influence where fragmentation can occur (Young & Clarke Reference Young and Clarke2016). Although, reductions in opacity can bring the fragmentation region closer to the central star (Meru & Bate Reference Meru and Bate2010; Rogers & Wadsley Reference Rogers and Wadsley2012), it’s still unlikely close to the central star. Essentially, it is very difficult for there to be conditions such that fragmentation could occur inside ~ 20 AU.
There are, however, still aspects that are uncertain. It might seem that the lack of convergence in self-gravitating disc simulations (Meru & Bate Reference Meru and Bate2011) is most likely numerical (Lodato & Clarke Reference Lodato and Clarke2011; Meru & Bate Reference Meru and Bate2012; Rice et al. Reference Rice, Forgan and Armitage2012, Reference Rice, Paardekooper, Forgan and Armitage2014), but this is still not fully resolved. Similarly, there is evidence for stochasticity (Paardekooper Reference Paardekooper2012) and a potential explanation for this (Young & Clarke Reference Young and Clarke2016), but there are still aspects of this that are uncertain. However, our basic understanding seems robust. Self-gravitating discs will tend to settle into a quasi-steady, self-regulated state in which heating balances cooling, in which the transport of angular momentum can be described as pseudo-viscous (Lodato & Rice Reference Lodato and Rice2004), and in which the amplitude of the perturbations depends on the magnitude of the effective viscosity α (Cossins et al. Reference Cossins, Lodato and Clarke2009). If the perturbations become very large, then the disc may undergo fragmentation.
It’s also unclear as to whether or not fragmentation ever actually occurs in protostellar discs. There are suggestions (Kratter et al. Reference Kratter, Murray-Clay and Youdin2010; Nero & Bjorkman Reference Nero and Bjorkman2009) that it may explain some of the directly imaged exoplanets on wide-orbits (Marois et al. Reference Marois, Macintosh, Barman, Zuckerman, Song, Patience, Lafreniére and Doyon2008) and suggestions that some may spiral inwards to smaller orbital radii than where they form (Nayakshin Reference Nayakshin2010). On the other hand, the paucity of directly imaged planets on wide-orbits (Biller et al. Reference Biller2013) and the properties of the known closer-in exoplanets, suggests that disc fragmentation may rarely form planetary-mass objects (Rice et al. Reference Rice, Lopez, Forgan and Biller2015).
However, disc self-gravity may still play a crucial role in angular momentum transport during the earliest stages of star formation, and the resulting density waves may then play a role in the growth of planet building material (Rice et al. Reference Rice, Lodato, Pringle, Armitage and Bonnell2004; Gibbons, Rice, & Mamatsashvili Reference Gibbons, Rice and Mamatsashvili2012). So, even if the gravitational instability does not play a significant role in forming planets directly, it may still play an important role in the growth of the building blocks crucial for planet formation.
ACKNOWLEDGEMENTS
The author would like to thank the organisers of the Disc Dynamics and Planet Formation meeting, held in Cyprus, that partly motivated this paper. The author would also like to thank the reviewer for a very thorough and detailed review that helped to improve the paper, and would like to acknowledge useful comments from Phil Armitage. KR gratefully acknowledges support from STFC grant ST/M001229/1. Some of the research leading to these results also received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement number 313014 (ETAEARTH).
