A&A 423, 169-182 (2004)
DOI: 10.1051/0004-6361:20040285

Simulating star formation in molecular cores

II. The effects of different levels of turbulence

S. P. Goodwin - A. P. Whitworth - D. Ward-Thompson

Dept. of Physics & Astronomy, Cardiff University, 5 The Parade, Cardiff CF24 3YB, UK

Received 17 February 2004 / Accepted 1 May 2004

Abstract
We explore, by means of a large ensemble of SPH simulations, how the level of turbulence affects the collapse and fragmentation of a star-forming core. All our simulated cores have the same mass ( $5.4 ~M_\odot$), the same initial density profile (chosen to fit observations of L1544), and the same barotropic equation of state, but we vary (a) the initial level of turbulence (as measured by the ratio of turbulent to gravitational energy, $\alpha_{\rm turb} \equiv U_{\rm turb}/\vert\Omega\vert = 0,~0.01,~0.025,~ 0.05,~0.10\;{\rm and}\;0.25$) and (b), for fixed $\alpha _{\rm turb}$, the details of the initial turbulent velocity field (so as to obtain good statistics).

A low level of turbulence ( $\alpha_{\rm turb} \sim 0.05$) suffices to produce multiple systems, and as $\alpha _{\rm turb}$ is increased, the number of objects formed and the companion frequency both increase. The mass function is bimodal, with a flat low-mass segment representing single objects ejected from the core before they can accrete much, and a Gaussian high-mass segment representing objects which because they remain in the core grow by accretion and tend to pair up in multiple systems.

The binary statistics reported for field G-dwarfs by Duquennoy & Mayor (1991, A&A, 248, 485) are only reproduced with $\alpha_{\rm turb} \sim 0.05$. For much lower values of $\alpha _{\rm turb}$ ($\la $0.025), insufficient binaries are formed. For higher values of $\alpha _{\rm turb}$ ($\ga $0.10), there is a significant sub-population of binaries with small semi-major axis and large mass-ratio (i.e. close binaries with components of comparable mass). This sub-population is not present in Duquennoy & Mayor's sample, although there is some evidence for it in the pre-Main Sequence population of Taurus analyzed by White & Ghez (2001, ApJ, 556, 265). It arises because with larger $\alpha _{\rm turb}$, more low-mass objects are formed, and so there is more scope for the binaries remaining in the core to be hardened by ejecting these low-mass objects. Hard binaries thus formed then tend to grow towards comparable mass by competitive accretion of material with relatively high specific angular momentum.

Key words: methods: numerical - stellar dynamics - stars: formation - ISM: general

1 Introduction

Turbulence appears to play a crucial role in the structure and evolution of molecular clouds, in the formation of star-forming cores within molecular clouds, and in the collapse and fragmentation of cores to form protostars.

The main evidence for turbulence in molecular clouds comes from their apparently fractal substructure (e.g., Elmegreen & Falgarone 1996; Elmegreen 2002), and from the almost universal power-law scaling relations between size (L) line-width ($\Delta v$) and mass (M)[*] and the almost universal power-law mass function[*] to which this substructure subscribes, over many orders of magnitude, from the largest molecular cloud complexes ( ${\sim}10^2~{\rm pc}$, ${\sim}10^6 ~M_\odot$), down to the smallest resolvable structures ( ${\sim}10^{-2}~{\rm pc}$, ${\sim}10^{-4} ~M_\odot$) (e.g., Larson 1981; Myers 1983; Stutzki & Güsten 1990; Hobson 1992; Hobson et al. 1994; Williams et al. 1994; Elmegreen & Falgarone 1996; Kramer et al. 1996; Kramer et al. 1998; Heithausen et al. 1998).

Until recently, it had been presumed that molecular clouds were long lived, being supported against collapse by their internal turbulence, and this was advanced as the reason for the low overall efficiency of star formation. However, it is now recognized that turbulence cannot support clouds for long, because - even with a frozen-in magnetic field - the turbulence dissipates on a dynamical timescale (Mac Low et al. 1998; Stone et al. 1998). Instead clouds are presumed to form and disperse on a dynamical timescale, without ever reaching equilibrium (Ballesteros-Paredes et al. 1999; Elmegreen 2000; Pringle et al. 2001; Hartmann et al. 2001).

In this highly dynamical scenario, cores form wherever a sufficiently dense and coherent converging flow is created by the turbulent velocity field (Elmegreen 1997; Padoan et al. 1997; Hartmann et al. 2001; Klessen & Burkert 2000, 2001; Klessen et al. 2000; Padoan & Nordlund 2002; Mac Low & Klessen 2004). Frequently these cores are not gravitationally bound, and therefore they disperse soon after they form. However, occasionally they are gravitationally bound, and in this case they are likely to proceed straight into gravitational collapse; these are the cores we identify as "prestellar''.

This scenario is in contrast to the idea that prestellar cores are supported magnetically, and evolve quasistatically by ambipolar diffusion, until they become magnetically supercritical and collapse (e.g., Basu & Mouschovias 1994, 1995a,b; Ciolek & Mouschovias 1993, 1994, 1995; Morton et al. 1994; Ciolek & Basu 2000). The main effects of the quasistatic ambipolar diffusion phase are (i) to give the core more time to lose angular momentum by magnetic braking; (ii) to organize the material so that its subsequent collapse is rather well focussed; and (iii) to allow turbulence to decay. All three effects mean that such cores are less likely to form multiple systems. Since most stars are observed to be in multiple systems (e.g., Duquennoy & Mayor 1991; Fischer & Marcy 1992; White & Ghez 2001), and since there is no observational evidence for magnetically subcritical cores (e.g., Crutcher 1999; Bourke et al. 2001; Crutcher et al. 2003), we shall assume that ambipolar diffusion does not play an important role in the evolution of prestellar cores. For simplicity, we ignore the magnetic field altogether.

In the highly dynamic scenario the collapse of a prestellar core is far more likely to lead to fragmentation and the formation of multiple systems (e.g., Whitworth et al. 1995; Turner et al. 1995; Whitworth et al. 1996; Klein et al. 2001, 2003; Bate et al. 2002a,b, 2003; Bonnell et al. 2003; Delgado-Donate 2003, 2004; Goodwin et al. 2004a,b; Hennebelle et al. 2003, 2004). This is because in the highly dynamic scenario prestellar cores are formed non-quasistatically and therefore (a) they are launched directly into the non-linear regime of gravitational collapse, and (b) they are likely to have retained some internal turbulence.

