Protostellar discs formed from rigidly rotating cores

Abstract

1 INTRODUCTION

Stars form from the gravitational collapse of pre-stellar cores. We expect that, to a first approximation, angular momentum is conserved minutely during the early stages of collapse, and this leads to the formation of a protostellar disc and/or to fragmentation into a multiple system. Thus, understanding the formation and evolution of protostellar discs is a critical stage linking the initial conditions within pre-stellar cores to the final outcome of the star formation process.

In this paper, we simulate the collapse of isolated, rigidly rotating cores, using the three-dimensional OpenMP-parallelized smoothed particle hydrodynamics (SPH) code vine (Nelson, Wetzstein & Naab 2009; Wetzstein et al. 2009). All cores have the same mass and the same initial density profile, given by a Bonnor–Ebert sphere (BES) (Ebert 1955; Bonnor 1956); they are distinguished only by their angular momentum. We study the influence of core angular momentum on the properties of the protostellar discs and protostars that condense out of such discs. In a companion paper (Walch et al. 2009), we explore and contrast the properties of protostellar discs and protostars condensing out of turbulent cores having the same net angular momenta as the rigidly rotating cores considered here.

Previous numerical work on the collapse and fragmentation of cores has demonstrated that the critical ingredients are the net amount of angular momentum and its distribution, the initial density profile and the treatment of thermodynamics (for a summary see Goodwin et al. 2006). For example, Miyama, Hayashi & Narita (1984) and Burkert & Bodenheimer (1993) examine the evolution of cores with a uniform initial density and an isothermal equation of state. They find that the outcome depends on the parameters α≡U_THERM/|U_GRAV| and β≡U_ROT/|U_GRAV|, where U_THERM, U_GRAV and U_ROT are the thermal, gravitational and rotational energies of the initial core. Specifically, Miyama et al. (1984) only find fragmentation for αβ < 0.12.

However, observations of pre-stellar cores suggest that they are centrally condensed and well represented by a Bonnor–Ebert profile (Alvés, Lada & Lada 2001; André et al. 2004), which is the equilibrium configuration for an isothermal cloud. Centrally condensed clouds appear to be more stable against fragmentation than cores with a uniform density (Lynden-Bell 1964; Boss 1993; Burkert & Bodenheimer 1996). For a review see Bodenheimer et al. (2000).

To date, only the parameter studies of the collapse of BES have been performed using grid codes. Matsumoto & Hanawa (2003) investigate the collapse of 3 M_⊙ critical BES with the density enhanced by 10 per cent. They invoke a barotropic equation of state and investigate different total angular momenta, different rotation laws and bar mode perturbations of different strengths. Their simulations use a nested grid with mirror symmetry about the equatorial plane, and fine resolution in the centre (the location of the primary protostar) but quite coarse resolution in the outer parts (the disc region), particularly in the vertical direction. They conclude that fragmentation requires Ω_Ct_FF > 0.05, where Ω_C and t_FF are the initial angular speed and freefall time at the centre of the core; this is equivalent to a lower limit on β. The study of disc fragmentation by Nelson (2006) suggests that these results may be compromised by poor vertical resolution.

Banerjee, Pudritz & Holmes (2004) investigate the collapse of a 2.1 M_⊙ core and a 168 M_⊙ core, using the flash amr code. The cores are initially BES with the density enhanced everywhere by 10 per cent, and then a 10 per cent m= 2 bar perturbation superimposed; they are initially in rigid rotation. In order to capture the thermodynamics more realistically, they use the tabulated molecular-line cooling prescription of Neufeld, Lepp & Melnick (1995). They perform only one simulation of a low-mass core, which has Ω_Ct_FF= 0.2 and produces a bar; they suggest that this bar might subsequently fragment if the simulation was followed further.

Boss et al. (2000) have demonstrated that core fragmentation is highly sensitive to the treatment of thermodynamics. However, although Krumholz, Klein & McKee (2007) have performed three-dimensional AMR simulations of 100 M_⊙ and 200 M_⊙ collapsing cores with radiative transfer, there has as yet been no comprehensive parameter study of low-mass collapsing cores using a realistic treatment of the gas thermodynamics. This paper is an attempt to perform such a study.

In Sections 2 and 3, we describe the initial conditions, the numerical code and the constitutive physics. In Sections 4 and 5, we present in detail the evolution of cores with, respectively, low and high Ω_O and in Section 6 we describe how the evolution changes as Ω_O is increased between these limits. In the Appendix, we develop an analytic model, and derive a condition for fragmentation based on the initial parameters of a core. In Section 7, this criterion is applied to a sample of observed cores, and it is shown that they are all stable against disc fragmentation. In Section 8, we summarize our main conclusions.

2 INITIAL CONDITIONS

2.1 The initial core density profile

The observed column-density profiles of isolated pre-stellar cores are often well matched by approximately critical BES, i.e. equilibrium isothermal gas clouds which are close to the critical state for gravitational condensation (see e.g. Alvés et al. 2001; Hennebelle et al. 2003). We have therefore modelled our initial cores as marginally supercritical BES, truncated at dimensionless radius ξ_B= 6.9 (the critical value is ξ_B= 6.45) and then with their density increased everywhere by 10 per cent, to ensure immediate contraction. The gas is assumed to be molecular hydrogen at 28 K, and therefore the isothermal sound speed is a_O= 0.34 km s⁻¹. The central density is ρ_C= 10⁻¹⁸ g cm⁻³, and so the mass and radius are

1

2

Here, ψ(ξ) is the Isothermal Function and ψ′ (ξ) its first derivative (see Chandrasekhar & Wares 1949), ψ_B≡ψ(ξ_B) and ⁠. The core must be contained by an external pressure

3

delivered by a hot rarefied medium. The freefall time at the core centre (hence the time-scale on which the central primary protostar starts to condense out) is t_FF= 67 kyr. The thermal energy of the core is U_THERM= 3M_Oa²_O/2, and the self-gravitational potential energy is

