On the fragmentation criteria of self-gravitating protoplanetary discs

Abstract

1 INTRODUCTION

The evolution of protoplanetary discs has been explored at great length in the past to understand the processes by which they form in the early Class 0 phase, the accretion from the molecular cloud core on to the disc in the Class I phase, the mass and angular momentum transfer in discs once they have formed as well as in the early Class II stage of an isolated disc, as well as the disc dispersal mechanisms. One such concept that has been considered is the importance of the disc self-gravity, particularly in the earlier period of a disc's lifetime. It is during the early stages when it is massive enough to be self-gravitating that the concepts of angular momentum transport and fragmentation in these discs are important. This is a particularly significant aspect when considering whether gas giant planets could form via this method. Historically, planet formation by gravitational instability has not been thought likely since the planets that form in this way are predicted to do so far out in a disc [≳O(100) au; e.g. Matzner & Levin 2005; Rafikov 2005; Stamatellos & Whitworth 2008; Boley 2009; Clarke 2009; Rafikov 2009; Kratter, Murray-Clay & Youdin 2010), whereas until recently, giant planets have been found only at small radii (≲10 au).

Recent advances in observations, such as the discoveries of planets at large radii [O(10–100) au] from the parent star (Kalas et al. 2008; Marois et al. 2008; Lafrenière, Jayawardhana & van Kerkwijk 2010), call for an alternative planet formation mechanism other than the standard core accretion method to be considered. Furthermore, theoretical advances, such as the increased likelihood that planets can form by gravitational instability closer to the parent star in low metallicity environments (Meru & Bate 2010) where core accretion finds it difficult (Kornet et al. 2005), call for planet formation by gravitational instability to be scrutinised in much more detail.

There are two quantities that have historically been used to determine whether a disc is likely to fragment. The first is the stability parameter (Toomre 1964),

1

where c_s is the sound speed in the disc, κ_ep is the epicyclic frequency, which for Keplerian discs is approximately equal to the angular frequency Ω, Σ is the surface mass density and G is the gravitational constant. Therefore, once the surface mass density and the rotation of the disc have been established, the stability is purely dependent on the disc temperature. Toomre (1964) showed that for an infinitesimally thin disc to fragment, the stability parameter must be less than a critical value, Q_crit≈ 1.

Gammie (2001) showed that in addition to the above stability criterion, the disc must cool at a fast enough rate. Using shearing sheet simulations, he showed that if the cooling time-scale can be parametrized as

2

where

3

u is the specific internal energy and du_cool/dt is the total cooling rate, then for fragmentation we require β≲ 3, for a ratio of specific heats γ = 2 (in two dimensions). Rice, Lodato & Armitage (2005) carried out three-dimensional simulations using a smoothed particle hydrodynamics (SPH) code and showed that this cooling parameter is dependent on the equation of state. They showed that fragmentation can occur if β < β_crit where β_crit≈ 6 –7 for discs with γ = 5/3 and β_crit∼ 12 –13 for discs with γ = 7/5. Gammie (2001) and Rice et al. (2005) also showed that in a steady-state disc where the dominant form of heating is that due to gravitational instabilities, since the gravitational stress in a disc can be linked to the cooling time-scale in the disc using

4

the rapid cooling required for fragmentation can also be interpreted as a maximum gravitational stress that a disc can support without fragmenting, which they show to be α_GI,max≈ 0.06.

The concept of a fast cooling needed for fragmentation is very clear from previous work. However, the value of the critical cooling time-scale, β_crit (and therefore by equation 4, α_GI,max), does not appear to be too clear-cut: Rice et al. (2003) found that for a 0.1 -M_⊙ disc with a surface mass density profile, Σ∝R^−7/4, extending to a radius, R_out = 25 au, around a 1 -M_⊙ star, the disc fragments using β = 3 but not for β = 5, whereas for a disc with mass M_disc = 0.25 M_⊙, the disc fragments for β = 5. On the other hand, Rice et al. (2005) suggest that the fragmentation boundary is independent of the disc to star mass ratio. Clarke, Harper-Clark & Lodato (2007) showed that the critical value of β (below which fragmentation will occur if the stability criterion is met) may depend on the disc's thermal history: if the time-scale on which the disc's cooling time-scale is decreased is slower than the cooling time-scale itself (i.e. a gradual decrease in β), then the critical value may decrease by up to a factor of 2. More recently, Cossins, Lodato & Clarke (2010) showed that the critical value varies with the temperature dependence of the cooling law. In addition, they carry out a simulation of a self-gravitating disc with a surface mass density profile, Σ∝R^−3/2[cf. Rice et al. (2005) who used Σ∝R⁻¹], with a ratio of specific heats, γ = 5/3, and show that the critical value β_crit≈ 4. Using equation (4), this is equivalent to α_GI,max = 0.1 which brings the result of α_GI,max = 0.06 described above into question. Yet a number of papers have been produced which base their work on the concept of a single critical value of β (or, equivalently, a maximum gravitational stress value; e.g. Clarke 2009; Rafikov 2009; Cossins et al. 2010; Kratter et al. 2010).

However, there is a more fundamental question that arises from previous numerical studies using a parametrized cooling method. In many simulations (e.g. Rice et al. 2003, 2005; Clarke et al. 2007), the fragments all form in the outer parts of the discs. If the fragmentation criterion for a self-gravitating disc only depends on β, then if fragments form, one would expect them to form at all radii since all radii have the same value of β. This implies that the cooling time-scale is not the only parameter on which fragmentation depends.

In this paper, we investigate the criteria for fragmentation in greater detail. In Section 2, we analytically investigate how fragmentation may be expected to depend on various disc parameters. We test this analytical theory by carrying out global three-dimensional simulations, the numerical set-up of which we describe in Section 3. In Section 4, we make initial comparisons between our simulations and previous simulations by Rice et al. (2005) as well as discuss the implications of the disc set-up. In Section 5, we present our simulations, the results of which we describe in Section 6. We finally discuss and make conclusions in Sections 7 and 8, respectively.

2 ANALYTICAL VIEW

As discussed in Section 1, for a disc to fragment, one criterion is that the Toomre stability parameter Q≲ 1. Making the approximation that and using H=c_s/Ω, where H is the isothermal scaleheight of the disc, the Toomre stability criterion becomes a condition on the aspect ratio, H/R:

5

Approximating the surface mass density as Σ≈M_disc/(πR²), equation (5) becomes (Gammie 2001)

6

The surface mass density of a disc can also be written in the form of a power law, Σ=Σ₀(R₀/R)^p, where Σ₀ is the surface mass density at radius R₀ and p is a constant for any one disc. Substituting this into equation (5), the condition for fragmentation becomes

7

Equations (6) and (7) show the following.

• Increasing the disc mass or decreasing the star mass is likely to promote fragmentation since a greater portion of the disc is likely to fulfil the above criteria.

• The surface mass density profile may play a part in the fragmentation of a disc. If p = 2, the disc is scale-free and each radius is equivalent to any other radius: the right-hand side (RHS) of equation (7) is constant with radius and the ratio of the cooling time-scale to the orbital time-scale, β, is also a constant. Therefore, if the disc settles into a quasi-steady state where the internal heating due to the gravitational instabilities balances the cooling, we expect Q also to be constant with radius [i.e. the left-hand side (LHS) of equation 7, H/R, is also a constant]. Consequently, if p = 2, either the entire disc should settle into a quasi-steady state or the entire disc should fragment. We note that in a fragmenting disc, because the heating, cooling and fragmentation time-scales are all proportional to the dynamical time-scale in the disc, the fragmentation should first occur (in absolute terms) at small radii.

  For p < 2, the RHS of equation (7) increases with radius. Although H/R is likely to increase with radius as well, since H/R will typically increase more slowly than the RHS of equation (7), this condition is more likely to be satisfied in the outer regions of a disc. Conversely, for p > 2, the RHS of equation (7) decreases with radius and hence the condition is more likely to be satisfied at small disc radii.

