Quantification of stochastic fragmentation of self-gravitating discs

Abstract

1 INTRODUCTION

It is widely believed that planets can form via fragmentation of disc gas only in the outermost regions of protoplanetary discs (i.e. >50 au; Rafikov 2005; Clarke 2009; Rice & Armitage 2009). This mechanism is thus favoured for explaining planets in wide orbits as detected via direct imaging (Marois et al. 2008, 2010) for which the time to form by the alternative core accretion model is uncomfortably long. Core accretion (i.e. the assembly of a rock core and subsequent accretion of disc gas) is therefore favoured for the formation of the vast majority of gas giant planets detected to date.

The argument against forming planets by fragmentation in the inner disc hinges on the fact that self-gravitating discs are very optically thick at smaller radii which renders β (the ratio of cooling time to dynamical time) very large. A large number of numerical experiments have shown that whether a self-gravitating disc fragments or whether it is maintained in a self-regulated (non-fragmenting) state depends on the value of β. We will not here repeat the considerable debate about the exact value of the critical β value for fragmentation, nor its dependence on numerical technique and resolution (Gammie 2001; Rice et al. 2003, 2014; Rice, Lodato & Armitage 2005; Lodato & Clarke 2011; Meru & Bate 2011, 2012; Rice, Forgan & Armitage 2012; Young & Clarke 2015). Suffice it to say that the range of values of β_crit claimed in the literature (3–30), corresponds to a rather narrow range of minimum radii for fragmentation because the radial dependence of β in the outer regions of a steady state self-gravitating disc is steep (β ∝ r^−4.5) (Rafikov 2005, 2009; Clarke 2009).

All of the above numerical experiments refer to prompt fragmentation (i.e. fragmentation that occurs within a cooling time of the disc being set up in a self-gravitating state). Paardekooper (2012) however, reported a distinct fragmentation mode in long time-scale integrations using the grid code fargo 2d. In this experiment, the disc apparently settled into the self-regulated state and then eventually formed a fragment after a large number of cooling times. He dubbed this alternative mode of fragmentation ‘stochastic fragmentation’.

The prospect of eventual fragmentation even at high β values could potentially overturn the conclusions drawn above about the impossibility of forming planets by gravitational fragmentation at small radii. Given that the self-gravitating lifetime of a disc is >10⁴ dynamical times at the radius of Jupiter, it follows that even an extremely small (but non-zero) probability of fragmentation per dynamical time could lead to planet formation by this route.

The possibility of stochastic fragmentation in self-gravitating media has also been discussed by Hopkins & Christiansen (2013). In their picture, stochastic fragmentation may result from the collapse of extremely rare non-linear density fluctuations on a scale much smaller than the Jeans length; their study quantifies the probability of this outcome using the statistics of isothermal turbulence. Although those studies that have attempted to quantify the power spectrum of the turbulence have found little power at sub-Jeans scales (Boley et al. 2007; Cossins, Lodato & Clarke 2009), we cannot rule out the possibility that such a mechanism is operative in self-gravitating discs (indeed no simulations performed to date have the extremely sub-Jeans length resolution required for this analysis). However, this is definitely not the mechanism that is operating in the simulations of Paardekooper (2012). Here, collapse is initiated on the Jeans scale (as expected) and the protoclump contracts on the (long) cooling time of the simulation. With such a very slow contraction, the expected outcome is that protoclumps are disrupted by collisions with the spiral arms. The ‘stochastic’ (and unexpected) outcome of the Paardekooper simulation was that clumps can occasionally avoid such collisions over a sufficiently long time window that they manage to collapse.

We have argued above that it is essential to quantify the probability of stochastic fragmentation before ruling out fragmentation of the inner disc. Such quantification is however extremely difficult through direct simulation. By definition, stochastic fragmentation simulations must employ large β values in order to avoid the well-studied prompt fragmentation regime. Spiral structure is weak at large β (Cossins et al. 2009) and so high resolution is required in order that weak spiral heating is not swamped by artificial viscosity (Lodato & Clarke 2011). The nature of stochastic fragmentation also implies that long integration times (many cooling times) are required. The combination of high resolution and long integration times implies that direct simulation is not likely to be a feasible strategy in the foreseeable future.