Burkert & Bodenheimer (2000) have pointed out that the internal turbulence in molecular cores can be represented by a Gaussian random velocity field having a power spectrum of the form $P(k) \propto k^\eta$, with $\eta \sim -3$ to -4. This not only reproduces the observed scaling relation between size and linewidth. It also reproduces the observed rotation of molecular cores. Thus there is no need to invoke ordered rotation as an additional source of support in molecular cores, and indeed there is no observational evidence for significant ordered rotation in prestellar cores (e.g., Jessop & Ward-Thompson 2001).

We have therefore undertaken a numerical study of the influence of turbulence on the collapse and fragmentation of prestellar cores. We have taken as our reference point a simple model of the core L1544, and in the first paper of this series (Goodwin et al. 2004a, hereafter Paper I) we have shown that cores with even a low level of turbulent energy routinely spawn multiple stellar systems. Specifically, in an ensemble of 20 simulations of the collapse of $5.4 ~M_\odot$ cores having an initial ratio of turbulent to gravitational energy

 \begin{displaymath} \alpha_{\rm turb} \equiv \frac{U_{\rm turb}}{\vert\Omega\vert} = 0.05 , \end{displaymath} (1)

$80\%$ of the cores form at least two, and in one case ten, objects (stars and brown dwarfs). In addition, the distributions of semi-major axis, mass ratio and eccentricity for the resulting multiple systems are consistent with the distributions for observed binary systems (e.g., Duquennoy & Mayor 1991, hereafter DM91). Paper I has also shown that while low levels of turbulence can easily produce multiple fragmentation, the outcome of any one simulation depends sensitively on the exact details of the turbulent velocity field. Consequently a statistical approach is essential in evaluating the influence of turbulence on multiple star formation in cores.

In this paper we extend the simulations of Paper I to examine the effect of different levels of turbulence on star formation within dense molecular cores. Using exactly the same core structure as in Paper I we simulate ensembles of between 10 and 20 cores with $\alpha_{\rm turb} = 0,~0.01,~0.025,~0.05,~0.10\;{\rm and}\;0.25$. We examine the numbers and masses of stars and brown dwarfs that form and the properties of the multiple systems to which some of them belong.

We note that these levels of turbulence involve much lower non-thermal velocities than the earlier work of Whitworth et al. (1995), Turner et al. (1995), Whitworth et al. (1996), Klein et al. (2001, 2003), Bate et al. (2002a,b, 2003), Bonnell et al. (2003), Delgado-Donate (2003, 2004). Consequently they may be applicable to scenarios in which instability develops more quasistatically due to ambipolar diffusion; provided that some turbulence can persist through (or be regenerated after) the ambipolar diffusion phase, and provided the subsequent collapse is sufficiently rapid.

In Sect. 2 we define the initial conditions for the simulations. In Sect. 3 we describe the code and the constitutive physics used. In Sect. 4 we outline the different ensembles of simulations performed with different values of $\alpha _{\rm turb}$, and in Sect. 5 we collate the statistics from the different ensembles. Section 6 discusses the statistics in terms of the underlying physics, and Sect. 7 gives our main conclusions.

  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig1.ps}\end{figure} Figure 1: The filled circles give estimated values of $\alpha _{\rm turb}$ (the ratio of turbulent to gravitational energy, see Eq. (1)) and $M_{\rm core}$ (core mass) for the starless cores in the Jijina et al. (1999) catalogue. The hatched region gives the range of $(\alpha _{\rm turb},M_{\rm core})$ values used in this paper (i.e. $0 \leq \alpha _{\rm turb} \leq 0.25$ and $M_{\rm core} = 5.4 ~M_\odot$). The cross (+) and the star (*) give - respectively - the $(\alpha _{\rm turb},M_{\rm core})$ values used by Bate et al. (2002a,b, 2003) and Delgado Donate et al. (2003a,b).
Open with DEXTER

2 Initial conditions

Molecular cores which are associated with IRAS sources are presumed to have already formed a protostar, and are classified as protostellar cores, whereas those which have no associated IRAS source are classified as starless cores (Beichman et al. 1986). The densest starless cores are presumed to be on their way to forming stars, and are therefore classified as prestellar cores (Ward-Thompson et al. 1994). We base our initial conditions on the observed properties of prestellar cores.

The density in a pre-stellar core is approximately uniform in the inner few thousand au, but further out it decreases as $r^{-\eta}$ with $2 \la \eta \la 5$, and eventually it merges with the background (e.g., Ward-Thompson et al. 1994; André et al. 1996; Ward-Thompson et al. 1999; André et al. 2000; Tafalla et al. 2004). A good fit to the density in a pre-stellar core is given by a Plummer-like profile,

\begin{displaymath}\rho(r) = \frac{\rho_{\rm kernel}}{(1 + (r/R_{\rm kernel})^2)^2} , \end{displaymath} (2)

where $\rho_{\rm kernel}$ is the central density and $R_{\rm kernel}$ is the radius of the region in which the density is approximately uniform. We set $\rho_{\rm kernel} = 3 \times 10^{-18}~{\rm g}~{\rm cm}^{-3}$ and $R_{\rm kernel} = 5 000~{\rm au}$, with the outer boundary of the core at $R_{\rm core} = 50~000~{\rm au}$, so the core has total mass $M_{\rm core} = 5.4 ~M_\odot$. The core is initially isothermal with $T = 10~{\rm K}$, and uniformly molecular with mean gas-particle mass $\bar{m} = 4 \times 10^{-24}~{\rm g}$. This means that the core has a ratio of thermal to gravitational energy of

\begin{displaymath}\alpha_{\rm therm} \equiv \frac{U_{\rm therm}}{\vert\Omega\vert} \simeq 0.3 . \end{displaymath} (3)

Thus far, these are exactly the same initial conditions as we used in Paper I.