• For a p < 2 disc with a low enough β such that it can fragment, the exact value of β may determine how much of the disc satisfies condition (7). If the cooling in a disc is fast such that β is small, the temperature and hence sound speed, c_s, will decrease more rapidly than in a disc where β is higher. Consequently, since H∝c_s, the aspect ratio will be lower at any particular radius and hence the disc is more likely to satisfy condition (7) for smaller β. Since gravitational instability typically develops on a dynamical time-scale, t_dyn∝ 1/Ω∝R^3/2, the instability will develop in the inner regions first, and therefore fragmentation will first occur as close to the inner regions as possible where the fragmentation criteria are satisfied. Since more of the disc satisfies the above criteria for decreasing β, the radius at which the first fragment forms will also decrease. We therefore expect the fragmentation radius to move inwards with more efficient cooling.

• Crucially, equation (7) shows that the radius is important and that for a shallow surface mass density profile (p < 2), there does not appear to be a limit for fragmentation (if the disc cools fast enough): provided the disc is large enough, condition (7) will be satisfied (since typically, an increase in H/R with radius will be much smaller than the increase in the RHS of equation 7).

3 NUMERICAL SET-UP

Our simulations are carried out using an SPH code originally developed by Benz (1990) and further developed by Bate, Bonnell & Price (1995) and Price & Bate (2007). It is a Lagrangian hydrodynamics code, ideal for simulations that require a large range of densities to be followed, such as fragmentation scenarios.

We include the heating effects in the disc due to work done on the gas and artificial viscosity to capture shocks. The cooling in the disc is taken into account using the cooling parameter, β (equation 2), which cools the gas on a time-scale given by equation (3).

In order to model shocks, SPH requires artificial viscosity. We use a common form of artificial viscosity by Monaghan & Gingold (1983), which uses the parameters α_SPH and β_SPH. A corollary of including artificial viscosity is that it adds shear viscosity and causes dissipation. If this viscosity is too large, the evolution of the disc may be driven artificially, while if it is too small, it will lead to inaccurate modelling of shocks (Bate 1995). As discussed in Meru & Bate (2010), various values of the SPH parameters have been tested and we find that a value of α_SPH∼ 0.1 provides a good compromise between these factors. Since typically, β_SPH≈ 2α_SPH, we choose α_SPH and β_SPH to be 0.1 and 0.2, respectively, which are fixed throughout the simulations. Furthermore, our work will begin with a comparison with Rice et al. (2005) and so we use the same values used by them. We use an adiabatic equation of state with a ratio of specific heats, γ = 5/3.

3.1 Numerical effects on fragmentation results

Rice et al. (2005) showed that for a disc with a ratio of specific heats, γ = 5/3, the critical value of the cooling time-scale in units of the orbital time-scale required for fragmentation is β_crit∼ 6 –7. The SPH code used for the simulations presented in this paper differs in the way the smoothing length of the particles is set from that of Rice et al. (2005): whilst their code sets the smoothing length by approximately fixing the number of neighbours that each particle has to ≈50, the current version uses a variable smoothing length which does not fix the number of neighbours but allows the smoothing length to be spatially adaptive whilst maintaining energy and entropy conservation (Monaghan 2002; Springel & Hernquist 2002), with our particular implementation described by Price & Bate (2007).

All explicit hydrodynamical simulations must limit the time-step based on the Courant condition. Our SPH code also applies a force condition and a viscous time-step condition (see Monaghan 1992 for a review). In the simulations presented here, since we apply an explicit cooling rate, it is important to ensure that the following time-step criterion is also met:

8

where C is a constant less than unity. We investigate the effects of varying the constant C on the critical cooling time-scale β_crit by testing values of C = 0.3, 0.03 and 0.003. We find that this does not have a significant effect on the fragmentation results and so we use C = 0.3 for the simulations presented here. However, the time-step criterion may become more important for small β or for particles at small radii. Therefore, for those simulations carried out with small values of β(≤3) or where fragmentation occurs at small radii (≲5 au), the simulations have been repeated with C = 0.03 to confirm that this does not play a part in the results.

4 BENCHMARKING SIMULATIONS

Table 1 summarizes the parameters and fragmentation results of the simulations presented here. Each simulation was run either beyond the point at which the disc attained a steady state [for greater than six outer rotation periods (ORPs)] or until it fragmented.

The simulations presented by Rice et al. (2005) also used an SPH code. However, since the way the smoothing length is set in our code differs from the code used by Rice et al. (2005), and since it is uncertain as to whether their time-stepping considered the cooling time-scale, we simulate the same disc that Rice et al. (2005) simulated in order to initially find the critical cooling time-scale in units of the orbital time-scale, β_crit. This is done by setting up a 1 -M_⊙ star with a 0.1 -M_⊙, disc made of 250 000 SPH particles, spanning 0.25 ≤R≤ 25 au. The initial surface mass density and temperature profiles of the disc are Σ∝R⁻¹ and T∝R^−1/2, respectively. The magnitudes of these are set such that the Toomre stability parameter (equation 1) at the outer edge of the disc Q_min = 2. This gives an aspect ratio H/R≈ 0.05. We model the 1 -M_⊙ star in the centre of the disc using a sink particle (Bate et al. 1995). At the inner disc boundary, particles are accreted on to the star if they move within a radius of 0.025 au of the star or if they move into 0.025 ≤R < 0.25 au and are gravitationally bound to the star. At the outer edge, the disc is free to expand.

The simulation was run using a ratio of specific heats, γ = 5/3, and hence according to Rice et al. (2005), 6–7. We find that the critical value is ≈5.6 since this is the lowest value of β that the discs can have without fragmenting (compare simulations Benchmark1–3). According to equation (4), this is equivalent to a critical value of the gravitational stress, α_GI,max∼ 0.07, which is similar to the value of ∼0.06 obtained by Rice et al. (2005). Given the differences between the codes, we consider this level of agreement acceptable. We therefore compare our remaining simulations to this value of β_crit.

Fig. 1 shows the initial Toomre stability profile of the Rice et al. (2005) disc (solid line) that is replicated here. As a simulation is started, the disc heats up due to the heating from gravitational instability, the resulting compression and viscous heating, and cools on the cooling time-scale defined by the cooling parameter, β. Consequently, it is expected that the initial disc temperature profile would not play a part in the resulting evolution of the disc. We therefore test this by setting up a disc with the same surface mass density profile, Σ∝R⁻¹, but with a temperature profile, T∝R, so that its initial Toomre stability profile is flat (i.e. constant over the entire disc) with Q = 2 (Fig. 1, dotted line). These discs were run for β = 5 and 6. Fig. 2 shows the images of the evolved disc with a decreasing Toomre stability profile and the flat Q disc, run with a cooling time, β = 5 (simulations Benchmark4 and Benchmark5, respectively). The first fragments form at ≈20 au in both discs irrespective of the initial temperature profile. Fig. 3 shows the final Toomre stability profiles of both discs run with β = 6 (simulations Benchmark1 and Benchmark6). Neither of these discs fragment, and both discs evolve into a steady state with very similar Toomre stability profiles. It can be seen that the change in the initial temperature profile does not make a difference to the final results since with β = 5 both discs fragment at the same radius and in the non-fragmenting cases, the final Toomre stability profiles are very similar. Since the temperature in a disc evolves, it is reassuring that the initial temperature profile of the disc does not play a part in the outcome.

Figure 1

Azimuthally averaged values of the Toomre parameter for the initial discs with a decreasing Toomre stability profile (simulations Benchmark1–4), set up in the same way as Rice et al. (2005, solid line), and with a flat Q profile with Q = 2 (simulations Benchmark5–7; dotted line). The critical value of Q_crit = 1 is also marked.

Open in new tabDownload slide

Figure 2

Surface mass density rendered image of the first fragments forming in the simulations using a decreasing Toomre stability profile (simulation Benchmark4, left-hand image) and a disc set up with a flat Q profile with Q = 2 (simulation Benchmark5, right-hand image). The discs were run with β = 5. In both cases, the discs first fragment at R_f≈ 20 au confirming that the initial temperature profile does not play a part in the evolution of the discs. The colour scale is a logarithmic scale ranging from log Σ=− 7 (dark) to −3 (light) M_⊙ au⁻².

Open in new tabDownload slide

Figure 3

Azimuthally averaged values of the Toomre parameter for the discs with initially decreasing (solid line) and flat (dotted line) Toomre stability profiles (simulations Benchmark1 and 6, respectively). The discs were run with β = 6. Despite having different initial temperature profiles, both discs reach a steady state with very similar Toomre stability profiles, confirming that the initial temperature does not play a part in the evolution of the discs. The critical value of Q_crit = 1 is also marked.