Hence, the ratio of thermal to gravitational energy is

4

we note that α > 1/2, because of the external pressure term in the Virial Theorem.

2.2 The initial velocity field

The only core property which we vary is the initial angular speed, Ω_O. For a core in rigid rotation about the unit vector ⁠, the initial velocity field is then

5

The moment of inertia of the core is

6

Hence, the total specific angular momentum is given by

7

and the initial ratio of rotational to gravitational energy by

8

Our choice of representative values of Ω_O is based on two constraints. (i) The total specific angular momenta in low-mass cores are typically of the order of 10²¹ cm² s⁻¹ (Goodman et al. 1993). (ii) The non-thermal velocity dispersions in low-mass cores are usually subsonic and at most transonic (e.g. Barranco & Goodman 1998; André et al. 2007). Therefore, the rotational velocity at the edge of the core should not be much greater than a_O. The values of Ω_O that we consider are listed in Table 1; the equator of the fastest spinning core rotates at Mach 1.5.

Table 1

Open in new tab

Column 1 gives the ID of the simulation, and Column 2 gives the initial rotational frequency, Ω_O. Column 3 gives the time at which the first protostellar object forms, t_O, and Column 4 gives the time at which 28 per cent of the mass is either in the protostellar disc or in protostars, t₂₈. Column 5 gives the mass of the primary protostar at this juncture, M_★, and Columns 6 and 7 give the mass, M_D, and radius, R_D, of the protostellar disc at the same time. Column 8 indicates whether the disc fragments.

Run	ΩO	tO	t28	M★	MD	RD	F
1	0.6	94.5	114.8	0.53	1.17	250	N
2	0.8	94.5	115.4	0.31	1.39	344	N
3	1.0	96.5	118.4	0.25	1.45	433	N
4	1.2	96.5	119.2	0.23	1.47	567	N
5	1.5	98.5	123.3	0.17	1.53	739	Y
6	1.6	102.4	126.7	0.24	1.46	1091	Y
7	1.8	104.4	130.4	0.16	1.54	952	Y

Run	ΩO	tO	t28	M★	MD	RD	F
1	0.6	94.5	114.8	0.53	1.17	250	N
2	0.8	94.5	115.4	0.31	1.39	344	N
3	1.0	96.5	118.4	0.25	1.45	433	N
4	1.2	96.5	119.2	0.23	1.47	567	N
5	1.5	98.5	123.3	0.17	1.53	739	Y
6	1.6	102.4	126.7	0.24	1.46	1091	Y
7	1.8	104.4	130.4	0.16	1.54	952	Y

Table 1

Open in new tab

Column 1 gives the ID of the simulation, and Column 2 gives the initial rotational frequency, Ω_O. Column 3 gives the time at which the first protostellar object forms, t_O, and Column 4 gives the time at which 28 per cent of the mass is either in the protostellar disc or in protostars, t₂₈. Column 5 gives the mass of the primary protostar at this juncture, M_★, and Columns 6 and 7 give the mass, M_D, and radius, R_D, of the protostellar disc at the same time. Column 8 indicates whether the disc fragments.

Run	ΩO	tO	t28	M★	MD	RD	F
1	0.6	94.5	114.8	0.53	1.17	250	N
2	0.8	94.5	115.4	0.31	1.39	344	N
3	1.0	96.5	118.4	0.25	1.45	433	N
4	1.2	96.5	119.2	0.23	1.47	567	N
5	1.5	98.5	123.3	0.17	1.53	739	Y
6	1.6	102.4	126.7	0.24	1.46	1091	Y
7	1.8	104.4	130.4	0.16	1.54	952	Y

Run	ΩO	tO	t28	M★	MD	RD	F
1	0.6	94.5	114.8	0.53	1.17	250	N
2	0.8	94.5	115.4	0.31	1.39	344	N
3	1.0	96.5	118.4	0.25	1.45	433	N
4	1.2	96.5	119.2	0.23	1.47	567	N
5	1.5	98.5	123.3	0.17	1.53	739	Y
6	1.6	102.4	126.7	0.24	1.46	1091	Y
7	1.8	104.4	130.4	0.16	1.54	952	Y

3 NUMERICAL METHOD AND CONSTITUTIVE PHYSICS

3.1 The vine code

The simulations have been performed with the Tree-SPH code vine. A full description of vine can be found in Wetzstein et al. (2009) and Nelson et al. (2009). vine is parallelized with OpenMP directives. It invokes a leapfrog integrator, and individual particle time-steps with a Courant-Friedrichs-Lewy (CFL) tolerance parameter of 0.1. Gravitational accelerations are estimated using a tree, with opening angle θ= 0.005 (Springel, Yoshida & White 2001). The gravitational softening length equals the hydrodynamical smoothing length (Bate & Burkert 1997), which is adapted so that each particle has neighbours. Hydrodynamical forces are treated with periodic boundary conditions, but gravitational forces are not. Artificial viscosity is treated using the standard prescription of Gingold & Monaghan (1983) with α_AV= 1 and β_AV= 2, plus the Balsara switch (Balsara 1995). We have experimented with time-varying viscosity (Morris & Monaghan 1997), but find that it makes no significant difference.

3.2 Equation of state and molecular-line cooling

For the purpose of calculating the gross thermodynamics, we assume that the gas is pure molecular hydrogen, with ratio of specific heats γ= 7/5 and molecular weight μ= 2; thus for simplicity we are neglecting the fact that the rotational degrees of freedom of H₂ are not fully excited at low temperatures, and we are neglecting the contribution of helium to μ. The equation of state is

9

and the isothermal sound speed is

10