Table 1: For each value of $\alpha _{\rm turb}$, we list the number of realizations simulated ( ${\cal N}_{\rm real}$), the mean mass in objects at the end of the simulations ( $M_{\rm tot}/M_\odot $), the average number of objects formed per simulation $\left (\left <{\cal N}_{\rm obj}\right >\right )$, the companion probability (cp), the companion frequency (cf), the multiplicity frequency (mf) and the pairing factor (pf). The companion frequency, cf, is given for all objects and then seperately for low-mass objects ( $M < 0.5 ~M_\odot$) and for high-mass objects ( $M > 0.5 ~M_\odot$).

Molecular cores show a significant non-thermal contribution to their line-widths, which is attributable to internal turbulence. Figure 1 shows the estimated ratios of turbulent to gravitational energy, $\alpha _{\rm turb}$ (see Eq. (1)), and the estimated masses, $M_{\rm core}$, for prestellar cores from the Jijina et al. (1999) catalogue. These cores were selected as prestellar on the basis of their having low temperature ( ${<}20~{\rm K}$), no associated IRAS source and no observed outflow. The shaded area in Fig. 1 shows the region of parameter space that the simulations in this paper cover, i.e. a $5.4 ~M_\odot$ core with a range of $\alpha _{\rm turb}$ from 0 to 0.25.

To model the turbulence, a divergence-free Gaussian random velocity field is superimposed on the core (cf. Bate et al. 2002a,b; Bate et al. 2003; Fisher 2004; Bonnell et al. 2003; Delgado-Donate et al. 2003, 2004). The power spectrum of the velocity field is $P(k) \propto k^{-4}$, so as to mimic the observed scaling between size and line-width in interstellar gas clouds (Larson 1981; Burkert & Bodenheimer 2000). The magnitude of the velocity field is normalized so that $\alpha_{\rm turb} = 0$, $0.01,~0.025,~0.05,~0.1\;{\rm or}\; 0.25~$, and at least ten realizations have been made with each value of $\alpha _{\rm turb}$ (as summarised in Table 1).

In each realization the random number seed for the turbulence is different, and hence the detailed structure of the velocity field is different. It is essential that many different realizations be performed for a given value of $\alpha _{\rm turb}$, because the mix of protostars and multiple systems which forms depends critically on the details of the velocity field. Therefore different values of $\alpha _{\rm turb}$ can only be compared in a statistical sense by performing an ensemble of different realizations for each representative value of $\alpha _{\rm turb}$.

Each simulation is evolved for $0.30~{\rm Myr}$. We have chosen this end-time for three reasons. (i) It facilitates comparison with the $\alpha _{\rm turb} = 0.05$ simulations reported in Paper I; (ii) in most simulations, object formation (i.e. sink creation) ceases around $0.15~{\rm to}~0.20~{\rm Myr}$, and in only one case is a new object formed after $0.25~{\rm Myr}$; (iii) by this stage, ${\sim} 66\%$ of the mass is already in stars and brown dwarves (this fraction decreases slightly with increasing $\alpha _{\rm turb}$), and so the remaining diffuse gas is likely to be affected by negative feedback from the existing stars and brown dwarves; feedback is not included in these simulations, although we are currently exploring its effect (Boyd et al., in preparation).

The properties of the resulting protostars in each ensemble are compared, both as a function of $\alpha _{\rm turb}$, and against the observed properties of young stellar objects. In this latter regard, the observational samples used for comparison are the local field G-dwarfs studied by DM91 and the multiple systems in Taurus studied by White & Ghez (2001, hereafter WG01).

3 Computational method and constitutive physics

The simulations are performed using SPH (Lucy 1977; Gingold & Monaghan 1977; Monaghan 1992). Our SPH code ( DRAGON) uses the standard M 4 kernel (Monaghan & Lattanzio 1985) and invokes variable smoothing lengths (so that each particle has ${\cal N}_{\rm neib} = 50 \pm 5$ neighbours). An octal tree is built to facilitate the computation of gravitational accelerations and the identification of neighbours. Gravity is kernel-softened with the particle smoothing length, and standard artificial viscosity (Monaghan 1992) is included, with $\alpha_v=1$ and $\beta_v=2$.

We use a barotropic equation of state,

 \begin{displaymath} \frac{P(\rho)}{\rho} \equiv c_{\rm s}^2(\rho) = c_0^2 \left[... ...eft ( \frac{\rho}{\rho_{\rm crit}} \right)^{2/3} \right] \cdot \end{displaymath} (4)

Thus at low densities ( $\rho < \rho_{\rm crit} = 10^{-13}~{\rm g}~{\rm cm}^{-3}$, where radiative cooling is optically thin) the gas is approximately isothermal with $c_{\rm s} \simeq c_0 = 0.2~{\rm km}~{\rm s}^{-1}$, corresponding to molecular gas at temperature $T = 10~{\rm K}$; and at high densities, ( $\rho > \rho_{\rm crit}$, where cooling radiation is trapped by the high optical-depth) the gas behaves adiabatically, with ratio of specific heats $\gamma \simeq 5/3$. ( $\gamma \simeq 5/3$ is the appropriate ratio of specific heats because the rotational degrees of freedom of H2 are not excited until the temperature rises above ${\sim} 300~{\rm K}$, and our simulations do not follow the gas dynamics to such high temperatures; see below). Equation (4) reproduces the thermal behaviour of the gas in detailed one-dimensional simulations of the collapse of $1 ~M_\odot$ cores (e.g., Larson 1969; Tohline 1982; Masunaga & Inutsuka 2000).

Whenever a gravitationally bound condensation forms and the density of an SPH particle within the condensation rises above $100 \rho_{\rm crit} = 10^{-11}~{\rm g}~{\rm cm}^{-3}$, all the particles within $5~{\rm au}$ of that particle are replaced with a sink particle having radius $R_{\rm sink} = 5~{\rm au}$. (From Eq. (4), we can estimate that the temperature of gas at $\rho = 100 \rho_{\rm crit}$ is ${\sim} 225~{\rm K}$, and so the introduction of sink particles makes it unnecessary to treat the thermal behavious of the gas at temperatures above ${\sim} 225~{\rm K}$.) Sink particles interact with the gas both gravitationally, and by accreting SPH particles that enter the sink radius and are bound to the sink (see Bate et al. 1995, for a detailed description of sink particles). As in Paper I, we refer to sink particles as "objects'', and then more specifically to "stars'' when the sink mass is ${>} 0.08 ~M_{\odot}$ and "brown dwarfs'' when the mass is lower than this.