Open in new tabDownload slide

As mentioned in Section 1, current wisdom is that according to fragmentation theory, if the Toomre stability parameter is below unity and the time-scale on which the disc cools is faster than a critical value, then the disc should fragment. Therefore, if a disc was set up so that its initial Toomre stability profile was flat with Q = 1, it would be expected that fragments would form everywhere in the disc soon after the simulation is started. Fig. 4 shows the results of this simulation (Benchmark7). It can be seen that despite starting the simulation in a marginally stable state where any part of the disc may fragment soon after the simulation is started, the disc only fragments in the outer regions. This illustrates that the disc fragmenting in the outer regions cannot be related to the initial value of the Toomre stability profile, Q, and, more importantly, fragmentation cannot be a function of β alone.

Figure 4

Surface density rendered image of the fragmenting disc in simulation Benchmark7 with an initial flat Q profile with Q = 1. The simulation was run with β = 5. Despite the initial disc being in a state of marginal stability such that, in theory, any part of the disc may fragment, the disc only fragments in the outer regions. The colour scale is a logarithmic scale ranging from log Σ=− 8 (dark) to −3 (light) M_⊙ au⁻².

Open in new tabDownload slide

5 MAIN SIMULATIONS

In this section, we describe the initial conditions for all the individual numerical simulations we have performed to test our analytical predictions from Section 2. Table 2 provides a summary of the initial conditions as well as the radius at which the first fragment forms in the discs that do fragment.

Table 2

Open in new tab

Summary of the main simulations. p describes the initial surface mass density profile, Σ∝R^−p, and Σ₀ is the normalization constant required to produce a disc with mass M_disc. The final column represents the RHS of equation (7) for the location at which the first fragment forms, R_f. The simulations were run with an initial flat Toomre stability profile, Q.

Simulation name	β	p	Σ0[M⊙ (au)−2]	Mdisc (M⊙)	M⋆ (M⊙)	Mdisc/M⋆	Q	Disc radius (au)	Rf (au)
Reference-beta6	6	1	6.4 × 10−4	0.1	1	0.1	2	25	–	–
Reference-beta5.5	5.5	1	6.4 × 10−4	0.1	1	0.1	2	25	22	1.4 × 10−2
Reference-beta5	5	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta4	4	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta3	3	1	6.4 × 10−4	0.1	1	0.1	2	25	8	5.1 × 10−3
Reference-beta2	2	1	6.4 × 10−4	0.1	1	0.1	2	25	3	1.9 × 10−3
Reference-beta1	1	1	6.4 × 10−4	0.1	1	0.1	2	25	2.5	1.6 × 10−3
p1.5-beta4	4	1.5	1.8 × 10−3	0.1	1	0.1	2	25	–	–
p1.5-beta3.5	3.5	1.5	1.8 × 10−3	0.1	1	0.1	2	25	18	7.5 × 10−3
p1.5-beta3	3	1.5	1.8 × 10−3	0.1	1	0.1	2	25	1.7	2.3 × 10−3
p2-beta3.5	3.5	2	3.5 × 10−3	0.1	1	0.1	2	25	–	–
p2-beta3	3	2	3.5 × 10−3	0.1	1	0.1	2	25	0.45	3.5 × 10−3
p2-beta2	2	2	3.5 × 10−3	0.1	1	0.1	2	25	0.3	3.5 × 10−3
p2.5-beta4	4	2.5	4.4 × 10−3	0.1	1	0.1	2	25	–	–
p2.5-beta3.5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.4	7.0 × 10−3
p2.5-beta3.5-Q5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	5	25	0.3	8.1 × 10−3
p2.5-beta3	3	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.3	8.1 × 10−3
p2.5-beta2	2	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.35	7.5 × 10−3
p1-Mstar2	5	1	6.4 × 10−4	0.1	2	0.05	2	25	–	–
p1-Mstar0.5	5	1	6.4 × 10−4	0.1	0.5	0.2	2	25	13	1.7 × 10−2
p1-Mdisc0.2	5	1	1.3 × 10−3	0.2	1	0.2	2	25	14	1.8 × 10−2
p1-Mdisc0.05	5	1	3.2 × 10−4	0.05	1	0.05	2	25	–	–
p1-beta1-Mdisc0.01	1	1	6.4 × 10−5	0.01	1	0.01	2	25	6.5	4.2 × 10−4
p1-beta2-Mdisc0.01	2	1	6.4 × 10−5	0.01	1	0.01	2	25	15	9.6 × 10−4
p1-beta2.5-Mdisc0.01	2.5	1	6.4 × 10−5	0.01	1	0.01	2	25	17	1.1 × 10−3
p1-beta3-Mdisc0.01	3	1	6.4 × 10−5	0.01	1	0.01	2	25	–	–
p1-beta8-Mdisc0.3	8	1	1.9 × 10−3	0.3	1	0.3	2	25	–	–
p1-beta10-Mdisc0.5	10	1	3.2 × 10−3	0.5	1	0.5	2	25	–	–
p1-beta5-Mdisc1	5	1	6.4 × 10−3	1	1	1	2	25	5.5	3.5 × 10−2
p1-beta7-Mdisc1	7	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta10-Mdisc1	10	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta15-Mdisc1	15	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta6-extended	6	1	6.4 × 10−4	0.2	1	0.2	2	50	24.5	1.6 × 10−2
p1-beta7-extended	7	1	6.4 × 10−4	0.2	1	0.2	2	50	29	1.9 × 10−2
p1-beta8-extended	8	1	6.4 × 10−4	0.2	1	0.2	2	50	30	1.9 × 10−2
p1.5-beta4-extended	4	1.5	1.8 × 10−3	0.15	1	0.15	2	50	33	1.0 × 10−2
p1-beta6-Mdisc0.1-extended	6	1	3.2 × 10−4	0.1	1	0.1	2	50	40	1.3 × 10−2
p1-beta6-Mdisc0.2-Mstar2-extended	6	1	6.4 × 10−4	0.2	2	0.1	2	50	34	1.1 × 10−2