In solving the energy equation we adopt the following procedure. At low densities, ρ < ρ_CRIT= 10⁻¹³ g cm⁻³, we use the cooling rates computed by Neufeld et al. (1995), with the constraint that the temperature is not allowed to fall below T= 9 K. At high densities, ρ > ρ_CRIT, we switch off the radiative cooling, and the gas then evolves adiabatically. Neufeld et al. (1995) have computed the equilibrium abundances of key molecules and atoms (CO, H₂O, H₂, HCl, O₂, C and O), and the resulting net cooling rate due to line emission, Λ, as a function of density and temperature, using the local velocity gradient method to treat optical-depth effects. Neufeld has kindly supplied these rates in the form of a look-up table for densities in the range and temperatures in the range 10 ≤ (T/K) ≤ 2500. Most previous work in this field has used a barotropic equation of state and therefore does not treat thermal inertia effects (e.g. Goodwin, Whitworth & Ward-Thompson 2004; Bate 1998; Matsumoto & Hanawa 2003), although Banerjee et al. (2004) use the same procedure as used by us. By invoking line cooling we hope to capture more realistically the thermal behaviour of low-density material as it falls on to the protostellar disc and then spirals into the central protostar. Ideally radiative cooling should also take account of continuum emission from dust, in the regime where the gas and dust are thermally coupled, for example by using the algorithm developed by Stamatellos et al. (2007b). Unfortunately, the computational machinery to treat both the line cooling regime and the dust-cooling regime does not exist yet. We are currently developing an algorithm to combine these regimes.

3.3 Resolution

The core is modelled with 430 000 SPH particles, each having mass m_SPH= 1.4 × 10⁻⁵ M_⊙. Following Bate & Burkert (1997), the minimum resolvable mass is therefore ⁠. For comparison, the minimum possible Jeans mass (corresponding to ρ=ρ_CRIT and T= 9 K) is 1.5 × 10⁻³ M_⊙. Therefore, the Jeans mass is always resolved. In addition, there are 24 000 ambient particles exerting the external pressure that contains the core (see equation 3). These particles also have mass m_SPH= 1.4 × 10⁻⁵ M_⊙, but they have much higher temperature, and their density is adjusted so that they fill the space between the core and the periodic boundaries. We do not use sink particles. Instead, we limit the CPU time by imposing a minimum smoothing length, h_MIN= 2 au. Consequently, we are able to resolve the details of the protostellar disc, both radially and vertically, but not the central protostar.

3.4 Definitions

Any SPH particle whose density rises above ρ_★= 10⁻¹¹ g cm⁻³ is presumed to be part of a protostar. The primary protostar is the one that forms first, at the centre of the core, and t_O is its time of formation. Values of t_O are given in Table 1. t_O increases monotonically with Ω_O, because more rapidly rotating cores take longer to reach the critical density, ρ_★. We define the disc to be all material with density in the range 10⁻¹⁶ g cm⁻³≤ρ≤ 10⁻¹¹ g cm⁻³, since in general this material has |v_φ| ≫ |v_r| (where v_φ and v_r are the azimuthal and radial components of velocity).

The evolution of a core can usefully be divided into two distinct phases. During the Isothermal Collapse Phase, radiative cooling is very efficient and the gas temperature is roughly constant. Throughout most of the isothermal phase centrifugal acceleration is a secondary effect, and the core collapses on an approximately freefall time-scale. The Protostellar Phase starts as soon as a disc forms. The pattern of evolution during the Protostellar Phase depends critically on Ω_O. In cores with high Ω_O, (i) formation of the primary protostar is delayed until a substantial disc has built up, and (ii) the outer disc is prone to fragment producing secondary protostars. Since the Protostellar Phase is shorter than the Isothermal Collapse Phase, we give times in the form t=t_O+Δt, so that Δt is time measured from the moment when the primary protostar forms.

Simulations are compared at time t₂₈, which is defined to be when 28 per cent of the core mass (i.e. ∼1.7 M_⊙) has density ρ > 10⁻¹⁶ g cm⁻³, i.e. it is either in the disc or in a protostar. 28 per cent is the typical efficiency of star formation observed in cores (e.g. Motte, Andre & Neri 1998; Alves, Lombardi & Lada 2007). In order to follow the simulations further, we would need to relax the imposition of adiabaticity at densities ρ > 10⁻¹³ g cm⁻³, so that the energy dissipation resulting from the redistribution of angular momentum can be radiated away. Without this radiation the long-term viscous evolution of the disc is inhibited.

4 A LOW ANGULAR MOMENTUM CORE

In this section, we describe the evolution of a low angular momentum core, having Ω_O= 0.6 × 10⁻¹³ s⁻¹, and hence j_O= 1.1 × 10²¹ cm² s⁻¹ and β= 0.014. This is Run 1.

4.1 The Isothermal Collapse Phase

Fig. 1 shows the density profile at various times during the Isothermal Collapse Phase of Run 1. The collapse is roughly spherically symmetric (at least until the late stages) and approximately self-similar: the central region with a flat density profile shrinks, both in radial extent and in mass. The infalling envelope with ρ∝ |r|⁻² extends ever further inwards. Density profiles are computed by arranging the particles in order of increasing |r|, then averaging the density of contiguous sets of 100 particles; evidently this works only as long as departures from spherical symmetry are small.

Figure 1

Radial-density profiles during the Isothermal Phase of Run 1. Reading from the bottom, profiles are shown at t= 0, 41, 51, 61, 71, 81 and 91 kyr. The axes are logarithmic, and the dotted line gives a slope of −2.

Open in new tabDownload slide

4.2 The Protostellar Phase

Fig. 2 shows velocity vectors superimposed on false-colour images of the density field on 400 × 400 au² slices through the centre of the core during Run 1. These images correspond to times from just before the formation of the primary protostar at t=t_O to t=t₂₈. They show the formation of the primary protostar and its attendant protostellar disc. The density field is computed on a 64 × 64 grid by first evaluating the column density, Σ, through a slice having thickness Δ= 10 au, and then setting ρ=Σ/Δ. The velocity field is computed in an analogous way, but on a coarser 20 × 20 grid.

Figure 2