4 Ensembles of simulations with different $\alpha _{\rm turb}$

4.1 ${\alpha _{turb} = 0}$ (no turbulence)

As reported in Paper I, when a core has no turbulence only one object forms, very close to the centre of the core. The evolution follows very closely the semi-analytic model of Whitworth & Ward-Thompson (2001). In particular, the accretion rate is very large early on, and then decreases. After $0.3~{\rm Myr}$ the stellar mass reaches $3.75 ~M_{\odot}$. (In this case the ten different realizations involve different initial SPH particle positions.)

4.2 ${\alpha _{turb} = 0.01}$

The ensemble of ten simulations with $\alpha_{\rm turb} = 0.01$ is virtually identical to that with $\alpha_{\rm turb} = 0.00$ (no turbulence). Only one star ever forms, and this happens about one free-fall time ( ${\sim} 0.055~{\rm Myr}$) after the start of the simulation, always close to the centre of mass of the core. This level of turbulence is apparently too low to induce multiple fragmentation.

4.3 ${\alpha _{turb} = 0.025}$

$\alpha_{\rm turb} = 0.025$ appears to be approximately the minimum level of turbulence required for multiple object formation. Of the ten simulations with $\alpha_{\rm turb} = 0.025$, eight produce only a single star (as with lower levels of turbulence), but one simulation produces six objects, and one produces nine objects.

Specifically, this last simulation produces three intermediate-mass stars in an hierarchical triple system embedded in the core (a tight binary with component masses $1.07 ~M_\odot$ and $0.88 ~M_\odot$ and semi-major axis $8.6~{\rm au}$, plus a third star with mass $1.65 ~M_\odot$ orbiting at ${\sim} 70~{\rm au}$). In addition, the simulation produces one very low-mass star ( $0.087 ~M_\odot$) and five brown dwarfs ( $0.034 ~M_\odot$ to $0.072 ~M_\odot$), all of which are ejected from the core.


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig2.ps}\end{figure} Figure 2: The average number of objects formed in a core, $\left <{\cal N}_{\rm obj}\right >$, as a function of the ratio of turbulent to gravitational energy $\alpha _{\rm turb}$. The error bars give the standard deviation (see Eq. (5)).
Open with DEXTER

4.4 ${\alpha _{turb} = 0.05}$

Of the twenty simulations performed with $\alpha _{\rm turb} = 0.05$, four produce just a single star, and the remaining sixteen produce 71 stars and 16 brown dwarfs between them (between 2 and 10 objects per simulation). Of these 71 stars, 44 remain in the core in multiple systems, and the rest are ejected from the core. Of the 16 brown dwarfs, 15 are ejected, and only one remains in a binary system in the core. The mean number of objects formed per simulation is 4.55. Further details of this ensemble of simulations are given in Paper I.

4.5 ${\alpha _{turb} = 0.10}$

Of the twenty simulations performed with $\alpha _{\rm turb}=0.10$, five produce just a single star, and the remaining fifteen produce 76 stars and 19 brown dwarfs between them (between 3 and 10 objects per simulation). Of these 76 stars, 53 remain in the core in multiple systems, and the rest are ejected from the core. Of the 19 brown dwarfs, 17 are ejected and only 2 remain in multiple systems in the core. The mean number of objects formed per simulation is 4.75.

4.6 ${\alpha _{turb} = 0.25}$

In the ten simulations with $\alpha _{\rm turb}=0.25$, a total of 60 objects are produced, 49 stars and 11 brown dwarfs, with each simulation producing between 3 and 10 objects. Of the 49 stars, 36 remain in the core in multiple systems, and the rest are ejected from the core. All 11 brown dwarfs are ejected from the core.

Table 2: A summary of the results of the simulations with $\alpha_{\rm turb} = 0.05,~0.10\;{\rm and}\;0.25~$, at time $t = 0.3~{\rm Myr}~$. Column 1 gives the simulation identifier and Col. 2 gives $\alpha _{\rm turb}$. Column 3 gives $M_{\rm obj}$, the total mass of objects formed (stars plus brown dwarfs), Col. 4 gives ${\cal N}_{\rm obj}$, the total number of objects formed, and Col. 5 gives ${\cal N}_{\rm BD}$, the total number of brown dwarfs formed. Column 6 gives the multiplicities of the multiple systems formed, and the final column (Col. 7) gives the mass of each individual object. Those objects which are part of a binary system are distinguished with b, those which are part of a triple system are distinguished with a t, and those which are part of a quadruple system (or in one case a quintuple system) are distinguished with q.

5 Statistics

Details of the numbers, masses and multiplicities of the objects produced in each of the fifty simulations with $\alpha_{\rm turb} = 0.05,~0.10\;{\rm and}\;0.25$ are shown in Table 2. The simulations with low turbulence ( $\alpha_{\rm turb} \leq 0.025$) are omitted because almost all of them produce only a single object, and therefore the discussion will now concentrate on the simulations with $\alpha_{\rm turb} = 0.05,~0.10\;{\rm and}\;0.25~$.

5.1 Numbers of objects and formation timescale

Increasing the level of turbulence has several effects. In the first instance, it delays somewhat the time at which objects are formed, and at the same time it increases the average number of objects formed.

With $\alpha _{\rm turb} = 0.05$, the first object (hereafter the primary protostar) forms $0.05~{\rm Myr}$ to $0.06~{\rm Myr}$ after the start of the simulation, and ${\sim} 95\%$ of all the other objects have formed by $0.12~{\rm Myr}$. With $\alpha _{\rm turb}=0.25$, the primary protostar forms $0.06~{\rm Myr}$ to $0.08~{\rm Myr}$ after the start of the simulation, and objects continue forming up to $0.15~{\rm Myr}$; this is because the extra turbulent energy gives the core extra support, and therefore delays its collapse.

The majority of secondary objects form in a dense disc-like slab around the primary protostar. The instabilities which produce these secondary objects are usually seeded - and propelled into the non-linear condensation regime - by the inhomogeneities in the accretion flow onto the slab. Higher levels of turbulence result in larger inhomogeneities, and hence in more secondary objects.