Simulation name	β	p	Σ0[M⊙ (au)−2]	Mdisc (M⊙)	M⋆ (M⊙)	Mdisc/M⋆	Q	Disc radius (au)	Rf (au)
Reference-beta6	6	1	6.4 × 10−4	0.1	1	0.1	2	25	–	–
Reference-beta5.5	5.5	1	6.4 × 10−4	0.1	1	0.1	2	25	22	1.4 × 10−2
Reference-beta5	5	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta4	4	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta3	3	1	6.4 × 10−4	0.1	1	0.1	2	25	8	5.1 × 10−3
Reference-beta2	2	1	6.4 × 10−4	0.1	1	0.1	2	25	3	1.9 × 10−3
Reference-beta1	1	1	6.4 × 10−4	0.1	1	0.1	2	25	2.5	1.6 × 10−3
p1.5-beta4	4	1.5	1.8 × 10−3	0.1	1	0.1	2	25	–	–
p1.5-beta3.5	3.5	1.5	1.8 × 10−3	0.1	1	0.1	2	25	18	7.5 × 10−3
p1.5-beta3	3	1.5	1.8 × 10−3	0.1	1	0.1	2	25	1.7	2.3 × 10−3
p2-beta3.5	3.5	2	3.5 × 10−3	0.1	1	0.1	2	25	–	–
p2-beta3	3	2	3.5 × 10−3	0.1	1	0.1	2	25	0.45	3.5 × 10−3
p2-beta2	2	2	3.5 × 10−3	0.1	1	0.1	2	25	0.3	3.5 × 10−3
p2.5-beta4	4	2.5	4.4 × 10−3	0.1	1	0.1	2	25	–	–
p2.5-beta3.5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.4	7.0 × 10−3
p2.5-beta3.5-Q5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	5	25	0.3	8.1 × 10−3
p2.5-beta3	3	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.3	8.1 × 10−3
p2.5-beta2	2	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.35	7.5 × 10−3
p1-Mstar2	5	1	6.4 × 10−4	0.1	2	0.05	2	25	–	–
p1-Mstar0.5	5	1	6.4 × 10−4	0.1	0.5	0.2	2	25	13	1.7 × 10−2
p1-Mdisc0.2	5	1	1.3 × 10−3	0.2	1	0.2	2	25	14	1.8 × 10−2
p1-Mdisc0.05	5	1	3.2 × 10−4	0.05	1	0.05	2	25	–	–
p1-beta1-Mdisc0.01	1	1	6.4 × 10−5	0.01	1	0.01	2	25	6.5	4.2 × 10−4
p1-beta2-Mdisc0.01	2	1	6.4 × 10−5	0.01	1	0.01	2	25	15	9.6 × 10−4
p1-beta2.5-Mdisc0.01	2.5	1	6.4 × 10−5	0.01	1	0.01	2	25	17	1.1 × 10−3
p1-beta3-Mdisc0.01	3	1	6.4 × 10−5	0.01	1	0.01	2	25	–	–
p1-beta8-Mdisc0.3	8	1	1.9 × 10−3	0.3	1	0.3	2	25	–	–
p1-beta10-Mdisc0.5	10	1	3.2 × 10−3	0.5	1	0.5	2	25	–	–
p1-beta5-Mdisc1	5	1	6.4 × 10−3	1	1	1	2	25	5.5	3.5 × 10−2
p1-beta7-Mdisc1	7	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta10-Mdisc1	10	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta15-Mdisc1	15	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta6-extended	6	1	6.4 × 10−4	0.2	1	0.2	2	50	24.5	1.6 × 10−2
p1-beta7-extended	7	1	6.4 × 10−4	0.2	1	0.2	2	50	29	1.9 × 10−2
p1-beta8-extended	8	1	6.4 × 10−4	0.2	1	0.2	2	50	30	1.9 × 10−2
p1.5-beta4-extended	4	1.5	1.8 × 10−3	0.15	1	0.15	2	50	33	1.0 × 10−2
p1-beta6-Mdisc0.1-extended	6	1	3.2 × 10−4	0.1	1	0.1	2	50	40	1.3 × 10−2
p1-beta6-Mdisc0.2-Mstar2-extended	6	1	6.4 × 10−4	0.2	2	0.1	2	50	34	1.1 × 10−2

Table 2

Open in new tab

Summary of the main simulations. p describes the initial surface mass density profile, Σ∝R^−p, and Σ₀ is the normalization constant required to produce a disc with mass M_disc. The final column represents the RHS of equation (7) for the location at which the first fragment forms, R_f. The simulations were run with an initial flat Toomre stability profile, Q.

Simulation name	β	p	Σ0[M⊙ (au)−2]	Mdisc (M⊙)	M⋆ (M⊙)	Mdisc/M⋆	Q	Disc radius (au)	Rf (au)
Reference-beta6	6	1	6.4 × 10−4	0.1	1	0.1	2	25	–	–
Reference-beta5.5	5.5	1	6.4 × 10−4	0.1	1	0.1	2	25	22	1.4 × 10−2
Reference-beta5	5	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta4	4	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta3	3	1	6.4 × 10−4	0.1	1	0.1	2	25	8	5.1 × 10−3
Reference-beta2	2	1	6.4 × 10−4	0.1	1	0.1	2	25	3	1.9 × 10−3
Reference-beta1	1	1	6.4 × 10−4	0.1	1	0.1	2	25	2.5	1.6 × 10−3
p1.5-beta4	4	1.5	1.8 × 10−3	0.1	1	0.1	2	25	–	–
p1.5-beta3.5	3.5	1.5	1.8 × 10−3	0.1	1	0.1	2	25	18	7.5 × 10−3
p1.5-beta3	3	1.5	1.8 × 10−3	0.1	1	0.1	2	25	1.7	2.3 × 10−3
p2-beta3.5	3.5	2	3.5 × 10−3	0.1	1	0.1	2	25	–	–
p2-beta3	3	2	3.5 × 10−3	0.1	1	0.1	2	25	0.45	3.5 × 10−3
p2-beta2	2	2	3.5 × 10−3	0.1	1	0.1	2	25	0.3	3.5 × 10−3
p2.5-beta4	4	2.5	4.4 × 10−3	0.1	1	0.1	2	25	–	–
p2.5-beta3.5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.4	7.0 × 10−3
p2.5-beta3.5-Q5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	5	25	0.3	8.1 × 10−3
p2.5-beta3	3	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.3	8.1 × 10−3
p2.5-beta2	2	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.35	7.5 × 10−3
p1-Mstar2	5	1	6.4 × 10−4	0.1	2	0.05	2	25	–	–
p1-Mstar0.5	5	1	6.4 × 10−4	0.1	0.5	0.2	2	25	13	1.7 × 10−2
p1-Mdisc0.2	5	1	1.3 × 10−3	0.2	1	0.2	2	25	14	1.8 × 10−2
p1-Mdisc0.05	5	1	3.2 × 10−4	0.05	1	0.05	2	25	–	–
p1-beta1-Mdisc0.01	1	1	6.4 × 10−5	0.01	1	0.01	2	25	6.5	4.2 × 10−4
p1-beta2-Mdisc0.01	2	1	6.4 × 10−5	0.01	1	0.01	2	25	15	9.6 × 10−4
p1-beta2.5-Mdisc0.01	2.5	1	6.4 × 10−5	0.01	1	0.01	2	25	17	1.1 × 10−3
p1-beta3-Mdisc0.01	3	1	6.4 × 10−5	0.01	1	0.01	2	25	–	–
p1-beta8-Mdisc0.3	8	1	1.9 × 10−3	0.3	1	0.3	2	25	–	–
p1-beta10-Mdisc0.5	10	1	3.2 × 10−3	0.5	1	0.5	2	25	–	–
p1-beta5-Mdisc1	5	1	6.4 × 10−3	1	1	1	2	25	5.5	3.5 × 10−2
p1-beta7-Mdisc1	7	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta10-Mdisc1	10	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta15-Mdisc1	15	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta6-extended	6	1	6.4 × 10−4	0.2	1	0.2	2	50	24.5	1.6 × 10−2
p1-beta7-extended	7	1	6.4 × 10−4	0.2	1	0.2	2	50	29	1.9 × 10−2
p1-beta8-extended	8	1	6.4 × 10−4	0.2	1	0.2	2	50	30	1.9 × 10−2
p1.5-beta4-extended	4	1.5	1.8 × 10−3	0.15	1	0.15	2	50	33	1.0 × 10−2
p1-beta6-Mdisc0.1-extended	6	1	3.2 × 10−4	0.1	1	0.1	2	50	40	1.3 × 10−2
p1-beta6-Mdisc0.2-Mstar2-extended	6	1	6.4 × 10−4	0.2	2	0.1	2	50	34	1.1 × 10−2