The density on 400 × 400 au² slices through the centre of the core during Run 1. The left- and right-hand columns show slices at y= 0, respectively. The times shown are (a) t_O− 1.9 kyr, (b) t_O+ 6.2 kyr and (c) t_O+ 18.2 kyr. The density scale ranges from 10⁻¹⁸ to 10⁻¹³ g cm⁻³, and contours are plotted at 10⁻¹⁶, 3 × 10⁻¹⁶ and 10⁻¹⁵ g cm⁻³. Velocity vectors are superimposed, and the insets show a 1 km s⁻¹ vector.

Open in new tabDownload slide

4.3 Density

Fig. 3(a) shows radial-density profiles on the equatorial plane, ρ(r, z= 0), from Run 1. Most of the mass interior to 10 au is in the primary protostar. Outside this the disc extends to beyond 100 au, with ρ∝r^−2.5. After t=t_O+ 8.1 kyr, the edge of the disc is very sharp.

Figure 3

Profiles from Run 1 at times t_O+ 0.2 kyr (i.e. just after the primary protostar appears), t_O+ 4.2 kyr, t_O+ 8.1 kyr (i.e. the stage at which the disc becomes well established), t_O+ 12.1 kyr, t_O+ 16.0 kyr and t_O+ 19.9 kyr (the end of the simulation). (a) The radial-density profile on the equatorial plane, ρ(r, z= 0). The axes are logarithmic, and a line with slope − 2.5 is shown to mimic the run of density in the disc between 10 and 200 au. The key gives the line-styles used to represent the different times. (b) The radial- temperature profile on the equatorial plane, T(r, z= 0). The axes are logarithmic, and a line with slope − 0.6 is shown to mimic the run of temperature in the disc between 10 and 200 au. (c) The radial-velocity profile on the equatorial plane, v_r(r, z= 0). The velocity axis is linear, and velocity is given in units of the local sound speed. The radius axis is logarithmic. (d) The azimuthal velocity profile on the equatorial plane, v_φ(r, z= 0). The axes are logarithmic, and lines with slopes of −1/2 and −1/4 are included to show that in the disc v_φ∝r^−1/4, corresponding to non-Keplerian differential rotation. Note that this plot has a larger radial extent than the others. (e) The radial profile of the disc aspect ratio, ⁠. Both axes are linear.

Open in new tabDownload slide

4.4 Temperature

Fig. 3(b) shows radial-temperature profiles on the equatorial plane, T(r, z= 0), from Run 1. Throughout most of the disc, the run of temperature on the mid-plane approximates to T∝r^−0.6. This is an intermediate between the scaling expected for a thin, non-self-gravitating, steady-state α-disc [which is heated by viscous dissipation, has a negative vertical temperature gradient, dT/d|z| < 0, and T_EFF∝r^−3/4 (Shakura & Syunyaev 1973)] and the scaling inferred from observations of protoplanetary discs around classical T Tauri stars [which are mainly heated by stellar irradiation, have a positive vertical temperature gradient, dT/d|z| > 0 and typically T_EFF∝r^−1/2 (Dullemond, Dominik & Natta 2001)]. Our simulations pertain to an earlier stage when the disc is self-gravitating and should be described as protostellar rather than protoplanetary. At this stage the infalling gas is heated impulsively when it encounters the accretion shock at the boundary of the disc; evidence for shock heating can be seen in the temperature spikes at the outer edge of the disc, in Fig. 3(b). Since the rate of infall on to the disc exceeds the rate of accretion through the disc and on to the primary protostar, the heating at the accretion shock dominates the energy budget of the disc. Subsequently, the gas is heated further by compression and viscous dissipation as it spirals inwards and on to the primary protostar, and therefore the vertical temperature gradient is negative, dT/d|z| < 0. Consequently, we anticipate that, as a disc adjusts from being protostellar to being protoplanetary, the 10 μm silicate feature will switch from absorption to emission, as the sign of the vertical temperature gradient changes from negative to positive.

4.5 Radial velocity

Fig. 3(c) shows radial-velocity profiles on the equatorial plane, normalized to the local isothermal sound speed, v_r(r, z= 0)/a(r, z= 0), from Run 1. Once the disc is well established (i.e. by t=t_O+ 8.1 kyr), the gas in the equatorial plane is in freefall, until it hits the accretion shock at the edge of the disc (100–200 au) at ∼ Mach 4. The gas away from the equatorial plane hits the accretion shocks bounding the disc above and below even faster, because these shocks are even deeper in the gravitational potential well of the core. Thus, the gas entering the disc is typically shock-heated to at least 100 K. The radial velocity of the gas in the disc is very tiny because the effective viscosity arising from gravitational torques is small, and therefore the gas takes many orbits to migrate inwards.

4.6 Azimuthal velocity