For $\alpha_{\rm turb} \leq 0.01$ the mean number of objects formed is $\left<{\cal N}_{\rm obj}\right> = 1$ (i.e. there are no secondary objects), whereas for $\alpha_{\rm turb} = 0.25~$, $~\left<{\cal N} _{\rm obj}\right> = 6$ (i.e. there are on average five secondary objects). Values of $\alpha _{\rm turb}$ and $\left <{\cal N}_{\rm obj}\right >$ are listed in Table 1, and plotted in Fig. 2, where the error bars give the standard deviation, i.e.

 \begin{displaymath} \pm\left[\left<{\cal N}_{\rm obj}^2\right> -\left<{\cal N}_{\rm obj}\right>^2\right]^{1/2} . \end{displaymath} (5)


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{02... ...ludegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig3c.ps}\end{figure} Figure 3: The normalized mass functions for $\alpha _{\rm turb} = 0.05$  (top), 0.10  (middle) and 0.25  (bottom). The shaded regions show stars in stable multiple systems, the hashed regions show stars in unstable multiple systems, and the open regions show single stars.
Open with DEXTER


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig4.ps}\end{figure} Figure 4: The combined (and un-normalized) mass function for all the simulations having $\alpha _{\rm turb} \geq 0.05$, (i.e. $\alpha_{\rm turb} = 0.05~(\times 20)$, $0.10~(\times 20)\;{\rm and}\;0.25~(\times 10)$), this being an approximation to the distribution of $\alpha _{\rm turb}$ values in the Taurus star formation region. The error bars give ${\pm}\!\sqrt{N}$ uncertainties. The dotted lines show separate fits to the unbound objects (the flat distribution at low mass) and the bound objects (the log-normal distribution at high mass); parameter fits are given in the main text). The dashed line shows the combination of these two separate fits.
Open with DEXTER

5.2 The mass function

The normalized mass functions (MF) for the ensembles with $\alpha_{\rm turb} = 0.05,~0.10\;{\rm and}\;0.25$ are shown in Fig. 3. In each case the shaded region shows the MF of objects in multiple systems and the open region shows the MF of single objects. The hashed region in the top figure shows an unstable triple that formed late on in one simulation, and is therefore likely to decay into a binary and an ejected single (see Paper I for more details).

The MFs are clearly very similar for each value of $\alpha _{\rm turb}$, viz. a high-mass peak of predominantly multiple stars and a low-mass tail of ejected stars and brown dwarfs. However, as $\alpha _{\rm turb}$ increases, the MF shifts slightly to lower masses. There are two reasons for this. (a) For higher $\alpha _{\rm turb}$ the overall collapse is delayed by the extra turbulent support, and therefore when the simulations are terminated at $0.3~{\rm Myr}$ less mass has been incorporated into objects (the mean mass incorporated into objects is given in the third column of Table 1); (b) for higher $\alpha _{\rm turb}$, the accretion flow onto the disc around the primary protostar is lumpier, so more objects are formed but individually they are less massive.

The combined mass function for all simulations with $\alpha _{\rm turb} \geq 0.05$ is shown in Fig. 4. The distribution of high-mass (predominantly bound) stars is well fitted by a log-normal distribution having mean $\left\langle \log_{10}~[M]\right\rangle = 0.05$ and standard deviation $\sigma_{\log_{10}~[M]} = 0.04$. The distribution of low-mass (predominantly unbound) objects is consistent with being flat in log-space from our resolution limit at ${\sim} 0.025 ~M_\odot$ up to ${\sim} 0.5 ~M_{\odot}$, above which it declines. The high-mass (predominantly bound) stars have an average mass of ${\sim} 1 ~M_{\odot}$ because, after ejections have removed some objects, there are usually two to four stars left in the core, and they are then able to accrete a total of ${\sim} 3 ~M_{\odot}$ between them. A more massive core would spawn more massive stars (Goodwin et al., in preparation).

The fraction of objects which are brown dwarfs is ${\cal N}_{\rm BD} / {\cal N}_{\rm obj} \sim 18\%$, and there does not appear to be a systematic dependence on the level of turbulence. This is somewhat higher than in Taurus ( ${\cal N}_{\rm BD}/{\cal N}_{\rm obj} \sim 13 \%$; Briceño et al. 2002), and somewhat lower than in Orion ( ${\cal N}_{\rm BD}/{\cal N}_{\rm obj}\sim 26 \%$; Muench et al. 2002). The fraction of low-mass objects ( $M < 0.5 ~M_\odot$) also appears to be independent of $\alpha _{\rm turb}$, and approximately  ${\sim} 50\%$.

5.3 Companion star frequencies

In analyzing the multiplicity of the objects formed in our simulations, we define "systems'' to include single objects, and "multiple systems'' to include only systems containing more than one object. The primary is the most massive star in a system; in a single it is the only star. Thus, if S is the number of single objects and B, T, Q and Q' are the numbers of binary, triple, quadruple and quintuple systems, respectively, the total number of objects is (S+2B+3T+4Q+5Q'+...), the total number of systems is (S+B+T+Q+Q'+...) (which is the same as the total number of primaries), and the total number of multiple systems is (B+T+Q+Q'+...).

Many different statistics have been introduced as measures of stellar (and brown dwarf) multiplicity (e.g., Reipurth & Zinnecker 1993), and they all reflect slightly different things[*]. They can be divided into two groups.

The first group of measures is normalized to the total number of objects, and is useful because it is straightforward to derive these measures for a subset of objects (for example, low-mass stars, in the range $0.08 ~M_\odot < M < 0.5 ~M_\odot$). The companion probability (Reipurth & Zinnecker 1993),

\begin{displaymath}{cp} = \frac{2B + 3T + 4Q + 5Q' + ... } {S + 2B + 3T + 4Q + 5Q' + ... } , \end{displaymath} (6)

gives the probability of an object having at least one companion, or equivalently the fraction of objects which is in multiple systems. However, it gives no indication of whether objects with companions are in binaries, or triples, or higher multiples. Therefore we prefer the companion frequency,


\begin{displaymath}{cf} = \frac{2B + 6T + 12Q + 20Q' + ... } {S + 2B + 3T + 4Q + 5Q' + ... } , \end{displaymath} (7)

which gives the mean number of companions per object.