Simulation name	β	p	Σ0[M⊙ (au)−2]	Mdisc (M⊙)	M⋆ (M⊙)	Mdisc/M⋆	Q	Disc radius (au)	Rf (au)
Reference-beta6	6	1	6.4 × 10−4	0.1	1	0.1	2	25	–	–
Reference-beta5.5	5.5	1	6.4 × 10−4	0.1	1	0.1	2	25	22	1.4 × 10−2
Reference-beta5	5	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta4	4	1	6.4 × 10−4	0.1	1	0.1	2	25	20	1.3 × 10−2
Reference-beta3	3	1	6.4 × 10−4	0.1	1	0.1	2	25	8	5.1 × 10−3
Reference-beta2	2	1	6.4 × 10−4	0.1	1	0.1	2	25	3	1.9 × 10−3
Reference-beta1	1	1	6.4 × 10−4	0.1	1	0.1	2	25	2.5	1.6 × 10−3
p1.5-beta4	4	1.5	1.8 × 10−3	0.1	1	0.1	2	25	–	–
p1.5-beta3.5	3.5	1.5	1.8 × 10−3	0.1	1	0.1	2	25	18	7.5 × 10−3
p1.5-beta3	3	1.5	1.8 × 10−3	0.1	1	0.1	2	25	1.7	2.3 × 10−3
p2-beta3.5	3.5	2	3.5 × 10−3	0.1	1	0.1	2	25	–	–
p2-beta3	3	2	3.5 × 10−3	0.1	1	0.1	2	25	0.45	3.5 × 10−3
p2-beta2	2	2	3.5 × 10−3	0.1	1	0.1	2	25	0.3	3.5 × 10−3
p2.5-beta4	4	2.5	4.4 × 10−3	0.1	1	0.1	2	25	–	–
p2.5-beta3.5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.4	7.0 × 10−3
p2.5-beta3.5-Q5	3.5	2.5	4.4 × 10−3	0.1	1	0.1	5	25	0.3	8.1 × 10−3
p2.5-beta3	3	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.3	8.1 × 10−3
p2.5-beta2	2	2.5	4.4 × 10−3	0.1	1	0.1	2	25	0.35	7.5 × 10−3
p1-Mstar2	5	1	6.4 × 10−4	0.1	2	0.05	2	25	–	–
p1-Mstar0.5	5	1	6.4 × 10−4	0.1	0.5	0.2	2	25	13	1.7 × 10−2
p1-Mdisc0.2	5	1	1.3 × 10−3	0.2	1	0.2	2	25	14	1.8 × 10−2
p1-Mdisc0.05	5	1	3.2 × 10−4	0.05	1	0.05	2	25	–	–
p1-beta1-Mdisc0.01	1	1	6.4 × 10−5	0.01	1	0.01	2	25	6.5	4.2 × 10−4
p1-beta2-Mdisc0.01	2	1	6.4 × 10−5	0.01	1	0.01	2	25	15	9.6 × 10−4
p1-beta2.5-Mdisc0.01	2.5	1	6.4 × 10−5	0.01	1	0.01	2	25	17	1.1 × 10−3
p1-beta3-Mdisc0.01	3	1	6.4 × 10−5	0.01	1	0.01	2	25	–	–
p1-beta8-Mdisc0.3	8	1	1.9 × 10−3	0.3	1	0.3	2	25	–	–
p1-beta10-Mdisc0.5	10	1	3.2 × 10−3	0.5	1	0.5	2	25	–	–
p1-beta5-Mdisc1	5	1	6.4 × 10−3	1	1	1	2	25	5.5	3.5 × 10−2
p1-beta7-Mdisc1	7	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta10-Mdisc1	10	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta15-Mdisc1	15	1	6.4 × 10−3	1	1	1	2	25	–	–
p1-beta6-extended	6	1	6.4 × 10−4	0.2	1	0.2	2	50	24.5	1.6 × 10−2
p1-beta7-extended	7	1	6.4 × 10−4	0.2	1	0.2	2	50	29	1.9 × 10−2
p1-beta8-extended	8	1	6.4 × 10−4	0.2	1	0.2	2	50	30	1.9 × 10−2
p1.5-beta4-extended	4	1.5	1.8 × 10−3	0.15	1	0.15	2	50	33	1.0 × 10−2
p1-beta6-Mdisc0.1-extended	6	1	3.2 × 10−4	0.1	1	0.1	2	50	40	1.3 × 10−2
p1-beta6-Mdisc0.2-Mstar2-extended	6	1	6.4 × 10−4	0.2	2	0.1	2	50	34	1.1 × 10−2

We set up a series of Reference discs with M_disc = 0.1 M_⊙ consisting of 250 000 SPH particles, spanning 0.25 ≤R≤ 25 au, surrounding a 1 -M_⊙ star, modelled using a sink particle. The inner and outer radial disc boundaries have been set up in the same way as the benchmarking discs in Section 4. All the discs in this section have been set up with a flat Q profile. Therefore, as the surface mass density profile is varied, the initial temperature profile, q, is also varied accordingly, where T∝R^−q. The Reference discs have been set up so that Σ∝R⁻¹ and T∝R, normalized so that Q = 2. We highlight where we have differed from these initial conditions in the remaining simulations.

The analytical work presented in Section 2 suggests that for shallow surface mass density profiles, p < 2, fragments would form in the outer regions of the discs, whilst for discs with steeper surface mass density profiles, p > rsim 2, the disc would fragment in the inner regions. We therefore test different values of the slope of the surface mass density profiles using p = 1 (simulations Reference-beta5.5 and Reference-beta6), 1.5 (simulations p1.5-beta3.5 and p1.5-beta4), 2.0 (simulations p2-beta2, p2-beta3 and p2-beta3.5) and 2.5 (simulations p2.5-beta2, p2.5-beta3.5 and p2.5-beta4). In addition, we also carry out a further simulation which is the same as simulation p2.5-beta3.5 but with an initial flat Q profile, Q = 5 (i.e. so that the initial temperature is 25/4 times hotter than the disc in simulation p2.5-beta3.5), to test the effects of an initially hotter disc on the location of fragmentation (simulation p2.5-beta3.5-Q5).

The analysis also suggests that for a disc with a fast enough cooling time-scale such that it would fragment, the location at which the first fragment would form would move inwards to smaller radii as the cooling time-scale is decreased. We therefore test the effect of decreasing β on the fragment location by simulating the Reference disc (i.e. with a surface mass density profile, p = 1) with different values of the cooling time-scale, β = 5.5, 5, 4, 3, 2 and 1 (simulations Reference-beta5.5, Reference-beta5, Reference-beta4, Reference-beta3, Reference-beta2 and Reference-beta1, respectively).

We also argued that varying the disc or star mass would affect fragmentation. We test the effects of doubling and halving the star mass in simulations p1-Mstar2 and p1-Mstar0.5, respectively, and compare these to the Reference-beta5 simulation. We also carry out extensive tests of the effects of varying the disc mass first by doubling and halving the disc mass (simulations p1-Mdisc0.2 and p1-Mdisc0.05, respectively) and secondly by considering more extreme disc masses of 0.01 M_⊙ (simulations p1-beta0.3-Mdisc0.01, p1-beta1-Mdisc0.01, p1-beta2-Mdisc0.01, p1-beta2.5-Mdisc0.01 and p1-beta3-Mdisc0.01), 0.3 M_⊙ (simulation p1-beta8-Mdisc0.3), 0.5 M_⊙ (simulation p1-beta10-Mdisc0.5) and 1 M_⊙ (simulations p1-beta5-Mdisc1, p1-beta10-Mdisc1 and p1-beta15-Mdisc1) whilst maintaining a central star mass of 1 M_⊙.

The analytical work presented in Section 2 also showed that for a shallow surface mass density profile (p < 2), the radius of the disc may be the limiting aspect that causes a disc not to fragment. We therefore test a series of extended discs which have outer radii, R_out = 50 au. Simulations p1-beta6-extended, p1-beta7-extended and p1-beta8-extended are set up so that Σ₀ and p are the same as in the Reference discs (Fig. 5). However, to extend the disc to R_out = 50 au, the disc masses are increased to M_disc = 0.2 M_⊙. We run this simulation for β = 6, 7 and 8 (i.e. values that are larger than the critical values identified in Section 4). (Note that in order to keep the mass of the individual SPH particles the same as in the Reference simulations, we use 500 000 particles to make up this disc.) In addition we also set up an extended disc with a surface mass density profile, p = 1.5, which has a disc mass of 1.5 M_⊙ (so that Σ₀ is the same as in simulation p1.5-beta4). We run this simulation using β = 4 (simulation p1.5-beta4-extended). (As before, since we wish to keep the mass of the SPH particles the same as in simulation p1.5-beta4, we use 375 000 particles in this disc.)

Figure 5

Initial surface mass density profiles of the discs used in simulations p1-beta6 and p1-beta6-extended. The extended disc has the same surface mass density profile as the smaller disc.

Open in new tabDownload slide

Furthermore, we progress the analysis of extended discs by simulating two further discs (using 500 000 particles) as follows. The first is the same as that in p1-beta6-extended but using a total disc mass of M_disc = 0.1 M_⊙ (simulation p1-beta6-Mdisc0.1-extended) so that the total disc mass is the same as in p1-beta6 but Σ₀ is smaller. The second is also the same as in p1-beta6-extended but the central star mass is also M_⋆ = 2 M_⊙ so that the disc to star mass ratio is kept constant (simulation p1-beta6-Mdisc0.2-Mstar2-extended). Both of these discs are run with β = 6.

6 RESULTS