Here we adopt a different approach. We argued in Young & Clarke (2015) that the β threshold for prompt fragmentation can be simply understood in terms of the mean spiral structure in self-gravitating discs. Cossins et al. (2009) presented a detailed characterization (wavenumbers, pattern speed, pitch angles) of the spiral modes present in such discs which allows the calculation of the mean time between spiral arm passages. Young & Clarke (2015) interpreted the fragmentation boundary as corresponding to the criterion that a protofragment can cool on a time-scale less than this mean interspiral time; if this condition is not fulfilled then protofragments are disrupted by spiral arm passage before they collapse to the small scales where they are immune to such disruption.

If this was the full story then there would be no possibility of stochastic fragmentation at larger β. However a casual inspection of simulations of self-gravitating discs is enough to show that individual spiral features form, dissolve and re-form on a time-scale that is close to dynamical. Although the Cossins et al. (2009) description is a reasonable representation of the time-averaged behaviour, it does not address the possible variability in the interspiral time. In the fluctuating density field encountered in the simulations it may be possible for individual protofragments to avoid spiral shocks over significantly longer periods.

In this paper, we therefore analyse the intermittent nature of spiral structure in self-gravitating discs, recording the distribution of interspiral times for fluid elements in the disc. We expect from Cossins et al. (2009) that the geometry (though not the amplitude) of spiral features is not a strong function of either β or resolution but test this hypothesis with a range of simulations. We emphasize that we are not seeking to model fragmentation itself (and choose simulation parameters with sufficiently high β that prompt fragmentation is avoided) but instead use the Lagrangian nature of smoothed particle hydrodynamics (SPH) to analyse particle histories and infer the distribution of wait times between spiral arm encounters. In going on to draw conclusions about the survival prospects of any protofragments that might eventually form in the disc, we will be assuming that such fragments would be sampling the same fluctuating density field that we analyse here.

2 DISC MODEL

All our simulations were run in 2D, with ξ = 5 and q = 0.2 to minimize computational expense (Young & Clarke 2015). Each simulation was run for at least 10 cooling times at the outer edge, ensuring that the disc reached a ‘settled’ state. We limited our analysis to the final 500 t_dyn of each simulation (200 t_dyn for the highest resolution run), where the disc was in the Q ∼ 1 gravoturbulent state.

3 RESULTS

Paardekooper (2012) and Young & Clarke (2015) interpreted stochastic fragmentation as arising from unstable overdensities that survive for long enough to become bound and resist disruption. That is, the gravoturbulent, Q ∼ 1 state is constantly producing gravitationally unstable overdensities. The majority of these overdensities are disrupted before they can become bound enough to survive disruption by the environment. If an overdensity is ‘lucky’ and can survive for long enough without being disrupted, it will form a fragment.

To quantify the likelihood of stochastic fragmentation, we measured the distribution of times between successive spiral shock wave encounters for many particles within our quasi-stable simulations. Because reliably detecting shock fronts is numerically challenging (Morris & Monaghan 1997; Cullen & Dehnen 2010), we only consider the relatively strong shocks in the disc. We tracked individual particles over roughly 45, 90, 180 or 360 dynamical times and marked all periods of increase in the particle's entropy and surface density (some of which spanned multiple snapshots) as potential encounters with shocks. We deemed a shock to have been encountered whenever a particle's entropy and surface density increased by more than the median increase across all potential shocks.¹ We limited our analysis to particles within R = 2–3 R_i to avoid corruption by edge effects.

Fig. 1 shows the resulting distributions of wait time between shocks for simulations at a range of resolutions, values of β and particle tracking times. There is no significant difference in the average wait times between simulations at different resolutions or with different cooling rates. Furthermore, the distributions are converged with respect to the length of time a particle is tracked for, indicating that neither increased resolution nor longer simulation runs would change the result.

Figure 1.

Distributions of wait times between shocks for simulations at different resolutions and β, inferred by tracking particles for different lengths of time. Particles were tracked across 128 (solid lines), 256 (dashed lines), 512 (dot–dashed lines) and 1024 (dotted lines) dynamical times at the disc's inner edge and all increases of entropy were recorded. A shock was deemed to have been encountered whenever the entropy and surface density increased by more than the median increase amongst all potential shocks (see text for details).

Open in new tabDownload slide

Most importantly, all our simulations show an exponential decay in the likelihood of a patch of disc remaining unshocked, with long periods of time (100 s of dynamical times or more) having effectively zero probability in all cases. This is shown more clearly in Fig. 2, which shows the log₁₀ cumulative probability of surviving for longer than a certain number of dynamical times without encountering a shock. The slope change at the extreme tail of each distribution is caused by low number statistics resulting from our finite sample of particles and finite tracking time. However, it is clear that the probability of a clump surviving for longer than ∼10 t_dyn without encountering a strong shock decreases exponentially.