The second group of measures is normalized to the total number of systems or the total number of multiple systems. The multiplicity frequency (RZ93),


\begin{displaymath}{mf} = \frac{B + T + Q + Q' + ... } {S + B + T + Q + Q' + ... } , \end{displaymath} (8)

gives the fraction of systems which is multiple. The pairing factor (RZ93),


\begin{displaymath}{pf} = \frac{B + 2T + 3Q + 4Q' + ... } {B + T + Q + Q' + ... } , \end{displaymath} (9)

gives the mean number of orbits per multiple system, or equivalently the mean number of companions per primary in multiple systems.

Table 1 records the values of all these measures for cores having $\alpha_{\rm turb} = 0.05,~0.10\;{\rm and}\;0.25$. We have also calculated the companion frequency, cf, separately for the low-mass objects ( ${<}0.5 ~M_{\odot}$) and the high-mass objects ( ${>}0.5 ~M_{\odot}$). The fraction of objects in multiple systems (i.e. the companion probability, cp) increases steadily with increasing $\alpha _{\rm turb}$. The mean number of companions (i.e. the companion frequency, cf) increases even more rapidly with increasing $\alpha _{\rm turb}$. cf is always larger for high-mass objects than for low-mass objects (i.e. there are fewer low-mass objects in multiples than high-mass objects), but the rate of increase of cf with increasing $\alpha _{\rm turb}$ is greater for the low-mass objects.

The smaller companion frequency for low-mass objects is due to the fact that low-mass objects have usually been ejected from the core before they could accrete much mass (that is why they have low mass), and ejected stars tend to be singles. Stars that end up in stable multiples also tend to remain in the core, and therefore they grow to larger masses by continuing to accrete.

The multiplicity of low-mass stars increases with increasing $\alpha _{\rm turb}$, because a greater number of higher-order multiples is formed in cores with higher $\alpha _{\rm turb}$. For example, when $\alpha _{\rm turb} = 0.05$, only $10\%$ of simulations produce quadruples, but this fraction rises to $60\%$ for $\alpha _{\rm turb}=0.25$. In higher-order multiples, the low-mass objects tend to be outlying members. They are less able to accrete from the remaining gas, and they tend to remain low-mass.


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{02... ...ludegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig5c.ps}\end{figure} Figure 5: The histograms show the cumulative distribution functions of semi-major axes for the ensembles with $\alpha _{\rm turb} = 0.05$  (top), 0.10  (middle) and 0.25  (bottom). The dashed line shows the Gaussian fit to the DM91 period distribution.
Open with DEXTER


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.45cm,angle=270,clip]{0285fig6.ps}\end{figure} Figure 6: For each binary system, the semi-major axis, a is plotted against the number of objects formed in that simulation, $N_{\rm obj}$, with open circles for for $\alpha _{\rm turb} = 0.05$, filled circles for $\alpha _{\rm turb}=0.10$, and stars for $\alpha _{\rm turb}=0.25$. The average semi-major axis, $\bar{a}$, decreases with increasing $N_{\rm obj}$, due to dynamical hardening. To the left of the dashed line at $a = 5~{\rm au}$ the a-values are upper limits, due to gravity softening.
Open with DEXTER

5.4 Separations

Figure 5 shows the cumulative distribution functions (CDFs) of the semi-major axes of multiple systems, for different initial levels of turbulence, $\alpha_{\rm turb} = 0.05,\;0.10\; {\rm and}\;0.25$. For comparison the Gaussian fit to the DM91 local G-dwarf sample is plotted as a dashed line.

As noted in Paper I, the semi-major axis distribution for the $\alpha _{\rm turb} = 0.05$ ensemble is consistent with the DM91 observations. In contrast, the semi-major axis distributions for the $\alpha_{\rm turb} = 0.10\;{\rm and}\;0.25$ ensembles both have too many hard binaries (a < 20 au), and both are rejected by the KS test as being drawn from the DM91 fit, at ${>}90\%$ confidence. A similar excess of hard binaries is predicted by the core fragmentation simulations of Delgado-Donate et al. (2003, 2004), which invoke even higher levels of turbulence ( $\alpha_{\rm turb} = 1$), and by the N-body simulations of Sterzik & Durisen (2003).

As described in Paper I, hard binaries are formed primarily by few-body interactions, including those which eject low mass objects. Consequently, in simulations where larger numbers of objects are formed, the binaries are on average harder (regardless of $\alpha _{\rm turb}$). For example, when few objects are formed, say $N_{\rm obj} \leq 3$, the average separation of binaries is ${>}100~{\rm au}$, whereas when $N_{\rm obj} = 4$ the average separation is ${\sim} 30~{\rm au}$, and when $N_{\rm obj} \geq 5$ it is ${\sim} 20~{\rm au}$. Figure 6 shows the semi-major axes of all systems plotted against the number of objects formed in that simulation. There is clearly a trend of decreasing semi-major axis with increasing number of objects.

Since the gravitational forces between objects are kernel softened with a smothing length h equal to the sink radius $R_{\rm sink} = 5~{\rm au}$, orbits with small semi-major axes ( $a < 5~{\rm au}$) are also softened. It follows that the distribution of semi-major axes below $5~{\rm au}$ is distorted. Given that the code conserves angular momentum very accurately, we infer that these already hard orbits should be even harder. For $\alpha_{\rm turb} = 0.10\;{\rm and}\;0.25$ this would exacerbate the difference between the numerically derived distribution of semi-major axes and the observations of DM91. Conversely, for $\alpha _{\rm turb} = 0.05$, it would improve the agreement with the DM91 distribution.


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig7.ps}\end{figure} Figure 7: The distribution of mass ratios for simulations with $\alpha _{\rm turb} = 0.05$ (filled), $\alpha _{\rm turb}=0.10$ (hashed) and $\alpha _{\rm turb}=0.25$ (open). (Note that the number of simulations in the $\alpha _{\rm turb}=0.25$ ensemble is only half of the number in the other ensembles.)
Open with DEXTER

5.5 Mass ratios

Figure 7 shows the distribution of binary mass ratios, q = M2 / M1, for all simulations having $\alpha _{\rm turb} \geq 0.05$. Note that high q means $q \sim 1$, i.e. components of comparable mass.