For the analysis that follows, the key aspect about the fragments that will be considered will be the first fragment that forms. This is because subsequent evolution of the disc following an initial fragmentation stage may involve additional complexities that are beyond the scope of this paper. Table 2 summarizes the key fragmenting results. The radius at which the first fragments form, R_f, has been determined by eye from the disc images. It is important to note that as seen in past simulations (e.g. Lodato & Rice 2004; Meru & Bate 2010), the surface mass density profile does not change significantly during the simulations, particularly for low-mass discs (M_disc≲ 0.2 M_⊙). We highlight where the surface mass density profiles do change significantly and discuss the effects of this.

6.1 Fragmentation dependency on the surface mass density profile

Figs 6–9 show that as the surface mass density profile steepens, the location at which the first fragment forms moves to smaller radii in the disc. The analytical theory presented in Section 2 shows that for a shallow surface mass density fall-off where p < 2, the fragments form in the outer regions of the disc (provided the cooling criterion is also satisfied). This is indeed the case for simulations with p = 1 and p = 1.5 (simulations Reference-beta5.5 and p1.5-beta3.5, respectively) as Figs 6 and 7 show that the fragments form at R_f≈ 20 and ≈19 au, respectively. Figs 10 and 11 show the radial profile of the aspect ratios (calculated by azimuthally averaging the sound speed at each radius and dividing by ΩR at that radius) in the discs compared to the RHS of equation (7) at a time shortly before the discs fragment. It can be seen that condition (7) is satisfied at the region in which the first fragment forms shortly after. The oscillations in H/R are due to temperature fluctuations since although the cooling rate in the disc changes smoothly with radius, the heating of the disc occurs primarily in the spiral shocks. This therefore confirms the analytical predictions for shallow surface mass density profiles (p < 2) presented in Section 2. It is important to note that for a flat Q profile, the temperature profiles in the discs are an increasing function of radius (T∝R) for p = 1 and a constant temperature profile for p = 1.5, yet the discs still fragment in the outer regions (again, re-emphasizing that the initial temperature profile does not play a part in the disc evolution).

Figure 6

Surface mass density rendered image of the fragmenting disc with initial surface mass density profile Σ∝R⁻¹. The simulation (Reference-beta5.5) used β = 5.5. The fragment forms in the outer regions of the disc, confirming the analytical predictions in Section 2. The colour scale is a logarithmic scale ranging from log Σ=− 6 (dark) to −3 (light) M_⊙ au⁻².

Open in new tabDownload slide

Figure 7

Surface mass density rendered image of the fragmenting disc with initial surface mass density profile Σ∝R^−3/2. The simulation used β = 3.5. The fragment forms in the outer regions of the disc, confirming the analytical predictions in Section 2. The colour scale is a logarithmic scale ranging from log Σ=− 7 (dark) to −2 (light) M_⊙ au⁻².

Open in new tabDownload slide

Figure 8

Surface mass density rendered image of the fragmenting disc with initial surface mass density profile Σ∝R⁻² (simulation p2-beta3). The simulation used β = 3. The fragment forms in the inner regions of the disc as shown by the zoomed-in image of the disc, confirming the analytical predictions in Section 2. The colour scale is a logarithmic scale ranging from log Σ=− 11 (dark) to 2 (light) M_⊙ au⁻² in the zoomed-out image and from log Σ=− 3.5 (dark) to 1 (light) M_⊙ au⁻² in the zoomed-in image.

Open in new tabDownload slide

Figure 9

Surface mass density rendered image of the fragmenting disc with initial surface mass density profile Σ∝R^−5/2 (simulation p2.5-beta3.5). The simulation used β = 3.5. The fragment forms in the inner regions of the disc, confirming the analytical predictions in Section 2. The colour scale is a logarithmic scale ranging from log Σ=− 12 (dark) to 3 (light) M_⊙ au⁻² in the zoomed-out image and from log Σ=− 4 (dark) to 0.4 (light) M_⊙ au⁻² in the zoomed-in image.

Open in new tabDownload slide

Figure 10

Plot of the disc aspect ratio, H/R (solid line), against the RHS of equation (7) (dotted line) for simulation Reference-beta5.5. Condition (7) is satisfied at ≈20 au where the disc first fragments, confirming the analytical predictions in Section 2.

Open in new tabDownload slide

Figure 11

Plot of the disc aspect ratio, H/R (solid line), against the RHS of equation (7) (dotted line) for simulation p1.5-beta3.5. Condition (7) is satisfied at ≈19 au where the disc first fragments, confirming the analytical predictions in Section 2.

Open in new tabDownload slide

The analytical theory for p > rsim 2 suggested that if the disc was to fragment, it would do so in the inner regions of the disc. Figs 8 and 9 show that this is indeed the case. Fig. 12 shows that the analytical condition is just satisfied for simulation p2-beta3 at the location at which the fragment forms. Fig. 13 shows that for simulation p2.5-beta3.5, the analytical condition is also satisfied at the location at which fragmentation occurs soon after.

Figure 12

Plot of the disc aspect ratio, H/R (solid line), against the RHS of equation (7) (dotted line) for simulation p2-beta3 for the radial range of the entire disc (upper panel) as well as zoomed into the inner regions (lower panel). Condition (7) is marginally satisfied at ≈0.4 au where the disc first fragments, confirming the analytical predictions in Section 2.

Open in new tabDownload slide

Figure 13

Plot of the disc aspect ratio, H/R (solid line), against the RHS of equation (7) (dotted line) for simulation p2.5-beta3.5 for the radial range of the entire disc (upper panel) as well as zoomed into the inner regions (lower panel). Condition (7) is satisfied at ≈0.4 au where the disc first fragments, confirming the analytical predictions in Section 2.

Open in new tabDownload slide

We therefore show that as the surface mass density profile is steepened so that more of the mass is concentrated in the inner regions of the disc, fragmentation moves towards smaller radii. It is important to note that the trend that fragmentation moves to smaller radii for steeper surface mass density profiles is valid even when considering a uniform value of β (compare simulations Reference-beta3, p1.5-beta3, p2-beta3 and p2.5-beta3 which are run using β = 3 and fragment at ≈8, 1.7, 0.45 and 0.3 au, respectively).

In addition, the results summarized in Table 2 show that a single value of β_crit is not applicable over all surface mass density profiles since the minimum value of β that a disc can have without fragmenting varies with the surface mass density profile.

Simulation p2.5-beta3.5-Q5 was the same as simulation p2.5-beta3.5 but had an initial disc that was hotter by a factor of 25/4. The results of this simulation show that the disc still fragmented in the inner regions.

6.2 Effect of the cooling time-scale, β, on the fragment location

In Section 2 the analytical work presented suggested that for a fragmenting disc with a shallow surface mass density profile (p < 2), a decrease in the value of β would cause the location at which the first fragment forms to move inwards to smaller radii.

Fig. 14 shows the radius at which the first fragment forms for different values of β (simulations Reference-beta5.5, Reference-beta5, Reference-beta4, Reference-beta3, Reference-beta2 and Reference-beta1). We can see a clear trend showing that the radius of fragmentation moves inwards for more efficient cooling.

Figure 14

The radius at which the first fragment forms in the Reference simulations. The discs in these simulations are identical with a surface mass density profile, p = 1, but were run with different values of the cooling time-scale in units of the orbital time-scale, β. The radius at which the first fragment forms moves inwards with more efficient cooling.

Open in new tabDownload slide

6.3 The influence of star mass on fragmentation

In Section 2, we showed that decreasing the star mass is more likely to cause conditions (6) and (7) to be satisfied over a larger part of the disc, and hence the disc is more likely to fragment. We test three identical discs with star masses of 0.5, 1 and 2 M_⊙ (simulations p1-Mstar0.5, Reference-beta5 and p1-Mstar2, respectively) which are run with the same cooling time-scale, β = 5. It can be seen from Table 2 that the Reference-beta5 disc first fragments at R_f≈ 20 au. However, when the star mass is halved, the first fragment forms at R_f≈ 13 au. Since the RHS of condition (7) ∝R²/M_⋆, if the star mass is halved for the same value of β (and hence the same value of the LHS of condition 7), the radius at which the first fragment forms, R_f, decreases by a factor of ⁠. Conversely, doubling the star mass makes it harder for the condition to be satisfied and indeed the disc in simulation p1-Mstar2 does not fragment. Instead, it settles into a state of marginal stability with Q≈ 1.