Figure 2.

log₁₀ of the probability that the time between shocks will be at least the number of dynamical times shown on the x-axis. Particles were tracked across 128 (solid lines), 256 (dashed lines), 512 (dot–dashed lines) and 1024 (dotted lines) dynamical times at disc's the inner edge. See text and Fig. 1 for further details of how the time between shocks was calculated.

Open in new tabDownload slide

Fig. 3 shows the probability of fragmentation as a function of radius for a disc around a solar mass star with a self-gravitating lifetime of 100 000 yr, where the blue line adopts the assumption discussed thus far, i.e. that fragmentation occurs only for clumps that avoid a spiral arm encounter for t_cool. The green line in Fig. 3 is a less stringent requirement and instead involves the assumption that the clump has been able to contract to a density significantly larger than the mean post-shock density before it meets a spiral arm. At large radii, where β is low and the spiral features are of high amplitude this is equivalent to the criterion described above. In the inner disc, where the cooling time is long and the spiral features are of low amplitude it implies a more permissive criterion for collapse; in this case we use the arm amplitude measured in our simulations combined with the β dependence of arm amplitude derived by Cossins et al. (2009) in order to require that the clump's density has increased by a factor of |$20\sqrt{20/\beta }$| between spiral arm encounters before it is deemed to be able to collapse.

Figure 3.

The probability of stochastic fragmentation occurring during the self-gravitating lifetime of a disc at different radii calculated from equation (7) for a disc around a solar mass star with a self-gravitating phase lasting 100 000 yr. The solid and dashed lines show the probability of fragmentation assuming, respectively, that fragmentation occurs when a clump survives for more than t_cool and that fragmentation occurs when a clump exceeds a certain critical density.

Open in new tabDownload slide

It is clear from Fig. 3 that the transition between fragmenting and non-fragmenting takes place over a very narrow range of radii. This is true whether we require a clump to meet a density cut-off or to survive for a cooling time before we consider it to have fragmented: indeed the two criteria are visually nearly indistinguishable in Fig. 3 because by the time one enters the high β regime at small radius where the two criteria diverge, the probabilities are vanishingly small. We thus conclude that stochastic fragmentation cannot decrease the maximum radius at which fragmentation is possible by more than a few tens of per cent.

4 DISCUSSION

In this study, we have chosen to only perform 2D simulations where higher resolution in h/H can be obtained more easily. It has recently been shown that 2D simulations are problematic for studying fragmentation, due to the effect gravitational softening has on preventing pressure-supported collapse (Young & Clarke 2015). However, it is important to note that the suppression of fragmentation in 2D results entirely from the modification to the gravitational force on small scales. Our conclusions rely entirely on results derived from the large-scale spiral structure of the disc, which is unaffected by any small-scale differences in the gravitational force law between 2D and 3D. As such, we expect the distribution of times between spiral wave encounters to be unchanged in 3D.

Figs 1 and 2 show that the spiral morphology of the disc, and hence the time between successive shock fronts, does not vary significantly with β or resolution. This is consistent with the findings of Cossins et al. (2009), who investigated the spiral geometry of 3D gravoturbulent accretion discs and found no β or resolution dependence. To the extent to which there is any trend with β and resolution, it is towards times between shock fronts becoming shorter at higher resolution and/or β, which would only strength our conclusions (i.e. make stochastic fragmentation at small radii more difficult).

Naturally, the cooling rates in a realistic disc are not likely to be well described by cooling at constant β once the fragments enter the regime of non-linear collapse. This consideration is unlikely to have a major bearing on our analysis since what is relevant is the cooling regime, while the fragments’ overdensity with respect to the surroundings is still relatively modest and where the cooling rates are therefore not expected to deviate greatly from those appropriate to the background disc. Following Cossins et al. (2010), who examined the effect of temperature dependent cooling on prompt fragmentation, we would not expect this to facilitate fragmentation apart from close to strongly temperature-dependent features in the opacity law, e.g. due to ice sublimation (see also Johnson & Gammie 2003). Although this enhanced cooling of overdense regions may bring the fragmentation boundary inwards by a modest factor, it is unlikely to change our conclusions regarding the sharpness of this transition.