For $\alpha _{\rm turb} = 0.05$, the distribution of mass ratios is quite flat, and reminiscent of the observed distribution for local G dwarfs (DM91; Mazeh et al. 1992).

For $\alpha = 0.10\;{\rm and}\;0.25$, the distribution is dominated by high-q close binaries. All binaries in these ensembles have semi-major axes $a < 100~{\rm au}$ (the high-a tail is produced by wider orbits in higher-order systems), and $64\%$ of these have q > 0.8. This is very similar to the mass ratio distribution observed in Taurus-Auriga by WG01, who found that over ${\sim} 60\%$ of binaries with separations ${<} 100~{\rm au}$ had q > 0.8. In our simulations, systems with high mass ratio tend to be close (all systems with q > 0.7, and most systems with q > 0.4, are binaries with $a < 20~ {\rm au}$), but the reverse is not always true: in other words, there are a few close binaries with low mass ratios. Close binaries are presumed to acquire high mass ratios because the material accreting onto the system has high specific angular momentum, and is therefore more readily accommodated by the secondary (Whitworth et al. 1995; Bate & Bonnell 1997; Paper I).


  \begin{figure} \par\includegraphics[height=12.1cm,width=7.4cm,angle=270,clip]{0285fig8.ps}\end{figure} Figure 8: The cumulative distribution function of eccentricities for all simulations with $\alpha _{\rm turb} \geq 0.05$. The dashed line shows the DM91 fit of a linearly increasing eccentricity distribution.
Open with DEXTER

5.6 Eccentricities

Figure 8 shows the CDF of eccentricity for all the simulations having $\alpha _{\rm turb} \geq 0.05$ and the linear fit to the observed distribution proposed by DM91. The two distributions are consistent.

6 Discussion

6.1 Minimum level of turbulence for multiple formation

When the initial level of turbulence is low, $\alpha_{\rm turb} \leq 0.01$, it seems that a core can only spawn a single central star. Even for $\alpha_{\rm turb} = 0.025$, the core is unlikely to spawn a multiple system. We therefore focus our discussion on the higher levels of turbulence, $\alpha_{\rm turb} = 0.05,~0.10\;{\rm and}\;0.25$, for which multiple star formation is the norm. In this range a number of significant systematic trends are evident.

6.2 Time-scale for star formation

As $\alpha _{\rm turb}$ increases from 0.10 to 0.25, the average timescale for star formation increases somewhat, due to the extra support which turbulence affords the core. For $\alpha _{\rm turb} = 0.05$, the primary protostar forms after $0.05~{\rm Myr}$, and most of the secondary protostars have formed by $0.12~{\rm Myr}$. For $\alpha_{\rm turb} = 0.10\;{\rm and}\;0.25$, the primary protostar forms after $0.06~{\rm Myr}$, and most of the secondary protostars have formed by $0.15~{\rm Myr}$. In only one case (run 073) do objects form after 0.25 Myr, and so it appears that the fragmentation phase is almost always over by the end of the simulations at 0.3 Myr. After the end of the simulations accretion will be on-going. However, feedback from the protostars is expected to become very important, possibly dispersing a significant fraction of the gas not already in stars or discs around them.

6.3 Number of objects formed

As $\alpha _{\rm turb}$ increases from 0.10 to 0.25, the average number of objects formed, $\left <{\cal N}_{\rm obj}\right >$, increases from $4.6 \pm 2.6$ to $6.0 \pm 2.0$ (see Fig. 2). This is because a higher level of turbulence generates more density contrast - i.e. more numerous and more compressed lumps - and therefore more protostars.

6.4 Masses of objects formed

As $\alpha _{\rm turb}$ increases from 0.10 to 0.25, the average mass of the objects formed decreases slightly (see Fig. 3). There are two factors involved here. First, as noted above a higher level of turbulence means that more objects are formed. Second, a higher level of turbulence means that the core has more support, and therefore a smaller fraction of its mass has condensed out after $0.3~{\rm Myr}$.

6.5 The bimodal mass function

Apart from this slight decrease in average mass with increasing turbulence, the form of the mass function appears to be independent of $\alpha _{\rm turb}$. Specifically, the mass function is bimodal: the lower-mass stars (which tend to be single stars ejected from the core) subscribe to a flat segment of the mass function; and the higher-mass stars (which tend to be those remaining in the core and pairing up in multiple systems) subscribe to a Gaussian segment of the mass function (see Fig.  4).

The critical mass seperating the two modes in the mass function, $M_{\rm crit} \sim 0.5 ~M_\odot$, arises because of the interplay between ejection by dynamical interaction and growth by accretion. Anosova (1986) has shown that the decay time for small-N systems is