Fig. 3(d) shows azimuthal velocity profiles on the equatorial plane, v_φ(r, z= 0), from Run 1. Throughout the disc, the gas is in approximate centrifugal balance; because the disc is massive and self-gravitating, this results in a profile with v_φ(r, z= 0) ∝r^−1/4[significantly flatter than a purely Keplerian ⁠. Outside the disc, the infalling matter changes from its initial solid-body v_φ(r, z= 0) ∝r to v_φ(r, z= 0) ∝r⁻¹, as it falls inside ∼6000 au.

4.7 Scaleheight

At each radius r in the disc, we can define the vertical scaleheight as

11

where Σ(r) is the azimuthally averaged surface density (as seen from z=±∞) and ρ(r, z= 0) is the azimuthally averaged volume-density on the mid-plane(z= 0). Fig. 3(e) shows radial profiles of the disc aspect ratio, ⁠, from Run 1. The aspect ratio is almost constant, at a value ⁠. Thus, the disc is neither flat nor flaring; rather the iso-density surfaces are approximately conical. Since the mid-plane density varies as ρ(r, z= 0) ∝r^−5/2 and the scaleheight varies as ⁠, the surface density of the disc satisfies Σ(r) ∝r^−3/2. [This is very similar to what is inferred for protoplanetary discs by Dullemond et al. (2001), but our discs are much younger and more massive.] It follows that the mass of the disc out to radius r is roughly M_DISC(r) ∝r^1/2, and hence the equilibrium azimuthal velocity should approximate to

12

as indeed it does (see Fig. 3d).

4.8 Disc growth

Fig. 4 shows the growth of the primary protostar, M_★(t), and the growth of the disc, M_D(t). The growth rates are roughly constant at and ⁠; hence the disc grows faster than the primary protostar. The disc is marginally Toomre unstable (Toomre 1964, see Fig. 5), to the extent that weak spiral arms develop and contribute to the inward transport of material through the disc and on to the primary protostar. However, the disc is not self-regulating, in the sense that the rate of transport through the disc does not adjust so that ⁠. The disc is quite concentrated: typically more than 60 per cent of the disc mass has density ρ > 10⁻¹³ g cm⁻³. Radiative cooling times in the disc are sufficiently long so that the disc undergoes gravo-acoustic oscillations in the vertical direction, in response to the continuing infall of matter.

Figure 4

The masses of the primary protostar, M_★(t) (magenta), and the protostellar disc, M_D(t) (black), as a function of time, for Runs 1 and 7.

Open in new tabDownload slide

Figure 5

The Toomre parameter, Q_T, calculated for the face-on projections of Run 1 at t₀+ 18.3kyr (left-hand plot), and Run 7 at t₀+ 21.5 kyr (right-hand plot).

Open in new tabDownload slide

5 A HIGH ANGULAR MOMENTUM CORE

In this section, we describe the evolution of a high angular momentum core, having Ω_O= 1.8 × 10⁻¹³ s⁻¹, and hence j_O= 3.3 × 10²¹ cm² s⁻¹ and β= 0.13. This is Run 7.

5.1 The Isothermal Collapse phase

The Isothermal Collapse Phase of Run 7 proceeds in the same way as for Run 1, except that, due to the higher angular momentum, departures from spherical symmetry occur sooner and are more extensive. As a result, the Isothermal Collapse Phase ends with the formation of a protostellar disc, and the primary protostar does not form for another 15 kyr, by which stage the disc extends to 700 au.

5.2 The Protostellar Phase

Fig. 6 shows velocity vectors superimposed on false-colour images of the density field on 2000 × 2000 au² slices through the centre of the core during Run 7. These images correspond to times from just before the formation of the primary protostar at t=t_O to t=t₂₈. They show the growth and fragmentation of the protostellar disc. Densities and velocities are computed as in Fig. 2. The features that distinguish this high-Ω_O run (Run 7) from the low-Ω_O one discussed in the preceding section (Run 1) are (i) that the primary protostar grows much more slowly and acquires most of its mass via the disc (rather than by direct infall), and (ii) that the outer parts of the disc are able to fragment creating secondary protostars.

Figure 6

The density on slices through the centre of the core during Run 7. The left-hand (respectively right-hand) column shows slices at z= 0 (respectively y= 0). The times shown are (a) t_O+ 11.8 kyr; (b) t_O+ 29.7 kyr; (c) t_O+ 33.7 kyr; (d) t_O+ 33.7 kyr, zoomed to show the first long-lived fragment; (e) t_O+ 38.0 kyr and (f) t_O+ 40.0 kyr. In all cases the slices are 2000 × 2000 au². The density scale ranges from 10⁻¹⁷ to 10⁻¹² g cm⁻³, and contours are plotted at 10⁻¹⁶, 3 × 10⁻¹⁶ and 10⁻¹⁵ g cm⁻³. Velocity vectors are superimposed, and the insets show a 1 km s⁻¹ vector.

Open in new tabDownload slide

5.3 Density

Fig. 7(a) shows azimuthally averaged radial-density profiles on the equatorial plane, ρ(r, z= 0) from Run 7. Most of the mass interior to 10 au is in the primary protostar. Outside this the disc extends to ∼100 au with ρ∝r^−1.5, and then from ∼100 to ∼1000 au with ρ∝r^−2.0. The disc in Run 7 is far more extended than that in Run 1, and its density profile is significantly shallower. The spikes superimposed on the density profiles are due to spiral arms and fragments condensing out of the disc.

Figure 7

Profiles from Run 7 at times t_O (i.e. the moment at which the primary protostar forms), t_O+ 5.9 kyr, t_O+ 12.2 kyr, t_O+ 16.1 kyr, t_O+ 20.1 kyr, t_O+ 24.0 kyr, t_O+ 28.0 kyr and t_O+ 31.9 kyr (the end of the simulation). (a) The radial-density profile on the equatorial plane, ρ(r, z= 0). The axes are logarithmic, and lines with slopes −1.5 and −2.0 are shown to mimic the run of density in the disc from 10 to 100 au, and from 100 to 1000 au, respectively. The key gives the line-styles used to represent the different times. (b) The radial-temperature profile on the equatorial plane, T(r, z= 0). The axes are logarithmic, and a line with slope − 0.5 is shown to mimic the run of temperature in the disc between 10 and 1000 au. (c) The radial-velocity profile on the equatorial plane, v_r(r, z= 0). The velocity axis is linear, and velocity is given in units of the local sound speed. The radius axis is logarithmic. (d) The azimuthal velocity profile on the equatorial plane, v_φ(r, z= 0). The axes are logarithmic, and lines with slopes of −1/2 and −1/4 are included to show that in the disc v_φ∝r^−1/4, corresponding to non-Keplerian differential rotation. (e) The radial profile of the disc aspect ratio, ⁠. The axes are linear.

Open in new tabDownload slide

5.4 Temperature

Fig. 7(b) shows radial-temperature profiles on the equatorial plane, T(r, z= 0), from Run 7. Throughout most of the disc, the run of temperature on the mid-plane approximates to T∝r^−0.5, which is very similar to Run 1. However, in Run 7 the temperatures are lower. There are two reasons for this. First, because the disc is more extended, the material hitting the accretion shock at the boundary of the disc has fallen less far into the gravitational potential well of the core and is therefore moving more slowly; it is therefore not heated so much in the accretion shock. Secondly, again because the disc is more extended and therefore has lower surface density, but also because it is cooler, its cooling radiation can escape more readily. This in turn is the main reason why the disc in Run 7 is more prone to fragmentation. Since the rate of infall on to the disc greatly exceeds the rate of accretion on to the primary protostar (see below), the heating at the accretion shock again dominates the energy budget of the disc. Further heating by compression and viscous dissipation ensures that the vertical temperature gradient is negative, ⁠. However, because the outer disc is very cool, we only expect to see the 10 μm feature in absorption from the inner disc. The spikes superimposed on the temperature profiles are due to fragments condensing out of the disc and being heated by compression.

5.5 Radial velocity

Fig. 7(c) shows radial-velocity profiles on the equatorial plane, normalized to the local isothermal sound speed, v_r(r, z= 0)/a(r, z= 0), from Run 7. The gas in the equatorial plane is close to freefall until it hits the accretion shock at the outer edge of the disc. Because the disc is so extended, the material hitting the accretion shock is not moving very fast; for example, at the outer edge of the disc it is moving with v_r < Mach 2. Thus, the gas entering the disc is only heated mildly. The radial velocity of the gas in the disc is small, although there are significant perturbations in the vicinities of fragments.

5.6 Azimuthal velocity

Fig. 7(d) shows azimuthal velocity profiles on the equatorial plane, v_φ(r, z= 0), from Run 7. Throughout the disc the gas is in approximate centrifugal balance, but because the gravitational field is dominated by the mass of the disc, rather than the mass of the primary protostar, the profile is significantly flatter than Keplerian, i.e. 0 ≲ (− d ln v_φ/d ln r) ≲ 0.25, except for where it is perturbed by a forming fragment. Outside the disc, the infalling matter changes from its initial solid-body v_φ(r, z= 0) ∝r to v_φ(r, z= 0) ∝r⁻¹, as it falls inside ∼10 000 au.

5.7 Scaleheight

Fig. 7(e) shows radial profiles of the disc aspect ratio, ⁠, from Run 7. The aspect ratio is almost constant at a value within the inner disc (r < 200 au), and it falls off slightly at larger radii, where it is strongly perturbed by spiral arms and fragments. Again, it is neither flat nor flaring; the iso-density surfaces are approximately conical. Since the mid-plane density follows 1.5 ≲ (−d ln ρ/d ln r) ≲ 2.0 and the scaleheight varies as ⁠, the surface density of the disc follows 0.5 ≲ (−d ln Σ/d ln r) ≲ 1.0, and the disc mass follows 0.5 ≲ (d ln M_DISC/d ln r) ≲ 1.0. The azimuthal velocity should approximate to ⁠, hence 0 ≲ (− d ln v_φ/d ln r) ≲ 0.25, as indeed it does (see Fig. 7d), modulo the fluctuations due to forming fragments.

5.8 Disc growth

Fig. 4 shows the growth of the primary protostar, M_★(t), and the growth of the disc, M_D(t), in Run 7. Until the disc starts to fragment, and then a fragment merges with the primary protostar, the growth rates are roughly constant at and ⁠. Thus, the disc not only forms before the primary protostar, but also grows much faster.

5.9 Fragmentation

Because the disc is massive, extended and cool, it becomes Toomre unstable (Q_T≲ 1; see Fig. 5), and strong spiral arms develop. Where the cooling time in the disc is shorter than the dynamical times, the disc fragments (Gammie 2001; Rice et al. 2003). The disc is only weakly concentrated: typically less than 20 per cent of the disc mass has density ρ > 10⁻¹³ g cm⁻³, and most of this dense gas is located in the spiral arms and fragments.

5.10 Individual fragments

The first fragment condenses out at t∼t_O+ 30 kyr at ∼200 au. This fragment can be seen in Figs 6(b) and 8. At its inception, this fragment has mass ∼0.05 M_⊙. Subsequently, it grows in mass, but at the same time it spirals inwards and at t∼t_O+ 38 kyr it merges with the primary protostar. Two further fragments condense out at t∼t_O+ 40 kyr, at ∼400 au and ∼600 au (see the fourth and fifth panel of Fig. 6). The simulation is terminated soon after their formation, so it is unclear whether they also merge with the primary protostar, or whether they survive as secondary companions.

Figure 8

The density field on 140 × 140 au² slices through the centre of the first fragment to form during Run 7. The slice thickness is 6 au and the time shown is t_O+ 33.7 kyr. The density scale ranges from 10⁻¹⁷ to 10⁻¹² g cm⁻³, and contours are plotted at 10⁻¹⁶, 3 × 10⁻¹⁶ and 10⁻¹⁵ g cm⁻³. Velocity vectors are superimposed, and the insets show a 1 km s⁻¹ vector.

Open in new tabDownload slide

5.11 Circum-fragmentary discs

Fig. 8 shows velocity vectors superimposed on false-colour images of the density field in 140 × 140 au² slices through the first fragment to form. Densities and velocities are computed as in Fig. 2, but using co-ordinates measured relative to the centre of mass of the fragment. We see that the fragment is surrounded by a well-defined, rotationally supported circum-fragmentary disc (CFD). Very similar discs are found around the other fragments. Typically they have radii between ∼40 and ∼60 au.

5.12 Transport of angular momentum

Although most of the core mass infalls directly on to the disc, the spiral arms in the disc are effective at transporting angular momentum outwards, by means of gravitational torques, thereby allowing material to spiral inwards and accrete on to the primary protostar. Fragmentation further enhances the transport of angular momentum. This is because CFDs rotate in the same sense as the mother disc. Consequently, on the inside of a CFD (i.e. the side closest to the primary protostar) the shear between the material in the CFD and the material in the mother disc acts to slow down the material in the mother disc, thereby reducing its angular momentum. Conversely, on the outside of the CFD (i.e. further from the primary protostar), the shear between the material in the CFD and the material in the mother disc acts to speed up the material in the mother disc, thereby increasing its angular momentum.

5.13 Merging

Efficient transport of angular momentum by the CFD may cause a fragment to spiral inwards and merge with the primary protostar. However, it is possible that the inspiral is being artificially accelerated by our imposition of adiabaticity at densities ρ > 10⁻¹³ g cm⁻³, since this forces a fragment and its CFD to remain hot, extended and therefore very dissipative, rather than allowing them to continue contracting towards protostellar densities.

5.14 Disc mineralogy

Whereas in Run 7 the azimuthally averaged disc temperature never rises above ∼200 K, in Run 1 it rises above ∼500 K, and locally it reaches ∼800 K, sufficiently high to produce crystalline silicates. Thus, we might expect crystalline silicates to be present in discs forming from low-Ω_O cores, and hence to give rise to characteristic features in the SEDs of single protostars. Conversely, we might expect crystalline silicates to be rare in the cold, extended discs forming from high-Ω_O cores, and hence their characteristic features should be weak in – or even absent from – the SEDs of multiple protostars. The existence of a link between the amount of crystalline dust which can be found in protoplanetary discs and the initial conditions within the parental molecular cloud core has already been suggested by Dullemond, Apai & Walch (2006).

6 CORES WITH INTERMEDIATE ANGULAR MOMENTA

Referring now to the complete ensemble of Runs (1–7) (see Fig. 9), we note that, as Ω_O increases, there is an essentially monotonic shift in the outcome; small departures from monotonicity at high values of Ω_O can be attributed to the stochastic nature of fragmentation. Key parameters are presented in Table 1.

Figure 9

Disc formation as a function of core angular momentum. For all seven runs, we show the density of the discs within slices of thickness 10 au about the disc mid-plane (z= 0) at t_final on the same scale. From top (Run 1) to bottom (Run 7) the core rotation Ω₀ and thus the total initial core angular momentum are increasing.

Open in new tabDownload slide

For low Ω_O, a substantial primary protostar is formed, surrounded by a compact, thick, hot disc; the disc is unable to fragment and simply accretes on to the primary protostar. As Ω_O is increased, the growth rate of the primary protostar decreases, not least because it acquires an increasing fraction of its mass via the disc (rather than by direct infall). At the same time the disc becomes more massive, more extended and cooler. Eventually, for Ω_O > 1.35 × 10⁻¹³ s⁻¹ (Runs 5, 6 and 7), the disc becomes unstable against fragmentation, i.e. it satisfies both the Toomre and Gammie conditions (Toomre 1964; Gammie 2001). We are unable to ascertain whether all the resulting fragments merge with the primary protostar or whether some of them survive as secondary companions.

As Ω_O increases, the time, t_O, at which the primary protostar forms, increases, and its subsequent growth rate, ⁠, decreases. Consequently, the overall pattern of accretion shifts from one in which the primary protostar forms first, and then the protostellar disc (low Ω_O), to one in which the disc forms before the primary protostar (high Ω_O).

Fig. 5 shows the Toomre parameter, Q_T, calculated for the face-on projections of Run 1 at t₀+ 18.3 kyr (the left-hand plot), and Run 7 at t₀+ 21.5kyr (the right-hand plot). In the warm, compact disc from Run 1, the lowest values of Q_T are in the range 1 ≲Q_T≲ 3, so that the disc develops weak spiral arms that transport angular momentum by exerting gravitational torques, but it does not fragment (cf. Lin & Pringle 1990; Laughlin & Bodenheimer 1994; Lodato & Rice 2004; Mejía et al. 2005). In the cool, extended disc from Run 7, Q_T falls below 1, leading to the growth of massive spiral arms and then fragmentation. This is consistent with recent findings of Rafikov (2005), Matzner & Levin (2005), Whitworth & Stamatellos (2006), Stamatellos, Hubber & Whitworth (2007a), Kratter, Matzner & Krumholz (2008), Rice & Armitage (2009), Clarke (2009), who suggest that in the inner disc (r≲ 50 au) fragmentation is highly unlikely, whereas it is possible at larger radii.

7 COMPARISON WITH OBSERVATIONS

In the Appendix, we develop an analytic model for the early growth and fragmentation of a protostellar disc, in terms of its mass, M_O, its initial radius, R_O, and its initial angular rotation speed, Ω_O. We show that cores with Q_F < 1, where

13

should be prone to disc fragmentation, whilst material is still infalling from the core envelope. We can apply this criterion to the samples of cores observed by Goodman et al. (1993), Barranco & Goodman (1998) and Caselli et al. (2002), for which values of M_O, R_O and Ω_O are listed. The resulting values of Q_F are plotted against M_O in Fig. 10.

Figure 10

The parameter Q_F (see equation 13) is plotted against core mass, M_O, for the core samples of Goodman et al. (1993), Barranco & Goodman (1998) and Caselli et al. (2002). Cores with Q_F < 1 are predicted to undergo disc fragmentation.

Open in new tabDownload slide

Although a few cores have Q_F∼ 1, none has Q_F < 1. We infer that none of these cores is strongly pre-disposed to fragment whilst material is still infalling from the core envelope. The observed angular speeds are too small, and the strong heating at the accretion shock (where the infalling envelope material meets the disc) acts to stabilize the disc.

This inference must be qualified as follows. (i) The core parameters used are very uncertain, in particular the core masses. (ii) Many of the cores appear already to have formed protostars (as evidenced by their association with IRAS sources), and even those that have not may already be in an advanced state of collapse. Therefore, the observed values of R_O and Ω_O may not be ‘initial’ in the sense we have defined them. (iii) Our simulations do not take account of the low levels of turbulence observed in such cores. (iv) We do not rule out the possibility that at a later stage, once infall ceases, some of these protostellar discs cool and fragment, particularly if more efficient cooling due to dust operates.

8 CONCLUSIONS

We have performed simulations of the collapse of rigidly rotating pre-stellar cores. Each core is initially a mildly supercritical BES (ξ_B= 6.9), with its density increased by 10 per cent, and contained by external pressure. The gas is pure molecular hydrogen with T= 28 K (hence the isothermal sound speed is 0.34 × 10⁵ cm s⁻¹), the initial central density is ρ_C= 10⁻¹⁸ g cm⁻³, the core mass is M_O= 6.1 M_⊙, its initial radius is R_O= 17 000 au and its external pressure is P_EXT=k_B 5.5 × 10⁵ cm⁻³ K. The initial angular speed is varied from Ω_O= 0.6 × 10⁻¹³ s⁻¹ to Ω_O= 1.8 × 10⁻¹³ s⁻¹ (corresponding to 0.013 < β < 0.12).

A slowly rotating core collapses to form a relatively massive primary protostar, and then a protostellar disc grows around this protostar. The disc is massive, but quite hot and compact, and so it does not fragment; it simply accretes on to the primary protostar. There is some observational evidence for the existence of massive discs (e.g. Eisner et al. 2005; Rodríguez et al. 2005; Greaves et al. 2008; Mann & Williams 2009). However, the phase we have simulated is short lived and still deeply embedded in the parental core, so our discs would be rather hard to observe. We conjecture that such discs will show signatures of crystalline silicates out to disc radii of several au.

A rapidly rotating core collapses initially to form a protostellar disc, and then a protostar condenses out at the centre. The protostar acquires almost all its mass via the disc, and therefore it grows quite slowly. The disc is extended, cool and thin, and tends to fragment.

On the basis of an analytic model, we obtain a condition on the mass, initial radius and initial angular speed of a core (equation 13) such that its protostellar disc will fragment, whilst material is still infalling from the core envelope. Applying this condition to observed samples of cores, we predict that such fragmentation is unlikely.

In a companion paper (Walch et al., in preparation), we explore how these results are fundamentally changed if a pre-stellar core is turbulent, but has the same net angular momentum as the cores modelled here.

We thank D. Neufeld for providing the cooling tables. SW performed this work with support from the International Max-Planck Research School, the DFG Cluster of Excellence Origin and Structure of the Universe (www.universe-cluster.de) and the Marie Curie RTN Constellation (MRTN-CT-2006-035890). We thank the referee for his/her detailed and thoughtful comments, which were very helpful.

REFERENCES

Appendix

APPENDIX A: Analytic model

We develop here a simple model to explain why rapidly rotating cores tend to fragment (and thereby presumably produce multiple systems), whereas slowly rotating cores do not. For simplicity, we assume that the initial core is cylindrical, with height equal to its diameter and hence initial density ρ_O=M_O/2πR³_O; we make this assumption purely because it allows us to formulate analytically the structure of the resulting disc. The core is initially in solid-body rotation with angular speed, Ω_O, and we assume that during the Isothermal Collapse Phase, and during the early Protostellar Disc Phase, angular momentum is conserved minutely. Consequently, the material initially in an annulus at radius r_O has mass ‘interior to it’ in the disc, and ends up at radius r with angular speed Ω_D(r) =Ω_Or²_O/r². Centrifugal balance then requires that ⁠, and hence

(A1)

The edge of the disc is therefore at

(A2)

and

(A3)

It follows that the final surface density of the disc is

(A4)

Infalling material impinges on the disc at radius r with speed

(A5)

We shall assume that it starts arriving at time

(A6)

and continues arriving at a constant rate until time 2t_D(r). Hence, in the time interval t_D(r) ≤t≤ 2t_D(r), the flux of material on to one side of the disc at radius r is

(A7)

the instantaneous surface density of the disc is

(A8)

the ram pressure of the infalling material is

(A9)

and the flux of energy carried by the infalling material, and dissipated in the accretion shock, is

(A10)

The disc cools mainly by molecular-line radiation. Under the conditions obtained in the discs simulated here, the lines tend to be thermalized and optically thick. Hence, the cooling rate per unit area of the disc, due to one species, M, is

(A11)

where I_M is the moment of inertia of M, A_M is the Einstein-A coefficient of the J= 1 → 0 transition of M, Z_M is the fractional abundance of M (by mass). Since we are neglecting redistribution of angular momentum, we also neglect viscous heating of the disc; this is a reasonable assumption, given that the disc grows faster than the primary protostar (i.e. the viscous evolution of the disc proceeds on a longer time-scale than infall on to the disc). Hence the principal heating mechanism is dissipation of kinetic energy in the accretion shock, and therefore we estimate the temperature by setting F_D(r, t) ∼F_M(r, t). It follows that, during the period when the disc is growing,

(A12)

The Toomre and Gammie parameters for the actively accreting part of the disc are

(A13)

(A14)

To proceed, we consider cooling by CO, for which I_CO≃ 1.4 × 10⁻³⁹ g cm², A_CO≃ 10⁻⁷ s⁻¹ and Z_CO≃ 0.002; and we introduce dimensionless core parameters,

Then,

(A15)

(A16)

(A17)

The mass of the system (primary protostar plus protostellar disc) is given by

(A18)

It follows that at t=t_FF (i.e. ⁠), M_SYS≃ 0.29M_O, and since this is very close both to the juncture when the disc is most unstable and to the juncture at which we compare simulations with different Ω_O, we will evaluate the stability of the discs at this time too. At t=t_FF, the disc is complete out to R_D/16, and is just starting to form at R_D, so between these radii the disc is actively accreting. From equation (A17), we see that the accreting disc is most unstable (Toomre Q_T∼ 1.07) at R_D/16, and this is basically where the simulated discs with high Ω_O tend to fragment. The condition for the disc to be able to cool fast enough to fragment at R_D/16 (i.e. Gammie Q_G < 1) then reduces to

(A19)

Equation (A19) agrees with the numerical results, which give fragmentation for but not for ⁠. It also indicates how the criterion for fragmentation scales with core mass and radius, and we exploit this to analyse observed cores in Section 7. The analytic model also accords well with the simulations in predicting that fragmentation should occur near R_D/16, which for those discs that fragment is between ∼200 and ∼700 au (cf. Matzner & Levin 2005; Rafikov 2005; Whitworth & Stamatellos 2006).