6.4 The influence of disc mass on fragmentation

The analytics presented in Section 2 showed that increasing the disc mass (and hence increasing Σ₀) allows conditions (6) and (7) to be satisfied over a larger part of the disc and hence the disc is more likely to fragment. We initially test this by comparing the results of simulations p1-Mdisc0.05, Reference-beta5 and p1-Mdisc0.2 which are identical discs except that the disc masses are 0.05, 0.1 and 0.2 M_⊙, respectively. Table 2 shows that doubling the disc mass from 0.1 to 0.2 M_⊙ does indeed cause the fragmentation conditions derived in Section 2 to be satisfied over a larger portion of the disc since the first fragments form at R_f≈ 20 and ≈14 au, respectively. Since the RHS of condition (7) ∝ΣR², if the disc mass (and hence Σ) is doubled for the same value of β (and hence the same value of the LHS of condition 7), the radius at which the first fragment forms, R_f, decreases by a factor of ⁠. However, halving the disc mass makes it harder for the conditions to be satisfied and consequently the disc does not fragment.

In addition, we also simulate a very low-mass disc (M_disc = 0.01 M_⊙) and found that it fragments if β = 1, 2 and 2.5 but not for β = 3 (simulations p1-beta1-Mdisc0.01, p1-beta2-Mdisc0.01, p1-beta2.5-Mdisc0.01 and p1-beta3-Mdisc0.01, respectively). As β increases, the fragment location moves out in the disc, as found in Section 6.2. It is clear that, as with varying the star mass, the disc mass plays a crucial role in the fragmentation and the condition for fragmentation cannot simply be described using a single critical value of the cooling time-scale.

In addition, we simulate higher mass discs, M_disc = 0.3 and 0.5 M_⊙, which are run using β = 8 and 10, respectively (simulations p1-beta8-Mdisc0.3 and p1-beta10-Mdisc0.5, respectively) as well as discs with M_disc = 1 M_⊙, which we simulate using β = 5, 10 and 15 (simulations p1-beta5-Mdisc1, p1-beta10-Mdisc1 and p1-beta15-Mdisc1, respectively). We find that with the exception of simulation p1-beta5-Mdisc1, the discs do not fragment. Fig. 15 shows the surface mass density profile of simulation p1-beta7-Mdisc1 at the start and more than four ORPs after the start of the simulation. It can be seen that unlike the lower mass discs, the profile steepens and the value of Σ₀ increases. This is the case for all the non-fragmenting high-mass disc simulations. Fig. 16 shows the plot of the aspect ratio profile of this simulation against the RHS of condition (7) and shows that condition (7) is just satisfied in the inner regions. However, since during the simulation the disc mass redistributes itself, the surface mass density profile changes and consequently, using the newly obtained values of Σ₀ and p, condition (7) is not satisfied. We note that the only high-mass disc that does fragment (simulation p1-beta5-Mdisc1) does so because the cooling time is so rapid that fragmentation occurs before the disc has had the chance to restructure itself.

Figure 15

Surface mass density profiles for simulation p1-beta7-Mdisc1 at the start (solid line) and at a time more than four ORPs later (dotted line). Unlike the low-mass simulations whose surface mass density profiles do not change throughout the simulations, the profile for this disc steepens causing a change in the effective values of Σ₀ and p.

Open in new tabDownload slide

Figure 16

Plot of the disc aspect ratio, H/R (solid line), for simulation p1-beta7-Mdisc1, against the RHS of equation (7) using the initial values of Σ₀ and p (dashed line) and the new values of Σ₀ and p (dotted line) determined after the disc has evolved for greater than ORPs by which time its surface mass density profile has changed. The condition is satisfied using the initial values of Σ₀ and p but not using the new values, and hence the disc does not fragment.

Open in new tabDownload slide

6.5 The role of the disc radius on fragmentation

In Section 2, we showed that for shallow surface mass density profiles (p < 2), fragmentation might occur for any value of the cooling time-scale, if the disc is large enough.

Simulation Reference-beta6 (a disc with R_out = 25 au) did not fragment and though we did not run the same simulation with β = 7 or 8, we would expect that they would also not fragment. However, extended discs with the same values of Σ₀ and p as Reference-beta6 do indeed fragment for β = 6 (simulation p1-beta6-extended), β = 7 (p1-beta7-extended) and β = 8 (p1-beta8-extended; Fig. 17) with the first fragments forming at R_f≈ 25, 29 and 30 au, respectively. Similarly, we also simulate an extended disc with p = 1.5 and β = 4 (simulation p1.5-beta4-extended) and show that whilst the same disc truncated at R_out = 25 au does not fragment, the extended disc does indeed fragment (at R_f≈ 33 au).

Figure 17

Surface mass density rendered image of the fragmenting disc in simulation p1-beta8-extended with initial surface mass density profile Σ∝R⁻¹, but extending to 50 au rather than 25 au. This simulation was run with β = 8. According to Rice et al. (2005), this disc should not fragment since the cooling time-scale β is larger than the critical value previously obtained with a radius of 25 au. This simulation shows that the fragmentation criterion is more complex than a single critical cooling parameter. The colour scale is a logarithmic scale ranging from log Σ=− 8 (dark) to −2 (light) M_⊙ au⁻².

Open in new tabDownload slide

In addition, given that in Section 6.4 we showed that the disc mass plays a part in whether fragmentation occurs or not, we simulate a 0.1 -M_⊙ disc which extends to R_out = 50 au (simulation p1-beta6-Mdisc0.1-extended). The surface mass density profile is the same as in simulation p1-beta6-extended (p = 1) but Σ₀ is decreased. The results show that the disc fragments at au (cf. ≈25 au for simulation p1-beta6-extended). Therefore, whilst the disc mass affects where in the disc the first fragment forms, the conclusion that a small disc which does not fragment for a particular value of β may fragment at larger radii for the same value of β is still valid. Furthermore, we can see that if the disc to star mass ratio is kept constant at M_disc/M_⋆ = 0.1 for an extended disc, the disc fragments (at R_f≈ 34 au; simulation p1-beta6-Mdisc0.2-Mstar2-extended) whilst the small disc with the same disc to star mass ratio does not fragment, further corroborating the fact that the radius of the disc is important.

7 DISCUSSION

It has previously been accepted that for a self-gravitating disc whose only source of heating is internally generated from the gravitational instability, the disc will fragment if the cooling time-scale is short enough (Gammie 2001). However, we find that fragmentation at a given radius is not only dependent on the cooling time-scale, β, but also on the disc surface density (i.e. disc mass and profile) and the star mass.

This is in contrast to Rice et al. (2005) who suggested that the fragmentation criterion is independent of the disc mass though in agreement with Rice et al. (2003) who found that for a higher disc mass, fragmentation was easier: using β = 5, they found that a 0.25 -M_⊙ disc fragmented whilst a 0.1 -M_⊙ disc did not, but instead required a lower value of β. In particular, we highlight that in the past, it has been thought that a massive disc is required for fragmentation to occur. However, we show that it is indeed possible for low-mass discs [M_disc∼O(0.01) M_⊙] to fragment if the cooling in the discs is rapid enough. On the other hand, for high-mass discs (M_disc≥ 0.3 M_⊙ within R_out = 25 au), the discs do not fragment, unless the cooling time is fast, due to a steepening of the surface mass density profile and an increase in Σ₀ making condition (7) harder to satisfy. Lodato & Rice (2005) also find that with larger disc masses, matter is redistributed causing the surface mass density profile to steepen.

In addition, we find that the critical value of β found for one particular surface mass density profile is not applicable to another disc mass distribution. Our simulations show that for a steeper surface mass density profile, the cooling time-scale required for a disc to fragment is smaller. Cossins et al. (2010) found that for a disc with Σ∝R^−3/2, β_crit∼ 4 –4.5 whereas Rice et al. (2005) found that β_crit∼ 6 –7 for a disc with Σ∝R⁻¹. In addition, Rice et al. (2003) found that for a surface mass density profile of Σ∝R^−7/4, the fragmentation boundary for a 0.1 -M_⊙ disc around a 1 -M_⊙ star is between β = 3 and 5 which Rice et al. (2005) note is inconsistent with their results. However, our results explain these previous inconsistencies present in the literature.