\begin{displaymath}t_{\rm decay} \sim 100~t_{\rm cross} \sim 17~~ \left( \frac{R... ...ight)^{3/2}~ \left( \frac{M}{M_\odot} \right)^{-1/2}~{\rm yr}, \end{displaymath} (10)

and the ensemble of low-mass secondaries which forms in our cores typically has $R \sim 200~{\rm au}$ and $M \sim 1\;{\rm to}\;2 ~M_\odot$, so $t_{\rm decay} \sim 3\;{\rm to}\;5 \times 10^4~{\rm yr}$. The accretion rate is ${\sim} 10^{-5} ~M_\odot~{\rm yr}^{-1}$, and so the objects which get ejected have masses ${\la} 0.5 ~M_\odot$. The probability of ejection is only weakly dependent on mass, ${\propto} M^{-1/3}$ (Anosova 1986), and so the ejected objects have a flat mass function. In contrast, the two or three objects which survive the dynamical decay phase remain embedded in the centre of the core and compete for the gas which continues to fall into the centre, so they grow to ${\sim} M_\odot$.

6.6 Companion star frequency

As $\alpha _{\rm turb}$ increases from 0.10 to 0.25, the companion star frequency increases slightly for intermediate-mass stars ( $0.5 ~M_\odot$ to $5 ~M_\odot$), and quite markedly for low-mass objects ( ${<}0.5 ~M_{\odot}$).

6.7 Distribution of semi-major axes

For $\alpha _{\rm turb} = 0.05$, the distribution of semi-major axes is broad and indistinguishable from the distribution inferred for local G-dwarfs by DM91. There is a lack of very close systems with $a < 5~{\rm au}$. This is due to the fact that sinks have a finite size and their mutual gravity is softened; therefore our code cannot resolve very close systems.

In contrast, for $\alpha_{\rm turb} = 0.10\;{\rm and}\;0.25$, there are many more close systems ( $5~{\rm au} \la a \la 20~{\rm au}$) than in the DM91 sample (see Fig. 5), and this discrepancy would not be alleviated if the code were able to resolve very close systems.

Much of the hardening which produces close binaries is due to dynamical interactions with other objects, in particular with the low-mass objects which get ejected in the process. The excess of close systems produced by higher levels of turbulence may therefore be due to the greater number of objects formed, and hence the greater potential for dynamical interactions, as suggested by Fig. 6.

6.8 Distribution of mass-ratios

For $\alpha _{\rm turb} = 0.05$, the distribution of mass ratios is flat and indistinguishable from the distribution for local G-dwarfs reported by DM91. In contrast, for $\alpha_{\rm turb} = 0.10\;{\rm and}\;0.25$, there is an excess of systems having high mass-ratio, i.e. components of comparable mass (see Fig. 7). Many of these systems with high mass-ratio arise because the binary system has accreted material with relatively high specific angular momentum, and this material can more easily be accommodated by the secondary (e.g., Whitworth et al. 1995).

6.9 Close systems with comparable components

For $\alpha_{\rm turb} = 0.10\;{\rm and}\;0.25$, the systems with high mass-ratio tend also to be close, and it is this sub-population of high-mass-ratio close binaries which is the main difference between the distributions of semi-major axis and mass-ratio for the protostars formed in these simulations, and the distributions of semi-major axis and mass-ratio for local G-dwarfs as reported by DM91. A similar excess of close systems with comparable components was found by Delgado-Donate et al. (2003, 2004), who simulated the collapse and fragmentation of cores with even higher levels of turbulence ( $\alpha_{\rm turb} = 1$).

Taken at face value, this suggests that the local population of G-dwarfs must have been formed in cores with low turbulence ( $\alpha_{\rm turb} \sim 0.05$). However, this conclusion rests on the assumption that the spherically symmetric $5.4 ~M_\odot~$ core and the $P_k \propto k^{-4}$ turbulence spectrum which we have adopted, are representative of the cores forming G-dwarfs, and there is no firm basis for this assumption. An alternative explanation is that a significant population of close, high-mass-ratio systems has escaped detection, but we believe this to be unlikely.

A significant contrast to this is found in Taurus. Here WG01 find that binaries in the separation range $10~{\rm au} \la a \la 100~{\rm au}$ do indeed have significantly higher mass ratios than wider binaries. Therefore they are compatible with formation in cores having higher levels of turbulence, $\alpha_{\rm turb} = 0.10\; {\rm to}\;0.25$. We have discussed the origin of the mass function and the binary statistics in Taurus in Goodwin et al. (2004b).

7 Conclusions

We have explored the influence of turbulence on the fragmentation of dense molecular cores, by means of a large ensemble of simulations. In this ensemble, we consider a spherically symmetric $5.4 ~M_\odot$ core with a Plummer-like density profile; this is a good representation of observed cores like L1544. We seed the core with a turbulent velocity field having power spectrum  $P(k) \propto k^{-4}$. The number of objects that forms, and the properties of the resulting multiple systems depend both on the level of turbulence $\alpha _{\rm turb}$, and on the details of the turbulent velocity field. Therefore for each value of $\alpha _{\rm turb}$ we we have simulated many different realizations by changing the random number seed for the turbulent velocity field. The main conclusions are:

(i)
The formation of multiple systems requires $\alpha_{\rm turb} \ga 0.025$; a core with $\alpha_{\rm turb} = 0.025$ has a ${\sim} 20\%$ chance of forming a multiple system, and a core with $\alpha_{\rm turb} \ga 0.25$ almost always forms a multiple system.

(ii)
As $\alpha _{\rm turb}$ is increased, the average time-scale for object formation increases, the average number of objects formed increases, the companion frequency increases (particularly for the lower-mass objects), and the average mass of objects decreases.

(iii)
The mass function has a bimodal form. The low-mass objects, which are usually single because they have been ejected from the core before they could grow above $0.5 ~M_\odot$, subscribe to a flat segment of the mass function. The high-mass stars, which have usually stayed embedded in the core, grown by accretion and paired up with one another in multiple systems, subscribe to a Gaussian segment of the mass function. Typically $20\%$ of objects are brown dwarfs ( $M < 0.08 ~M_\odot$), and $50\%$ are low-mass stars ( $0.08 ~M_\odot < M < 0.5 ~M_\odot$).

(iv)
For $\alpha_{\rm turb} \ga 0.10$, there is a significant subpopulation of binary systems having small semi-major axes and high mass ratios, i.e. close systems with components of comparable mass. This subpopulation is also found in the simulations of Delgado-Donate et al. (2003, 2004) who treat the extreme case $\alpha_{\rm turb} = 1$. It is not present in the sample of local G-dwarfs observed by DM91, but there is some evidence for it in the Taurus pre-Main-Sequence sample observed by WG01.

(v)
The ensemble of simulations for cores with $\alpha _{\rm turb} = 0.05$ reproduces the binary statistics of the DM91 sample (companion-star frequency and distributions of semi-major axis, eccentricity and mass-ratio) very well. Therefore if the other core parameters, which we have not varied (e.g. mass), are representative of the cores forming local G dwarfs, we infer that these cores must have had finite but low levels of turbulence, $\alpha_{\rm turb} \sim 0.05$.

(vi)
Both the mass function, and the binary statistics, for the WG01 sample of pre-MS stars in Taurus are reproduced by a mix of simulations with $\alpha_{\rm turb} =$ 0.05 (20 realizations), 0.10 (20 realizations), and 0.25 (10 realizations), as shown by Goodwin et al. (2004b).

Acknowledgements

S.P.G. acknowledges support of PPARC grant PPA/G/S/1998/00623 and is now a UKAFF Fellow. We thank B. Sathyprakash and R. Balasubramanian for allowing us extensive use of the Beowulf cluster of the gravitational waves group at Cardiff; and Matthew Bate for helpful discussions and for providing the code to generate turbulence.

References

 
Copyright ESO 2004