In focusing on the time between spiral wave passages, we have assumed that contracting clumps are disrupted by shocks. In apparent conflict with this assumption, it has been shown that when fragmentation is immediate, fragments arise out of the dense post-shock gas that makes up the spiral arms (Rogers & Wadsley 2012). However, in this scenario the post-shock gas is the seed of the instability (which then rapidly contracts to form a fragment), but does not promote the contraction of an existing clump. It is possible that the dense post-shock gas is the location where the gravitational instability is most commonly triggered. None the less, for stochastic fragmentation to occur, this unstable patch of disc must contract significantly from its ‘birth density’. Paardekooper (2012) found no evidence that the stochastic fragments formed in his simulations were driven to fragment by spiral wave encounters. Additionally, if it was the case that spiral waves promoted clump contraction, it is difficult to see how any gravitationally unstable patch of disc could fail to progress into a fragment.

Fig. 1 shows that the most common wait time between shocks is ∼4–5 t_dyn. This suggests that mξ ≈ 1 in equation (5), in good agreement with the findings of Cossins et al. (2009). Note that the ‘fragmentation boundary’ measured by Meru & Bate (2011, 2012) and Rice et al. (2012, 2014) corresponds to the point where the probability of surviving for more than β t_dyn becomes low, not to the most common wait time. Given this, our results suggest the fragmentation boundary should be in the range β = 7–12, although the uncertainty in the number of overdense regions formed per dynamical time (see below) prevents us from placing stronger constraints.

Paardekooper (2012) found that clumps formed in the gravoturbulent state for all values of β up to β ∼ 50, but was only able to find stochastic fragmentation up to β = 20. The probability of a simulation fragmenting in its lifetime is approximately the probability shown in Fig. 2 (which is the probability of fragmentation per clump) multiplied by the number of overdense clumps that form in the simulations lifetime. Paardekooper (2012) found roughly one clump per ∼50 t_dyn in his β = 9 simulations and found two simulations fragmented after ∼600 and ∼800 t_dyn, while two simulations remained stable for 1000 t_dyn. This implies a probability of fragmentation per clump of ∼2/68 ≈ 3 per cent (3400/50 = 68 clumps in 3400 t_dyn), which is broadly consistent with the ∼10 per cent predicted by Fig. 2. By β = 20, Fig. 2 predicts that only ∼0.5 per cent of clumps will fragment, suggesting fragmentation is unlikely but not impossible.

We therefore conclude that stochastic fragmentation only modifies the ‘fragmentation boundary’ in radius by ∼20 per cent (see Fig. 3) and rules out the formation of giant planets via any type of gravitational collapse (stochastic or otherwise), except in the outermost parts of protoplanetary discs. This agrees with the conclusions of earlier disc fragmentation studies (Rice et al. 2005, 2014; Paardekooper, Baruteau & Meru 2011; Meru & Bate 2012). The only remaining avenue for fragmentation at low radii is via extremely large sub-Jean's length density fluctuations (Hopkins & Christiansen 2013), which cannot presently be directly probed computationally.

5 CONCLUSION

In this paper, we have used simple, 2D SPH simulations of self-gravitating discs to place hard, quantitative limits on the values of the cooling time parameter, β, for which massive discs can stochastically fragment. We find that the average time any patch of disc spends between spiral wave encounters is broadly consistent with the findings of Young & Clarke (2015) and that the stochastic nature of the disc's spiral structure can only increase this time by at most a factor of a few.

Consistent with earlier studies (Cossins et al. 2009), we find that the distribution of wait times between strong shocks is independent of both resolution and β. Furthermore, the distribution of wait times decreases exponentially for wait times beyond ∼15 t_dyn. Combining this finding with the expected β ∝ R^−9/2 scaling of cooling time with radius, we find that stochastic fragmentation is unable to modify the radius at which fragmentation is possible by more than a few tens of per cent. This finding restricts direct collapse as a mechanism for giant planet formation to the outer parts of protoplanetary discs.

6 MATERIALS AND METHODS

In the interests of reproducibility and transparency, all code and data used in performing this work have been made freely available online at https://bitbucket.org/constantAmateur/discfragmentation.

This work benefited greatly from discussions with Giuseppe Lodato, Deborah Sijacki, Richard Nelson, Farzana Meru, Richard Booth and Daniel Price on issues of both physics and numerics. Sijme-Jan Paardekooper deserves special thanks for having provided advice and feedback throughout the course of this project.

MY gratefully acknowledges the support of a Poynton Cambridge Australia Scholarship. This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG.

This work used the DIRAC Shared Memory Processing system at the University of Cambridge, operated by the COSMOS Project at the Department of Applied Mathematics and Theoretical Physics on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/J005673/1, STFC capital grant ST/H008586/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure.

REFERENCES