We also find that for p < 2, if a disc does not fragment for a particular cooling time-scale in units of the orbital time-scale, β, a larger disc with the same surface mass density profile may well fragment (compare simulations Reference-beta6 and p1-beta6-extended or p1.5-beta4 and p1.5-beta4-extended, and also p1-beta7-extended or p1-beta8-extended). Therefore, a critical cooling time-scale can only be specified for a particular surface mass density at a particular disc radius and can therefore not be a general rule. The previous fragmentation criterion found by Rice et al. (2005) was for a disc to star mass ratio of 0.1 with M_⋆ = 1 M_⊙, a surface mass density profile of Σ∝R⁻¹ and R_out = 25 au.

Finally, we also show that the mass of the star plays a part in whether a disc fragments or not. Therefore, as shown by the RHS of condition (7), whether a disc fragments or not is essentially a ‘trade-off’ between the local surface mass density and the star mass.

7.1 The link between β, M_disc, M_★, the surface mass density profile, p, and fragmentation

Section 6 shows that the fragmentation criterion is clearly a complex problem which cannot simply depend on a single critical cooling time-scale as has previously been thought to be the case. Equation (7) and the results presented here clearly show that there is a link between the cooling time-scale (in terms of β), the disc mass (or, more accurately, the local surface mass density) and the star mass.

Such a link can be explained physically. As a disc cools, gravitational instability develops resulting in density fluctuations above and below the unperturbed density, δρ/ρ. The spiral structures involve shocks which produce heat that may balance the disc's cooling, thus reaching an equilibrium state. If the disc mass was irrelevant, for any particular star mass, and β, one would expect the fluctuations, δρ/ρ, to be the same in all discs with the same surface density profile. However, comparing a low-mass disc with a high-mass disc with the same relative density fluctuations (δρ/ρ), the density enhancement, δρ, in the higher mass disc will clearly be greater. At some disc mass, this enhancement will be self-gravitating (i.e. it will become a fragment), while in a lower mass disc the density fluctuation will not form a fragment (unless the value of β is lowered). Similarly, if the disc is kept the same and the star's mass is increased, a given density enhancement may be sheared apart by the differential rotation so that a low value of β will be required for fragmentation.

Fig. 18 shows a graph of β against ΣR²_f/M_⋆, where R_f is the radius at which the first fragment forms, and includes all the fragmenting simulations presented in this paper (including the Benchmarking simulations). When interpreting these results, one should also note that an increase in ΣR²/M_⋆ does not necessarily imply an increased disc to star mass ratio: it is possible to have a low disc to star mass ratio but using a steeper surface mass density profile. We can see that there is a clear trend that as the RHS of equation (7) increases, so too does the value of β that will allow a fragment to form. The trendline presents the relation

9

where δ≈ 1/2 and the constant of proportionality η≈ 47 which we find using a least-squares fit. We have performed a few calculations with higher resolution (2 × 10⁶ SPH particles) and find that this trend of β increasing with ΣR²/M_⋆ is maintained, although the exact values of η and δ may change slightly if we were able to perform all the calculations with high resolution.

Figure 18

Logarithmic graph showing the trend between β and ΣR²_f/M_⋆ determined by considering the location at which the first fragment forms in the discs, R_f. The results include those simulations with a surface mass density profile, p = 1 (filled triangles), p = 1.5 (open triangles), p = 2 (open squares) and p = 2.5 (crosses). It is clear that a single critical value of β is not the case for all discs and that there is a relation between β, M_disc, M_⋆ and the surface mass density profile, p, that determines whether fragmentation occurs or not. The trendline has been determined by considering discs with shallow surface mass density profiles, p < 2, only as those discs with p > rsim 2 will always fragment in the innermost regions first. The grey shaded region is where we expect that subsequent fragmentation may take place in discs with p < 2.

Open in new tabDownload slide

It is also important to note that the trendline has only been produced using the results of simulations with p < 2. This is because for p > rsim 2, the disc will always fragment in the innermost regions first. Consequently, if the results for p > rsim 2 were included, this would cause the trendline to be somewhat skewed.

The trendline can give some very useful information for those simulations with p < 2. The grey region is where we predict that subsequent fragmentation is feasible. Traversing the plot in a vertical direction downwards from the trendline into the grey region, at any particular value of ΣR²/M_⋆ the disc will fragment at all values of β less than the limit given by the trendline (though the radius being considered will not necessarily be the first location at which fragmentation occurs). Similarly, traversing the plot in a horizontal direction from the trendline to the right-hand side of the plot into the grey region, for a particular value of β, fragmentation is possible at a particular radius if the disc mass is increased or star mass is decreased. Similarly, for a particular value of β and any one combination of the disc to star mass ratio, fragmentation is also possible at larger radii than the location specified by the trendline, i.e. to the right-hand side of the trendline. Therefore, the trendline predicts the minimum possible radius at which fragments could in theory form in discs with shallow surface mass density profiles, p < 2.

For discs with p > rsim 2 which fragment in the inner regions, we expect that subsequent fragmentation may take place further out in the disc as far out as given by the trendline in Fig. 18. In other words, for these discs, we expect there to be a maximum radius outside of which fragmentation will not occur (since the surface mass density fall-off is steep, the outer regions of these discs may struggle to have enough mass for gravitational instability to be significant). However, since we stop the simulations soon after the first fragment forms due to the increased computational resources required to follow the simulations further, we are unable to test this prediction. Further work needs to be done in this area and is beyond the scope of this paper. We also note that in real discs, for fragmentation to occur in the inner regions, the cooling time would have to be very small since the dynamical time at small radii would also be very small. Such short cooling times may not be possible in real discs.

It is also important to note that these simulations have been carried out using a ratio of specific heats, γ = 5/3. As shown by Rice et al. (2005), the ratio of specific heats plays a key role in the fragmentation boundary. We therefore also anticipate a further dependency on the equation of state. Furthermore, for high-mass discs (M_disc > rsim 0.3 M_⊙), we find that the initial surface mass density conditions (i.e. Σ₀ and p) cannot be used to determine whether the disc will fragment or not. This is because as the disc evolves, the surface mass density profile steepens causing Σ₀ and p to change. Consequently, though parts of the disc may start off in the grey shaded region of Fig. 18 and hence may be expected to fragment, the disc may restructure itself on a time-scale faster than the cooling time-scale such that it moves out of this region and hence does not fragment.

Following the work of Gammie (2001) and Rice et al. (2005), a number of authors have used the critical cooling time-scale, β_crit (or gravitational stress, α_GI,max), to predict fragmentation in realistic discs (e.g. Clarke 2009; Rafikov 2009; Cossins et al. 2010; Kratter et al. 2010). In light of the new results presented here, we would encourage previous conclusions based on these critical values to be revisited.

8 CONCLUSIONS

We present an analytical approach to examine the fragmentation of self-gravitating protoplanetary discs and confirm the results using global three-dimensional numerical simulations. Our key result is that fragmentation does not simply depend on the disc cooling time-scale, β, but also on the ratio of the surface mass density at radius, R, to the stellar mass, i.e. ΣR²/M_⋆. We find that fragmentation occurs when

10

where δ≈ 1/2 and η≈ 47. For a power-law surface mass density, Σ∝R^−p, this relation predicts the innermost radii at which subsequent fragments can form in a disc with shallow surface mass density profiles (p < 2) as well as the radius outside of which fragments cannot form in discs with steep surface mass density profiles (p > rsim 2). Generally, we find that an increase in the steepness of the disc surface mass density profile promotes fragmentation at smaller radii.

We thank the anonymous referee for their thoughtful comments on the paper. We also thank Jim Pringle and Daniel Price for useful discussions. The calculations reported here were performed using the University of Exeter's SGI Altix ICE 8200 supercomputer and Intel Nehalem (i7) cluster. The disc images were produced using SPLASH (Price 2007). MRB is grateful for the support of a EURYI Award which also funded FM. This work, conducted as part of the award ‘The formation of stars and planets: Radiation hydrodynamical and magnetohydrodynamical simulations’ made under the European Heads of Research Councils and European Science Foundation EURYI (European Young Investigator) Awards scheme, was supported by funds from the participating organizations of EURYI and the EC Sixth Framework Programme.

REFERENCES