format: text

text:
# The formation of a star cluster: predicting the properties of stars and brown dwarfs

## 1 Introduction

The collapse and fragmentation of molecular cloud cores to form bound multiple stellar systems has been the subject of many numerical studies (e.g. Boss & Bodenheimer 1979; Boss 1986; Bonnell et al. 1991; Nelson & Papaloizou 1993; Bonnell 1994; Burkert & Bodenheimer 1993; Bate, Bonnell & Price 1995; Truelove et al. 1998). These calculations have resulted in the adoption of fragmentation as the favoured mechanism for the formation of binary and multiple stars, since it can produce a wide range of binary properties through simple variations of the pre-collapse initial conditions.

However, while individual binary systems can be reproduced by such fragmentation calculations, it is extremely difficult to use these calculations to predict the statistical properties of the stellar systems that should result from the fragmentation model. Quantities that we may wish to determine include the initial mass function (IMF), the relative frequencies of single, binary and multiple stars, the properties of multiple stars, the properties of circumstellar discs, and the efficiency of star formation. In order to predict these statistical properties, we need to produce a large sample of stars. There are two possibilities.

We could perform many calculations of isolated cloud cores using a representative sample of initial conditions. However, this has two disadvantages. First, the conditions in molecular clouds are not sufficiently well understood to be able to select a representative sample of cloud cores for the initial conditions. Second, the production of isolated stellar systems neglects interactions between systems that may be important in determining stellar properties, especially in young star clusters. Examples of such interactions include binary formation via star-disc capture (Larson 1990; Clarke & Pringle 1991a,b; Heller 1995; McDonald & Clarke 1995; Hall, Clarke & Pringle 1996), truncation of protostellar discs (Heller 1995; Hall 1997), and competitive accretion leading to a range of stellar masses (Larson 1978; Zinnecker 1982; Bonnell et al. 1997, 2001a,b).

The second possibility is to perform a calculation of the collapse and fragmentation of a large-scale molecular cloud to form many stars simultaneously. This is the approach taken in this paper. Interactions between stars are automatically allowed for. We must still specify global initial conditions, but the formation of individual cores within the cloud occurs self-consistently; we do not have to select initial conditions for each core arbitrarily. The only disadvantage is that such a calculation is extremely computationally intensive.

Several such global calculations have been performed in the past. The earliest was that of Chapman et al. \[Chapman et al. 1992\] who followed the collapse and fragmentation of a shock-compressed layer of molecular gas between two colliding clouds. The calculation produced many single, binary, and multiple protostars, but there was no attempt to derive the statistical properties of these systems. Klessen, Burkert & Bate \[Klessen, Burkert, & Bate 1998\] followed the collapse of a large-scale clumpy molecular cloud to form $`60`$ protostars. They found that the mass function of the protostars could be fit by a lognormal mass function that has a similar width to the observed stellar initial mass function. The protostellar masses were set by a combination of the initial density structure, competitive accretion, and dynamical interactions. Further calculations in which the global initial conditions were varied confirmed the lognormal form of the mass function and showed that the mean mass of the protostars was similar to the mean initial Jeans mass in the cloud (Klessen & Burkert 2000, 2001; Klessen 2001). These calculations enabled us to identify some of the processes that may help to determine the initial mass function. However, they did not have the resolution to follow the collapsing molecular gas all the way down to the opacity limit for fragmentation (Low & Lynden-Bell 1976; Rees 1976; Silk 1977a), or even to resolve the median separation of binary systems of $`30`$ AU \[Duquennoy & Mayor 1991\]. Thus, we cannot determine the total number of stars that will form or the stellar initial mass function from these calculations, let alone the frequency of binary stars or the importance of star-disc encounters.

This paper presents results from the first calculation to follow the collapse and fragmentation of a large-scale turbulent molecular cloud to form a stellar cluster while resolving beyond the opacity limit for fragmentation. Thus, assuming that fragmentation does not occur at densities greater than those at which the opacity limit sets in (Section 2), it resolves all potential fragmentation, including that which produces binary systems. This allows us to predict a wide variety of stellar properties. Two papers that contain results from this calculation have already been published. They concentrate on the formation mechanisms of brown dwarfs (Bate, Bonnell & Bromm 2002a) and close binaries (Bate, Bonnell & Bromm 2002b). In this paper, we consider how the dynamics of star formation determine the properties of stars and brown dwarfs, and we compare these properties with observations.

The outline of this paper is as follows. In Section 2, we briefly review the opacity limit for fragmentation. The computational method and the initial conditions for our calculation are described in Section 3. Section 4 discusses the evolution of the cloud and the star formation that occurs during the calculation. The properties of the resulting stars and brown dwarfs are compared with observations of star-forming regions in Section 5. Finally, in section 6, we give our conclusions

## 2 The opacity limit for fragmentation

When a molecular cloud core begins to collapse under its own gravity the gravitational potential energy that is released can easily be radiated away so that the collapsing gas is approximately isothermal (e.g. Larson 1969). Thus, the thermal pressure varies with density $`\rho `$ as $`p\rho ^\eta `$ where the effective polytropic exponent $`\eta \mathrm{d}\mathrm{log}[p]/\mathrm{d}\mathrm{log}[\rho ]1`$. This allows the possibility of fragmentation because the Jeans mass decreases with increasing density if $`\eta <4/3`$.

The opacity limit for fragmentation occurs when the rate at which energy is released by the collapse exceeds the rate at which energy can be radiated away (Rees 1976; Low & Lynden-Bell 1976; Masunaga & Inutsuka 1999). The gas then heats up with $`\eta >4/3`$, the Jeans mass increases, and a Jeans-unstable collapsing clump quickly becomes Jeans-stable so that a pressure-supported fragment forms. The density at which this occurs depends on the opacity (and the initial temperature) of the gas (Masunaga & Inutsuka 1999), hence the term ‘opacity limit for fragmentation’. For standard molecular gas at an initial temperature of 10 K the gas begins to heat significantly at a density of $`10^{13}\mathrm{g}\mathrm{cm}^3`$ (Larson 1969; Masunaga & Inutsuka 2000).

The pressure-supported fragment initially contains several Jupiter-masses (MJ) and has a radius of $`5`$ AU \[Larson 1969\]. Such a fragment is expected to be embedded within a collapsing envelope. Thus, its mass grows with time and its central temperature increases. When its central temperature reaches 2000 K, molecular hydrogen begins to dissociate. This provides a way for the release of gravitational energy to be absorbed without significantly increasing the temperature of the gas. Thus, a nearly isothermal second collapse occurs within the fragment that ultimately results in the formation of a stellar core with radius $`1R_{}`$ \[Larson 1969\].

Several studies have investigated the possibility of fragmentation during this second collapse (Boss 1989; Bonnell & Bate 1994; Bate 1998, 2002). Boss \[Boss 1989\] found that fragmentation was possible during this second collapse, but that the objects spiralled together and merged due to gravitational torques from non-axisymmetric structure. Bonnell & Bate \[Bonnell & Bate 1994\] found that fragmentation to form close binaries and multiple systems could occur in a disc that forms around the stellar core. However, both these studies began with somewhat arbitrary initial conditions and modelled only the pressure-supported fragment. Bate \[Bate 1998\] performed the first three-dimensional calculations to follow the entire collapse of a molecular cloud core through the formation of the pressure-supported fragment, the second collapse phase, and the formation of the stellar core and its surrounding disc. In these and subsequent calculations (Bate 2002), it was found that the second collapse did not result in sub-fragmentation due to the high degree of thermal pressure and angular momentum transport via gravitational torques from non-axisymmetric structure.

Therefore, it appears that the opacity limit for fragmentation is real and that fragmentation cannot occur at densities exceeding $`10^{13}\mathrm{g}\mathrm{cm}^3`$. The implications are that there should be a minimum ‘stellar’ mass of $`10`$ MJ (Low & Lynden-Bell 1976; Silk 1977a, Boss 1988) and that protobinaries should have a minimum initial separation of $`10`$ AU due to the sizes of the pressure-supported fragments. The exact value of the minimum mass is uncertain with theoretical values in the literature ranging from $`110`$ MJ (Low & Lynden-Bell 1976; Silk 1977a; Boss 1988; Masunaga & Inutsuka 1999; Boss 2001). Surveys of young star clusters are beginning to probe masses down to this theoretical minimum mass (Zapatero Osorio et al. 1999; Lucas & Roche 2000; Béjar et al. 2001; Martín et al. 2001b; Lucas et al. 2001), with the masses of some objects estimated to be as low as 3 Jupiter-masses (Zapatero Osorio et al. 2002a; McCaughrean 2003). However, the observational uncertainties are still large and a cutoff in the mass function has yet to be detected.

## 3 Computational method

The calculation presented here was performed using a three-dimensional, smoothed particle hydrodynamics (SPH) code. The SPH code is based on a version originally developed by Benz (Benz 1990; Benz et al. 1990). The smoothing lengths of particles are variable in time and space, subject to the constraint that the number of neighbours for each particle must remain approximately constant at $`N_{\mathrm{neigh}}=50`$. The SPH equations are integrated using a second-order Runge-Kutta-Fehlberg integrator with individual time steps for each particle \[Bate, Bonnell, & Price 1995\]. Gravitational forces between particles and a particle’s nearest neighbours are calculated using a binary tree. We use the standard form of artificial viscosity (Monaghan & Gingold 1983; Monaghan 1992) with strength parameters $`\alpha __\mathrm{v}=1`$ and $`\beta _\mathrm{v}=2`$. Further details can be found in Bate et al. \[Bate, Bonnell, & Price 1995\]. The code has been parallelised by M. Bate using OpenMP.

### 3.1 Equation of state

To model the opacity limit for fragmentation, discussed in Section 2, without performing radiative transfer, we use a barotropic equation of state for the thermal pressure of the gas $`p=K\rho ^\eta `$, where $`K`$ is a measure of the entropy of the gas. The value of the effective polytropic exponent $`\eta `$, varies with density as

$$\eta =\{\begin{array}{cc}\begin{array}{cc}\hfill 1,& \rho 10^{13}\mathrm{g}\mathrm{cm}^3,\hfill \\ \hfill 7/5,& \rho >10^{13}\mathrm{g}\mathrm{cm}^3.\hfill \end{array}\hfill & \end{array}$$
(1)
We take the mean molecular weight of the gas to be $`\mu =2.46`$. The value of $`K`$ is defined such that when the gas is isothermal $`K=c_\mathrm{s}^2`$, with the sound speed $`c_\mathrm{s}=1.84\times 10^4`$ cm s-1 at 10 K, and the pressure is continuous when the value of $`\eta `$ changes.

This equation of state has been chosen to match closely the relationship between temperature and density during the spherically-symmetric collapse of molecular cloud cores as calculated with frequency-dependent radiative transfer (Masunaga, Miyama, & Inutsuka 1998; Masunaga & Inutsuka 2000). A comparison of our simple parameterisation with Masunaga and Inutsuka’s temperature-density relation is given in Figure 1. Our parameterisation reproduces the temperature-density relation to an accuracy of better than 20% in the non-isothermal regime up to densities of $`10^8`$ g cm-3. Test calculations of the spherically-symmetric collapse of a molecular cloud core using this equation of state have been performed (Bate 1998, 2002) and give excellent agreement with the results of Larson \[Larson 1969\] and Winkler & Newman (1980a,b) for the mass and size of the pressure-supported fragment that forms. Thus, our equation of state should model collapsing regions well, but may not model the equation of state in protostellar discs particularly accurately due to the departure from spherical symmetry.

### 3.2 Sink particles

The opacity limit for fragmentation results in the formation of distinct pressure-supported fragments in the calculation. As these fragments accrete, their central density increases, and it becomes computationally impractical to follow their internal evolution until they undergo the second collapse to form stellar cores because of the short dynamical time-scales involved. Therefore, when the central density of a pressure-supported fragment exceeds $`\rho _\mathrm{s}=10^{11}\mathrm{g}\mathrm{cm}^3`$, we insert a sink particle into the calculation Bate et al. \[Bate, Bonnell, & Price 1995\].

In the calculation presented here, a sink particle is formed by replacing the SPH gas particles contained within $`r_{\mathrm{acc}}=5`$ AU of the densest gas particle in a pressure-supported fragment by a point mass with the same mass and momentum. Any gas that later falls within this radius is accreted by the point mass if it is bound and its specific angular momentum is less than that required to form a circular orbit at radius $`r_{\mathrm{acc}}`$ from the sink particle. Thus, gaseous discs around sink particles can only be resolved if they have radii $`\text{ }>10`$ AU. Sink particles interact with the gas only via gravity and accretion.

Since all sink particles are created from pressure-supported fragments, their initial masses are $`10`$ MJ, as given by the opacity limit for fragmentation (Section 2). Subsequently, they may accrete large amounts of material to become higher-mass brown dwarfs ($`\text{ }<75`$ MJ) or stars ($`\text{ }>75`$ MJ), but all the stars and brown dwarfs begin as these low-mass pressure-supported fragments.

The gravitational acceleration between two sink particles is Newtonian for $`r4`$ AU, but is softened within this radius using spline softening \[Benz 1990\]. The maximum acceleration occurs at a distance of $`1`$ AU; therefore, this is the minimum separation that a binary can have even if, in reality, the binary’s orbit would have been hardened. Sink particles are not permitted to merge in this calculation.

Replacing the pressure-supported fragments with sink particles is necessary in order to perform the calculation. However, it is not without an element of risk. If it were possible to follow the fragments all the way to stellar densities (as done by Bate 1998) while continuing to follow the evolution of the large-scale cloud over its dynamical time-scale, we might find that a few of the objects that we replace with sink particles merge together or are disrupted by dynamical interactions. We have tried to minimise the degree to which this might occur by insisting that the central density of the pressure-supported fragments exceeds $`\rho _\mathrm{s}`$ before a sink particle is created. This is two orders of magnitude higher than the density at which the gas is heated and ensures that the fragment is self-gravitating, centrally-condensed and, in practice, roughly spherical before it is replaced by a sink particle (see, for example, the fragments that have not yet been replaced by sink particles in Figure 5 at $`t=1.337`$). In theory, it would be possible for a long collapsing filament to exceed this density over a large distance, thus making the creation of one or more sink particles ambiguous. However, the structure of the collapsing gas that results from the turbulence prohibits this from occurring; no long, roughly uniform-density filaments with densities $`\rho _\mathrm{s}`$ form during the calculation. Furthermore, each pressure-supported fragment must undergo a period of accretion before its central density exceeds $`\rho _\mathrm{s}`$ and it is replaced by a sink particle. For example, it is common in the calculation to be able to follow a pressure-supported fragment that forms via gravitational instability in a disc for roughly half an orbital period before it is replaced. Thus, the fragments do have some time in which they may merge or be disrupted. Only occasionally are low-mass pressure supported fragments disrupted during the calculation; most are eventually replaced by sink particles.

### 3.3 Initial conditions

The initial conditions consist of a large-scale, turbulent molecular cloud. The cloud is spherical and uniform in density with mass of 50 M and a diameter of $`0.375`$ pc (77400 AU). At the temperature of 10 K, the mean thermal Jeans mass is 1 M (i.e. the cloud contains 50 thermal Jeans masses). The free-fall time of the cloud is $`t_{\mathrm{ff}}=6.0\times 10^{12}`$ s or $`1.90\times 10^5`$ years.

Although the cloud is uniform in density, we impose an initial supersonic turbulent velocity field on it in the same manner as Ostriker, Stone & Gammie \[Ostriker, Stone, & Gammie 2001\]. We generate a divergence-free random Gaussian velocity field with a power spectrum $`P(k)k^4`$, where $`k`$ is the wavenumber. In three dimensions, this results in a velocity dispersion that varies with distance, $`\lambda `$, as $`\sigma (\lambda )\lambda ^{1/2}`$ in agreement with the observed Larson scaling relations for molecular clouds \[Larson 1981\]. This power spectrum is slightly steeper than the Kolmogorov spectrum, $`P(k)k^{11/3}`$. Rather, it matches the amplitude scaling of Burgers supersonic turbulence associated with an ensemble of shocks (but differs from Burgers turbulence in that the initial phases are uncorrelated). The velocity field is generated on a $`64^3`$ uniform grid and the velocities of the particles are interpolated from the grid. The velocity field is normalised so that the kinetic energy of the turbulence equals the magnitude of the gravitational potential energy of the cloud (i.e. initially, the cloud has more than enough turbulent energy to support itself against gravity). The initial root-mean-square Mach number of the turbulence is $`=6.4`$.

### 3.4 Resolution

The local Jeans mass must be resolved throughout the calculation to model fragmentation correctly (Bate & Burkert 1997; Truelove et al. 1997; Whitworth 1998; Boss et al. 2000). Bate & Burkert \[Bate & Burkert 1997\] found that this requires $`\text{ }>2N_{\mathrm{neigh}}`$ SPH particles per Jeans mass; $`N_{\mathrm{neigh}}`$ is insufficient. We have repeated their calculation using different numbers of particles and find that $`1.5N_{\mathrm{neigh}}=75`$ particles is also sufficient to resolve fragmentation (see Appendix A). The minimum Jeans mass in the calculation presented here occurs at the maximum density during the isothermal phase of the collapse, $`\rho =10^{13}`$ g cm-3, and is $`0.0011`$ M (1.1 MJ). Thus, we use $`3.5\times 10^6`$ particles to model the 50-M cloud. This SPH calculation is one of the largest ever performed. It required approximately 95000 CPU hours on the SGI Origin 3800 of the United Kingdom Astrophysical Fluids Facility (UKAFF).

## 4 The evolution of the cloud

Although the initial turbulent velocity field is divergence-free, hydrodynamic evolution of the cloud soon results in the formation of shocks. Small-scale shocks form first, followed by larger structures as the calculation proceeds (Figure 2, $`t=00.8`$ tff). Kinetic energy is lost from the cloud in these shocks. Therefore, although the cloud initially has more than enough turbulent energy to support itself, gravity soon begins to dominate and collapse occurs in parts of the cloud. The turbulence decays on the dynamical time-scale of the cloud, and star formation begins after just one global free-fall time at $`t=1.037t_{\mathrm{ff}}`$ (i.e. $`t=1.97\times 10^5`$ yrs). By this time, the root-mean-square Mach number of the turbulence has fallen from its initial value of $`=6.4`$ to only $`=3.8`$. This rapid decay of the turbulence is consistent with other numerical studies of turbulence in molecular clouds, both with and without magnetic fields and self-gravity (e.g. MacLow et al. 1998; Stone et al. 1998; Ostriker et al. 2001).

Regions of overdensity form from converging gas flows during the decay of the turbulence collapse (Figure 2, $`t0.8`$ tff), producing three dense star-forming cores within the cloud (Figure 2, $`t=1.0`$ and $`t=1.2`$ tff; Figures 3, 4, 5, 6). The properties of these cores and the numbers of stars and brown dwarfs they produce during the calculation are summarised in Table 1. Figure 7 gives the times at which each star or brown dwarf is formed (i.e. the time at which the pressure-supported fragments are replaced by sink particles) in terms of the number of free-fall times from the start of the calculation (and also in years).

The most massive core begins forming stars first (Figures 3, 4). It contains $`3.0`$ M when star formation begins ($`6`$% of the mass of the cloud), although this figure depends on the density threshold that is used (Table 1). The core undergoes an initial burst of star formation lasting $`18000`$ years followed by a quiet period of $`24000`$ years during which only three brown dwarfs form (see Figure 7). This pattern then repeats with another burst lasting $`18000`$ years and a quiet period that lasts the remaining $`9000`$ years until the end of the calculation, during which time only one brown dwarf forms in the most massive core. At this point we had exhausted our allocation of computer time and the calculation was stopped at $`t=1.40`$ tff ($`t=2.66\times 10^5`$ years). In all, 38 stars and brown dwarfs formed in the most massive core.

The second and third dense cores begin forming stars at $`t=1.296`$ tff ($`t=2.47\times 10^5`$ years) and $`t=1.318`$ tff ($`t=2.51\times 10^5`$ years), respectively. They each contain around 1 M of gas when the star formation begins. They produce 7 and 5, respectively, objects during the calculation.

The evolution of the velocity field and the density structure in the cloud will be discussed in more detail in a later paper. In this paper, we concentrate on the process of star formation and the properties of the resulting stars and brown dwarfs. When the calculation is stopped, the cloud has produced 23 stars and 18 brown dwarfs. An additional 9 objects have substellar masses but are still accreting. Three of these have very low masses and accretion rates and therefore would probably end up with substellar masses if the calculation were continued. The other six already have masses near the stellar/substellar boundary and are therefore likely to become stars.

### 4.1 The star formation process in the dense cores

Snapshots of the process of star formation in the most massive core are shown in Figures 3 and 4. The star formation in the two low-mass cores is depicted in Figures 5 and 6. A true appreciation of how dynamic and chaotic the star formation process is can only be obtained by studying an animation of the simulation. The reader is encouraged to download an animation of the simulation from one of the two internet sites http://www.astro.ex.ac.uk/people/mbate/Research/Cluster/ or http://www.ukaff.ac.uk/starcluster .

The gravitational collapse of the most massive dense core produces filamentary structures which fragment (e.g. Bastein 1983; Bastien et al. 1991; Inutsuka & Miyama 1992) to form a combination of single objects and multiple systems (Figure 3, $`t=1.061.10`$ tff). Many of the multiple systems result from the fragmentation of massive circumstellar discs that form around single objects that fragment out of the filaments (e.g. Bonnell 1994; Bate & Bonnell 1994; Whitworth et al. 1995; Burkert, Bate & Bodenheimer 1997). Subsequently, most of these objects fall together into the gravitational potential well to form a small stellar cluster (Figure 3, $`t=1.121.14`$ tff). The cluster only contains $`13`$ objects and, thus, dynamical interactions quickly result in its dissolution (Figure 3, $`t=1.161.22`$ tff). The formation of stars is essentially halted during this dissolution phase because a significant fraction of the gas was used up in producing the stars and brown dwarfs. Later, when more gas has fallen into the potential well of the main dense core, a new burst of star formation occurs (Figures 3 and 4, $`t=1.261.34`$ tff) producing a second cluster that contains $`16`$ objects at any one time. Again, this small cluster quickly disperses due to dynamical interactions (Figure 4, $`t=1.341.40`$ tff), although at the end of the calculation it has not yet fully dissolved and still contains $`11`$ objects.

The second and third dense cores produce only a small number of objects compared to the most massive core (Table 1). The dynamical interactions in these two stellar groups work to arrange dynamically unstable systems into more stable configurations (Figures 5 and 6). During the calculation, neither of these low-mass cores ejects any objects, although ejections are likely in the long-term. The formation of the multiple systems in these low-mass cores results almost exclusively from the fragmentation of massive circumstellar discs (Figures 5 and 6). Only in the third core, near the end of the calculation, is an object formed in a separate filament (Figure 6, $`t=1.40`$ tff).

The formation of the stars and brown dwarfs is not a Poisson process. In Figure 7, along with the two overall bursts of star formation, it appears to be common for pairs of objects or several objects to be formed within a short time of each other. This can be investigated by considering the distribution of time intervals between successive star formation events. In Figure 8, we plot the cumulative distribution of these time intervals (solid line) and compare it with the distribution that would be expected if the events were Poissonian (dashed line). A Kolmogorov-Smirnov test of the two distributions gives only a 0.7% probability of the time intervals being distributed in a Poisson manner. Instead, there is an excess of short time intervals. This is due to multiple fragmentation events in gravitationally-unstable discs (e.g. Figure 5, $`t=1.331.34`$ tff), along with the fragmentation of collapsing gas filaments to form binaries (e.g. the first two stars to form, caption of Figure 3). Thus, in addition to the bursts of star formation on the scales of molecular cloud cores, fragmentation events on small-scales are correlated.

### 4.2 Comparison with observed star-forming regions

Comparing our calculation to real star-forming regions that have been well studied, this calculation could be regarded as modelling part of the $`\rho `$ Ophiuchus dark cloud. The main cloud of Ophiuchus contains around 550 M of gas within an area on the sky of around $`2\times 1`$ pc (Wilking & Lada 1983). It is centrally-condensed and contains 6 main dense cores. The gas masses of the cores range from $`862`$ M with mean densities of $`n(\mathrm{H}_2)10^510^6`$ cm-3 (Motte, Andre & Neri 1998). The calculation presented here contains a total of 50 M and forms one large and two smaller dense cores in a region $`0.4`$ pc across. Just before the star formation begins, the large core contains $`5.3`$ M and measures $`0.1`$ pc in diameter (counting gas with $`n(\mathrm{H}_2)>3\times 10^5`$ cm-3). Its mass, mean density, and size are similar to those of the Ophiuchus-F core which contains $`8`$ M (Motte et al. 1998). Thus, the calculation presented here is comparable to modelling a region containing one of the main dense cores within the $`\rho `$ Ophiuchus dark cloud.

Alternatively, the calculation can be viewed as modelling part of the cloud that formed the Orion Trapezium Cluster. Although the total mass of our cloud is far less than the progenitor of the Trapezium Cluster and our cloud produces only low-mass stars, it is similar in terms of the resulting stellar density. The stellar densities in the centre of the Trapezium Cluster are $`24\times 10^4`$ pc-3 (McCaughrean & Stauffer 1994; Hillenbrand & Hartmann 1998; Bate, Clarke & McCaughrean 1998) and the density falls off with radius as $`n2\times 10^4(r/0.07\mathrm{pc})^2`$ pc-3 (Bate et al. 1998). Our calculation produces $`30`$ stars from a cloud with an initial diameter of 0.375 pc giving a stellar density of $`1\times 10^3`$ pc-3. Thus, the overall stellar densities in our calculation are similar to those in the Trapezium Cluster 0.3 pc from the centre (a calculation of the radius at which the stellar densities fall to $`10^3`$ pc-3 using Hillenbrand & Hartmann’s King model fit with $`r_0=0.16`$ pc gives a radius of $`0.38`$ pc). This is within the Trapezium Cluster’s half-mass radius of 0.8 pc (Hillenbrand & Hartmann 1998). Thus, for processes that depend on stellar densities, such as the dynamical truncation of circumstellar discs, our calculation can also be compared with the Orion Trapezium Cluster.

## 5 Properties of the stars and brown dwarfs

For the remainder of the paper, we examine the properties of the stars and brown dwarfs that form in the calculation and compare them to the observed properties of stars and brown dwarfs. In this way, we aim to test whether our current understanding of the fragmentation of a turbulent molecular cloud is realistic and to predict some properties of stars and brown dwarfs that have not yet been determined observationally.

### 5.1 Star formation efficiency and timescale

As described above, the star formation occurs in one large and two small dense cores that form within the initial cloud due to the decay of the initial supersonic turbulence. The properties of the cores and the stars and brown dwarfs they produce are given in Table 1, along with the totals for the cloud as a whole. For each of the cores, the local star formation efficiency is high. In fact, it is this high efficiency that is responsible for the bursts of star formation in Figure 7. The most massive core undergoes a rapid burst of star formation (Figure 3, $`t=1.061.12`$ tff) lasting $`18000`$ years that severely depletes its reservoir of gas, temporarily halting the star formation. However, the core (including those stars and brown dwarfs that are not rapidly ejected) still dominates the local gravitational potential and, therefore, attracts more gas (Figure 3, $`t=1.141.24`$ tff). When the gas becomes sufficiently dense, a new burst of star formation occurs (Figure 3, $`t=1.261.34`$ tff) and the process is repeated.

Although the local star formation efficiency is high in the dense cores, most of the gas in the cloud is in low-density regions where no star formation occurs. Thus, the overall star formation efficiency of the cloud is low, $`10`$%. Due to computational limitations, we have not been able to follow the cloud until star formation ceases entirely. However, by the end of the calculation a large fraction of the initial cloud has drifted off to large distances due to a combination of the initial velocity dispersion and pressure gradients and is not gravitationally unstable. Thus, the global star formation efficiency is unlikely to rise above a few tens of percent. Furthermore, the calculation neglects all feedback processes. The most massive star to form is only 0.73 M so feedback may not be very important. Even so, jets, outflows and heating of the gas may be expected to reduce the efficiency of star formation still further. One way in which the star formation efficiency might be increased is if the initial turbulence in the cloud was reduced since then less of the gas would be able to drift away from the cloud.

Observations show that star formation efficiencies vary widely across star-forming regions. Efficiencies may be high locally, but are generally low globally. For example, in the $`\rho `$ Ophiuchus cloud, Wilking and Lada \[Wilking & Lada 1983\] find the overall star formation efficiency to be 20-30% but may be as high as 47% locally. On a grander scale, star formation in the Orion B molecular cloud does not occur uniformly but is concentrated in massive ($`M>200`$ M) dense cores that have high star formation efficiencies \[Lada 1992\]. Similarly, the Orion Trapezium Cluster has a high star formation efficiency \[Hillenbrand & Hartmann 1998\], but it is just one small part of the Orion A molecular cloud. We note that the star formation efficiency is much easier to measure from simulations than from observations because stars can quickly disperse from the dense cores in which they formed (see Section 5.2).

The timescale on which star formation occurs is the dynamical one. Star formation begins after one global free-fall time of the cloud. Similarly, the two bursts of star formation in the most massive core, each lasting for $`18000`$ years, and the time between the bursts (i.e. the time-scale for gas to fall into the core and replenish it) of $`24000`$ years are all roughly equal to the dynamical timescale of the massive core. This is consistent with the idea that star formation is a highly dynamical process (Pringle 1989; Elmegreen 2000; Hartmann, Ballesteros-Paredes & Bergin 2001). Although the calculation described here does not include magnetic fields, calculations by other authors including magnetic fields also find that the dissipation of turbulent support and gravitational collapse occur on the dynamical timescale of the clouds they model (e.g. Ostriker et al. 2001; Li et al. in preparation). Hence, the fact that we have neglected magnetic fields does not invalidate this conclusion.

### 5.2 Stellar velocity dispersion and distribution

The main dense core forms a small cluster of stars and brown dwarfs. This cluster never contains more than $`16`$ objects, so its timescale for dissolution is similar to that on which the objects are forming. The objects undergo chaotic interactions and eventually disperse with stars and brown dwarfs ejected in random directions from the core. The velocities of the stars and brown dwarfs relative to the centre-of-mass velocity of all the objects are given in Figure 9. The root mean square velocity dispersion is 2.1 km/s in three dimensions (3-D) or 1.2 km/s in one dimension (1-D) (using the centre-of-mass velocity for close binaries). This is roughly a factor of three greater than the 3-D velocity dispersion of the gas when the stars begin to form ($`=3.8`$ giving a 3-D velocity dispersion of 0.7 km/s). Thus, dynamical interactions are the primary source of the overall stellar velocity dispersion, but they do not often result in extreme ejection velocities.

The velocity dispersion of the stars and brown dwarfs is similar to that observed in low-mass star-forming regions. For example, the 1-D velocity dispersion in Taurus-Auriga has been measured at $`\text{ }<2`$ km/s using proper motions (Jones & Herbig 1979; Hartmann et al. 1991; Frink et al. 1997). The radial velocity dispersion of stars in Chamaeleon I has been measured at $`3.6`$ km/s \[Joergens & Guenther 2001\], but it is thought that this value is high due to the radial velocity ‘noise’ exhibited by T Tauri stars \[Guenther et al. 2001\]. For the Orion Trapezium Cluster, the 1-D velocity dispersion is $`2.3`$ km/s (Jones & Walker 1988; Tian et al. 1996). This is somewhat greater than in our calculation, but the stellar cluster that we form also has a much lower mass.

We find the velocity dispersion is independent both of stellar mass and of binarity (Figure 9). The lack of dependence on mass has been found in N-body simulations of the break up of small-N clusters with $`N>3`$ \[Sterzik & Durisen 1998\] and in SPH calculations of N=5 clusters embedded in gas (Delgado-Donate, Clarke & Bate 2003). Essentially, as $`N`$ increases, the escape speed is determined by the gravitational potential well of the group rather than by individual stellar masses. However, both Sterzik & Durisen and Delgado-Donate et al. found that binaries should have a smaller velocity dispersion than single objects due to the recoil velocities of binaries being lower. This is not what we find. Part of this is due to the fact that the stellar velocities in our calculation are contributed to by the motions of the dense cores (i.e. turbulent motions in the gas), something that was not considered in the above studies because they modelled only isolated small-N clusters. Another factor is that the most massive core in our calculation produces a large number of objects ($`N=38`$) and the presence of gas allows dissipative interactions between stars. Thus, the core forms several binaries, some of which are ejected, essentially eliminating any difference in the velocity dispersions of binaries and single objects. The formation of several binaries is difficult in N-body simulations (e.g. Sterzik & Durisen 1998) due to the absence of dissipative interactions between the stars and gas, while in $`N=5`$ clusters with gas (Delgado-Donate et al. 2003) it is unlikely simply due to the small number of objects.

Both Sterzik & Durisen \[Sterzik & Durisen 1999\] and Reipurth & Clarke \[Reipurth & Clarke 2001\] suggested that brown dwarfs might have large radial velocities due to their formation in and subsequent ejection from multiple systems. However, although all the brown dwarfs formed in our calculation do result from the ejection of low-mass objects from unstable multiple systems (see Section 5.4.1), our results show that this mechanism does not generally produce large ejection velocities. This is confirmed observationally by Joergens & Guenther \[Joergens & Guenther 2001\] who studied the radial velocities of brown dwarfs in the Chamaeleon I dark cloud and found a velocity dispersion of only 2 km/s. One potential area of concern is that we soften the gravitational fields of stars and brown dwarfs on scales less than 4 AU (Section 3.2). Thus, ejection velocities from very close encounters will be underestimated. However, we have plotted the final velocities of stars and brown dwarfs versus their closest encounter distance and find no correlation (either for objects with closest encounters greater than or less than 4 AU). The reason for this is probably the same as that which results in the independence of the velocity dispersion on mass, namely that the most likely ejection velocity from a stellar group depends on the escapee having sufficient kinetic energy to escape the gravitational potential well of the group rather than the details of individual interactions. Furthermore, the escape velocity of two brown dwarfs separated by 4 AU is $`7`$ km/s which is significantly higher than the velocity dispersion we find. Thus, we can safely conclude that our results are not significantly affected by our gravitational softening length.

Although we find that the 3-D velocity dispersion is quite low, the objects are still able to travel $`0.2`$ pc (the radius of the initial cloud) in $`10^5`$ years (the timescale over which star formation occurs). Therefore, at the end of the calculation, the cores are surrounded by a halo of objects. An observer looking at the positions of these stars on the sky would have no idea whether these stars formed in their current locations or in the dense cores. Only by determining their 3-D space velocities to high precision could our prediction that stars form in dense cores and are ejected from them be tested observationally.

Nevertheless, there is some observational support for the concept of dense cores resulting in small, expanding associations of stars. For example, in Taurus the 6 Gomez groups \[Gomez et al. 1993\] have radii of $`0.51.0`$ pc and contain $`1020`$ stars each. With 1-D expansion velocities of $`1`$ km/s they need only be $`0.51`$ million years old. A more recent survey of pre-main-sequence stars in Ophiuchus \[Allen et al. 2002\] finds that the dense cores are surrounded by young stars within one or two tenths of a parsec. It is intriguing to speculate that these stars formed within the observed dense cores and were ejected to occupy their current locations while more stars are currently being formed in the dense cores in new ‘bursts’ of star formation (as discussed in Section 5.1). The stars observed by Allen et al. could have travelled to their current locations in $`\text{ }<2\times 10^5`$ years.

### 5.3 Initial mass function

In Figure 10 we plot the initial mass function (IMF) obtained from the calculation. This is the first initial mass function to be determined from a hydrodynamical calculation that resolves objects down to (and beyond) the opacity limit for fragmentation. Hence, it predicts both the stellar and sub-stellar IMF. It is worth pointing out that the number of objects formed (50 stars and brown dwarfs) is larger than the number of objects observed in some surveys of star-forming regions.

We obtain a mass function that is consistent with

$$\frac{\mathrm{d}N}{\mathrm{dlog}M}M^\mathrm{\Gamma }$$
(2)
where

$$\mathrm{\Gamma }=\{\begin{array}{cc}\hfill 1.35& \text{for }M\text{ }>0.5\text{ M}\text{}\hfill \\ \hfill 0.0& \text{for }0.006<M\text{ }<0.5\text{ M}\text{}\hfill \end{array}$$
(3)
and there are no objects below the opacity limit for fragmentation ($`0.005`$ M). The Salpeter slope is $`\mathrm{\Gamma }=1.35`$ \[Salpeter 1955\].

In the present calculation, the effect of the opacity limit is mimicked by the switch from an isothermal equation of state to $`\eta =7/5`$ at a density of $`10^{13}`$ g cm-3 (equation 1). This is only an approximation of the true behaviour of the gas. In reality, the opacity limit depends on the heating and cooling rates of collapsing gas (Section 2). In order to get the exact value of the opacity limit for fragmentation correct, calculations including full radiative transfer are required. Thus, while the precise mass at which the cutoff in the IMF occurs cannot be derived from the present calculation, the qualitative result that there is a sharp cutoff in the mass function due to the opacity limit at $`0.01`$ M is expected.

#### 5.3.1 Comparison with the observed IMF

How does this theoretical IMF compare with the observed IMF? It is generally agreed that the observed IMF has a Salpeter-type slope above $`0.5`$ M and a flatter slope below this mass (e.g. Kroupa 2001; Luhman et al. 2000). Our theoretical IMF is consistent with these observations, although the relatively small number of stars and the fact that we don’t form any stars with masses $`\text{ }>0.75`$ M means that we cannot say much about the slope above 0.5 M. On the other hand, there is no a priori reason that the calculation could not have produced, say, one 10-M star (or even a 2-M star) and a cluster of brown dwarfs, something that would have been inconsistent with the observed IMF.

The real prediction of this calculation is for the low-mass stellar and substellar portions of the IMF where, unlike observations, we can detect all the brown dwarfs that form. This region of the IMF, particularly the substellar portion, is currently attracting an enormous amount of observational effort, and its exact form is still open to debate. Kroupa \[Kroupa 2001\] considered many observational surveys and concluded that the IMF rises slowly from 0.5 M to the stellar/substellar boundary ($`\mathrm{\Gamma }=0.3`$) and decreases with $`\mathrm{\Gamma }0.7`$ in the brown dwarf regime (Figure 10). The most recent surveys appear to be converging to a substellar slope of $`\mathrm{\Gamma }0.5`$ (Hillenbrand & Carpenter 2000; Najita, Tiede & Carr 2000; Moraux, Bouvier & Stauffer 2001). With our small number of objects, the best we can say at the moment is that our theoretical IMF is roughly flat in the range $`0.010.5`$ M and is consistent with these observational results.

A flat slope from $`0.010.5`$ M implies that there are roughly equal numbers of stars and brown dwarfs (Figure 10; Table 1). Observational surveys in the local solar neighbourhood agree with this prediction; both Reid et al. \[Reid et al. 1999\] and Chabrier \[Chabrier 2002\] estimate that the numbers of stars and brown dwarfs are roughly equal.

In Section 5.1, we likened our calculation to modelling part of the $`\rho `$ Ophiuchus star-forming region. Several papers have considered the low-mass IMF in $`\rho `$ Ophiuchus. Comeron et al. \[Comeron et al. 1993\] and Strom, Kepner & Strom \[Strom, Kepner & Strom 1995\] both estimated $`\mathrm{\Gamma }=0`$ in the range $`0.051`$ M from luminosity functions. Luhman & Rieke (1999) find a break in the IMF at $`0.5`$ M and estimate a lower limit of $`0.5`$ for the slope in the range $`0.020.5`$ M. Finally, the HST/NICMOS survey of Allen et al. \[Allen et al. 2002\] estimates that 30% of the objects in $`\rho `$ Ophiuchus may be brown dwarfs. Thus, our results are consistent with the IMF in Ophiuchus.

#### 5.3.2 The origin of the IMF

There is an ongoing discussion regarding the origin of the IMF (see the review by Kroupa 2002). Some suggest that it is due to stellar feedback (e.g. Silk 1977b; Adams & Fattuzzo 1996), others that it is due to the density structure in molecular gas (e.g. Larson 1978; Motte et al. 1998; Elmegreen 2002, Padoan & Nordlund 2002), and yet others attribute it to competitive accretion of gas and dynamical interactions (e.g. Larson 1978; Zinnecker 1982; Bonnell et al. 1997, 2001a,b; Klessen et al. 1998; Klessen & Burkert 2000).

The calculation presented here does not include feedback processes, yet still produces a realistic IMF. Moreover, there is no direct mapping between the instantaneous density structure in the gas and the masses of the stars that result. Although the stars form within dense cores, their individual masses are determined by accretion over the same time-scale on which the gas distribution in the core evolves. Furthermore, with each core containing multiple stars, there is no one-to-one link between core masses and final stellar masses. Rather, the origin of the IMF is best described as resulting from the combination of cloud structure, competitive accretion, and dynamical encounters that eject the stars from the dense cores and terminate their accretion (Bonnell et al. 1997, 2001a,b; Klessen et al. 1998; Klessen & Burkert 2000, 2001). It is a highly dynamical, chaotic process that needs to be understood in a statistical manner.

To determine how the initial mass function depends on environment will require further large-scale calculations in which the initial conditions such as the thermal Jeans mass and the density and velocity structure of the gas are varied. From the current calculation, we find that, if we take our ‘best guess’ values for the initial conditions in local star-forming molecular clouds, we do reproduce the observed IMF.

### 5.4 Brown dwarfs

Even though the initial thermal Jeans mass in the cloud is 1 M, the calculation produces many brown dwarfs with masses as low as $`0.005`$ M. Even if we take into account the initial velocity dispersion of the gas, with Mach number $`=6.4`$, the local Jeans mass in isothermal shocks is still $`0.16`$ M and by the time the stars and brown dwarfs begin to form the Mach number has dropped to $`=3.8`$. Therefore, it is of interest to determine how the brown dwarfs form. In particular, do they form in the same manner as the stars or in some other way?

#### 5.4.1 The formation of brown dwarfs

The formation mechanism and resulting properties of the brown dwarfs in this calculation have been studied in detail in the companion paper, Bate et al. \[Bate, Bonnell, & Bromm 2002a\]. We find that all 18 of the definite brown dwarfs (i.e. those that are not accreting significantly at the end of the calculation) form in dynamically-unstable multiple systems and are ejected from the regions of dense gas in which they form before they can accrete enough gas to become stars. We emphasise that all objects begin as opacity-limited fragments containing only a few Jupiter masses (Section 3.2). Those that subsequently become stars accrete large quantities of gas from the dense cores in which they form, while those that remain as brown dwarfs do not because they are ejected.

Watkins et al. \[Watkins et al. 1998\] speculated that brown dwarfs may be ejected from unstable multiple systems, but they did not recognise the importance of these ejections in halting accretion. Reipurth & Clarke \[Reipurth & Clarke 2001\] proposed that brown dwarfs form because dynamical ejections halt the accretion onto objects that would otherwise become stars. This is the mechanism by which all our brown dwarfs form. However, without numerical calculations, the main processes involved in the formation of these unstable multiple systems, the efficiency of the ejection mechanism, and the resulting properties of the brown dwarfs can only be conjectured.

We find (Bate et al. 2002a) that roughly three quarters of the brown dwarfs (14 of the 18) form the via the fragmentation of gravitationally-unstable circumstellar discs (Bonnell 1994; Whitworth et al. 1995; Burkert et al. 1997). The remaining brown dwarfs (4 of the 18) form via the fragmentation of filaments of molecular gas (e.g. Bonnell et al. 1991). These objects either form in, or quickly fall into, unstable multiple systems are subsequently ejected from the dense gas, limiting their masses to be substellar.

Examining the origins of the 23 stars formed in the calculation, we find that only one third (7 of the 23) form via disc fragmentation. The majority (16 of the 23) form directly from the collapse of the cloud in filaments of molecular gas. Thus, both stars and brown dwarfs can form by both mechanisms, but the primary mechanism by which stars form is the fragmentation of collapsing gas while the primary mechanism for brown dwarf formation is disc fragmentation. As proposed by Reipurth & Clarke \[Reipurth & Clarke 2001\], the main difference between stars and brown dwarfs is that the brown dwarfs are ejected soon after they are formed, before they have been able to accrete to stellar masses. This is more likely if a disc fragments into multiple objects than if an object forms on its own and has to collide with a nearby multiple system before it is ejected from the dense gas. Thus, the different fractions of stars and brown dwarfs that form from cloud and disc fragmentation are explained.

The ease with which disc fragmentation occurs depends primarily on the rate at which it accretes mass from the surrounding cloud (Bonnell 1994) and the disc’s equation of state (e.g. Pickett et al. 2000). The density of the gas in the discs that fragment to form brown dwarfs is high enough that the gas is in the $`\eta =7/5`$ regime (Section 2.1). Thus, the gas resists fragmentation far more than it would if an isothermal equation of state were used, although our equation of state does not include heating from shocks. On the other hand, because the flattened disc geometry may allow more rapid cooling than a spherically-symmetric geometry, real discs may be cooler and more unstable than those we model here. Hence, the number of brown dwarfs may be even greater. Another factor is that stars and brown dwarfs cannot merge in our calculation. Mergers may reduce the final number of brown dwarfs. In summary, a more definitive prediction will have to wait until a large-scale calculation is performed with radiative transfer and even higher resolution. For the present, we have demonstrated that the fragmentation of a turbulent molecular cloud is capable of forming similar numbers of brown dwarfs and stars as well as an IMF that is in agreement with observations (Section 5.3).

#### 5.4.2 The frequency of binary brown dwarfs

We find a low frequency ($`5`$%) of binary brown dwarf systems (Bate et al. 2002a). This is primarily due to the closeness of the dynamical encounters that eject the brown dwarfs from the dense gas in which they form before they can accrete to stellar masses. The minimum separations during the encounters are usually less than 20 AU (Section 5.6), so any wide systems are usually disrupted. However, another type of dynamical interaction also plays a role. Several binary brown dwarf systems that form during the calculation are destroyed by exchange interactions where one or both of the brown dwarfs are replaced by stars.

These two effects result in the low frequency of binary brown dwarfs. Of the 18 definite brown dwarfs, none are in binaries. However, at the end of the calculation there is a close binary brown dwarf (semimajor axis 6 AU) within an unstable multiple system (Figure 5, $`t=1.40t_{\mathrm{ff}}`$; Table 2, binary 39,41). This system will undergo further dynamical evolution, and it is still accreting. Since the binary brown dwarf is very close, it is possible that it will survive the dissolution of the multiple system and be ejected before it has become a stellar binary. Even so, this would result in only one binary brown dwarf system and $`20`$ single brown dwarfs. Thus, the formation of close binary brown dwarfs is possible, but the fraction of brown dwarfs with a brown dwarf companion should be low ($`5`$%).

Observationally, the frequency of brown dwarf binaries is not yet clear. Both Reid et al. (2001) and Close et al. (2002) observed 20 brown dwarf primaries and found that 4 have companions giving binary frequencies of $`20`$%. However, as discussed by Close et al., these surveys are magnitude limited rather than volume limited and therefore are likely to overestimate the true frequency of brown dwarf binaries. On the other hand, these surveys cannot detect faint companions or those at very small separations, which would boost the frequency. Our calculation favours a lower frequency, but due to our small number of objects, we cannot exclude a frequency of 20% (there would be a probability of $`6`$% of finding 1 binary out of 20 systems). It is important to note that none of the binary brown dwarf systems currently known have projected separations $`>15`$ AU (Reid et al. 2001; Close et al. 2002). This is consistent with their having survived dynamical ejection from unstable multiple systems.

#### 5.4.3 The sizes of discs around brown dwarfs

Another result of the close dynamical encounters that occur during the ejection of the brown dwarfs from the dense gas is that most brown dwarfs do not have large discs. We find that only one ($`5`$%) of the definite brown dwarfs is ejected with a large (radius $`\text{ }>20`$ AU) circumstellar disc. Of the 18 definite brown dwarfs, 14 have dynamical encounters at separations $`<23`$ AU (Section 5.6). However, this is not the only reason that few brown dwarfs have large discs. To become brown dwarfs they must, by definition, be ejected from the dense gas before they have been able to accrete much gas. Three of the definite brown dwarfs are ejected so soon after their formation that they have not had time to accrete gas with the high specific angular momentum required to form large discs. Only one of the 18 definite brown dwarfs is ejected with a resolved disc (radius $`60`$ AU). The other brown dwarfs will possess smaller discs, but in this calculation we are unable to resolve discs with radii $`\text{ }<10`$ AU (Section 3.2). In any case, truncation radii $`10`$ AU are not observationally relevant since the viscous timescale at this radius is

$$\begin{array}{cc}t_{\mathrm{visc}}\hfill & \frac{1}{\alpha }\left(\frac{R}{H}\right)^22\pi \sqrt{\frac{R^3}{GM}}\hfill \\ & \\ & =1\times 10^6\left(\frac{R}{10\mathrm{A}\mathrm{U}}\right)^{3/2}\left(\frac{0.075\mathrm{M}_{}}{M}\right)^{1/2}\mathrm{yr}\hfill \end{array}$$
(4)
where we have assumed the standard Shakura-Sunyaev (1973) viscosity parameter $`\alpha =0.01`$ (Hartmann et al. 1998) and the disc height to radius $`H/R=0.1`$ (e.g. Burrows et al. 1996; Stapelfeldt et al. 1998). Thus, any disc that is truncated to a radius $`10`$ AU around a brown dwarf is likely to evolve viscously to a radius of $`10`$ AU before it is observed. Of course, such a disc may have a very low mass, first because of the truncation and second because much of the remaining gas would be accreted during its viscous evolution back out to $`10`$ AU radius. Armitage & Clarke \[Armitage & Clarke 1997\] studied the evolution and observability of discs truncated to between 1 and 10 AU. They found that the K-band excess becomes too low to detect within a few times $`10^5`$ years, but that infrared emission at wavelengths $`\text{ }>5`$ $`\mu `$m is detectable for much longer.

Muench et al. \[Muench et al. 2001\] find $`65`$% of brown dwarfs in the Orion Trapezium cluster have infrared excesses indicative of discs. Spectral energy distributions indicative of discs have also been found around young brown dwarfs in low mass star-forming regions (Natta & Testi 2001; Testi et al. 2002) as have H$`\alpha `$ and Br$`\gamma `$ accretion signatures (Martin et al. 2001a; Zapatero Osorio et al. 2002b). Our results show that if brown dwarfs are formed by the ejection mechanism, few of these discs should have radii $`>20`$ AU and most discs should have radii $`10`$ AU, although those that were severely truncated and evolved back out to this size might be difficult to detect.

#### 5.4.4 Brown dwarf companions to stars

As we mention in Section 4, when the calculation is stopped, 9 objects have substellar masses but are still accreting. All of these potential brown dwarfs are members of unstable multiple systems (Table 2). Along with the one close binary brown dwarf system mentioned above, there is one close stellar/substellar binary (26,40 in Table 2). This system is part of the same unstable septuple system that contains the binary brown dwarf system. It has a semi-major axis of 7 AU. It is known that brown dwarf companions to FGK stars with semi-major axes $`\text{ }<4`$ AU are very rare ($`\text{ }<1`$%; Halbwachs et al. 2000; Zucker & Mazeh 2001; Pourbaix & Arenou 2001). This is the so-called brown dwarf desert. Our stellar/substellar system is of interest because it indicates that there may not be a brown dwarf desert for M stars. Of course, it has several hurdles to pass before it becomes a stable stellar/substellar system. It must avoid the accretion of much more gas, which would tend to equalise the masses \[Bate 2000\]; it must survive the break up of the septuple system; and the brown dwarf must avoid migrating into the star due to its interaction with the disc \[Armitage & Bonnell 2002\]. We note that the ejection hypothesis for the formation of brown dwarfs is not necessarily inconsistent with the detection of close brown dwarf companions to stars. The hypothesis relies on the assertion that brown dwarfs are ejected from the regions of dense gas in which they form before they are able to accrete to stellar masses. However, there is no reason why they cannot be in a close orbit around a star when the system is ejected. More difficult would be the formation of wide ($`\text{ }>100`$ AU) brown dwarf companions to stars by this mechanism since it is difficult to see how such a wide system could be ejected dynamically. Such systems are known to exist (Nakajima et al. 1995; Rebolo et al. 1998; Lowrance et al. 2000; Gizis et al. 2001; Els et al. 2001). In these cases, it may simply be that the specific angular momentum of the gas accreted by the system was low enough to be preferentially accreted by the primary, leaving the secondary with a substellar mass (Bate & Bonnell 1997; Bate 2000).

### 5.5 The formation of binary and multiple systems

The favoured mechanism for binary star formation is fragmentation, either of collapsing molecular cloud cores (Boss & Bodenheimer 1979; Boss 1986; Bonnell et al. 1991; Nelson & Papaloizou 1993; Burkert & Bodenheimer 1993; Bate, Bonnell & Price 1995; Truelove et al. 1998) or of massive circumstellar discs (Bonnell 1994; Whitworth et al. 1995; Burkert et al. 1997). However, it has been pointed out that, especially in small groups of stars, star-disc encounters may form binaries (Larson 1990; Clarke & Pringle 1991a,b; Heller 1995; McDonald & Clarke 1995). In a star-disc encounter resulting in capture, one star passes through the disc of another, dissipating enough kinetic energy to form a bound system \[Hall, Clarke & Pringle 1996\]. The calculation presented here can be used to examine which of these formation mechanisms is most prevalent.

We find that the dominant mechanism for the formation of the binary and multiple systems in the calculation is fragmentation, either of gaseous filaments or of massive circumstellar discs. Although star-disc encounters do occur during the calculation, most of these serve only to truncate the circumstellar discs (see Section 5.6) and do not result in bound stellar systems (c.f. Clarke & Pringle 1991a). There are a few exceptions. One star-disc encounter between a binary and a single star results in the formation of an unstable triple system which subsequently breaks up into a binary and a single star. Another star-disc capture is depicted in Figure 11. It is also important to note that, although star-disc encounters do not usually form simple bound systems directly, they do result in dissipation, which aids in the formation of the small-N bound groups that later dissolve and produce binary and multiple systems. Thus, dissipative encounters do play an important role in the star formation process, if not in the simple picture of star-disc capture.

#### 5.5.1 Multiplicity

At the end of the calculation, there exist 4 multiple systems or stellar groups. Their properties are displayed in Table 2 and Figure 12. Two of these systems originate in the first dense core. They are a close binary system that was ejected from the cloud (7,8) and the remains of a small-N group consisting of 11 objects (see also the last panel of Figure 4). Each of the other two dense cores contains a multiple system, one an unstable quadruple system and the other an unstable system consisting of seven objects (see the last panels of Figures 5 and 6).

All these systems, except the ejected close binary, will undergo further evolution. It is likely that most of the close binary systems and some of the triple systems will survive, but it is not possible to determine the final binary and multiple frequencies or properties of the wide systems that would eventually be formed from the simulation if it were continued. The best we can do is provide an upper limit on the final companion star frequency

$$CSF=\frac{B+2T+3Q+\mathrm{}.}{S+B+T+Q+\mathrm{}.}$$
(5)
where $`S`$ is the number of single stars, $`B`$ is the number of binaries, $`T`$ is the number of triples, etc. We have 26 single objects, 1 binary, 1 quadruple, 1 septuple and 1 system of 11 objects. This gives a companion star frequency of $`20/30=67`$%. Alternately, the number of companions divided by the total number of objects is $`20/50=40`$%. These high frequencies are in broad agreement with the large fractions of binary and multiple systems found in young star-forming regions (Ghez, Neugebauer & Matthews 1993; Leinert et al. 1993; Richichi et al. 1994; Simon et al. 1995; Ghez et al. 1997; Duchêne 1999). One potential difficulty with the calculation is that there are no wide binary systems when it is stopped. The multiple systems are expected to evolve into stable configurations that are most likely to be binaries or triples. However, it is not clear that these systems will be wide. Furthermore, although there are many low-mass objects in the multiple systems, these are the most likely components to be ejected in subsequent dynamical interactions. Thus, any wide binaries that do form may not have low-mass secondaries, whereas we know from observations that most wide binaries do have unequal mass components \[Duquennoy & Mayor 1991\]. Similar results are obtained by Delgado-Donate et al. (2003) in their simulations of $`N=5`$ clusters embedded in molecular cloud cores. These issues need to be addressed in future calculations.

#### 5.5.2 The formation of close binary systems

As reviewed in section 2, the opacity limit for fragmentation sets a minimum initial binary separation of $`10`$ AU. However, at the end of the calculation, there exist 7 close binary systems (separations $`<10`$ AU). The mechanisms by which these close binaries form and their properties have been discussed in detail in the second companion paper to this, Bate et al. (2002b). We only briefly summarise the conclusions here.

Bate et al. (2002b) analyse the fragmentation that occurs in the hydrodynamical calculation and find that the three smallest separations between an existing object and a forming fragment are 9, 21, and 22 AU (Figure 13, open circles). These separations are consistent with the lower limit set by the opacity limit for fragmentation. Only the last of these ends up in a close binary system; the other systems are disrupted by dynamical encounters. Bate et al. (2002b) find that, rather than forming directly via fragmentation, the close binary systems form through a combination of accretion, the interaction of binaries and triples with circumbinary and circumtriple discs, and dynamical interactions. Accretion onto a binary from an envelope decreases the binary’s separation unless the specific angular momentum of the accreted material is significantly greater than that of the binary (Artymowicz 1983; Bate 1997; Bate & Bonnell 1997; Bate 2000). It can also change the relative separations of a triple system, destabilising it and forcing dynamical interactions (Smith, Bonnell & Bate 1997). Circumbinary discs can remove large amounts of orbital angular momentum from an embedded binary system via gravitational torques, thus tightening its orbit (Pringle 1991; Artymowicz et al. 1991). However, although both of these processes play a role, the most important ingredient for the formation of the close binaries in our calculation is stellar dynamical interaction (Bate et al. 2002b). Dynamical interactions can harden existing wide binaries by removing angular momentum and energy from their orbits. They also produce exchange interactions in which a temporary unstable multiple system decays by ejecting one of the components of the original binary. Usually the lowest-mass object is ejected. Although such dynamical interactions play the dominant role in forming the close binaries, they cannot produce the observed close binary frequency on their own. Kroupa & Burkert \[Kroupa & Burkert 2001\] find that N-body calculations which begin with star clusters (100 to 1000 stars) consisting entirely of binaries with periods $`4.5<\mathrm{log}(P/\mathrm{days})<5.5`$ produce almost no binaries with periods $`\mathrm{log}(P/\mathrm{days})<4`$. Similarly, the dissolution of small-N clusters typically results in binaries with separations only an order of magnitude smaller than the size of the initial cluster (Sterzik & Durisen 1998). The main difference between these calculations and ours is that the dynamical interactions in our calculation are usually dissipative. Along with the effects of accretion and the interaction of multiple systems with circumstellar discs described above, the presence of gas allows dynamical encounters to be dissipative and to transport angular momentum. Such dissipative encounters include star-disc encounters (Larson 1990; Clarke & Pringle 1991a,b; McDonald & Clarke 1995; Heller 1995; Hall, Clarke & Pringle 1996) and other tidal interactions (Larson 2002).

These processes of accretion, interaction with circumbinary and circumtriple discs and dissipative dynamical interactions result in the formation of seven close binaries by the time the calculation is stopped. If all of these survive the break up of the remaining multiple systems, the resulting frequency of close binaries would be $`7/4316`$%. This is in good agreement with the observed frequency of close (separation $`<10`$ AU) binaries of $`20`$% \[Duquennoy & Mayor 1991\], demonstrating that close binaries need not be created by fragmentation in situ.

#### 5.5.3 The properties of close binary systems

The formation mechanisms discussed above lead to several consequences for the properties of close binaries (Bate et al. 2002b). We find a preference for equal masses (Table 2, Figure 12), with all close binaries having mass ratios $`q\text{ }>0.3`$ and most having $`q>1/2`$. This is due to the mass-equalising effect of the accretion of gas with high specific angular momentum (Artymowicz 1983; Bonnell & Bastien 1992; Whitworth et al. 1995; Bate & Bonnell 1997; Bate 2000) and dynamical exchange interactions that usually result in the ejection of the least massive component. These processes give a natural explanation for the observation that close binaries (periods $`\text{ }<10`$ years) tend to have higher mass ratios than wider binaries (Mazeh et al. 1992; Halbwachs, Mayor & Udry 1998; Tokovinin 2000). The wider multiple systems in the calculation have a larger range of mass ratios (Figure 12), although it is unclear how this would change with further evolution.

Successive dynamical exchanges also lead to a dependence of the close binary fraction on primary mass, since each time a binary encounters a star more massive than the primary, the most massive star will usually become the new primary. Of the $`20`$ brown dwarfs there is only one close binary brown dwarf system (Section 5.4.2), whereas 5 of the 11 stars with masses $`>0.2`$ M are members of close binary systems. While it is difficult to extrapolate these results to larger star clusters and more massive stars, this trend of the frequency of close binaries increasing with stellar mass is supported by observational surveys (Garmany et al. 1980; Abt et al. 1990; Morrell & Levato 1991; Mason et al. 1998).

At the end of the calculation, only one of the close binaries has been ejected from the cloud and is on its own. The others remain members of larger-scale bound groups and three are members of hierarchical triple systems (Table 2). This large number of wider companions is yet another indication of the importance of multiple systems in producing close binaries. Even allowing for the eventual dissolution of the bound groups, it seems likely that some of the hierarchical triple systems will survive. Although the true frequency of wide companions to close binaries is not yet well known, many close binaries do have wider components (e.g. Mayor & Mazeh 1987; Tokovinin 1997, 2000). Indeed, it was this observation that led Tokovinin (1997) to propose that dynamical interactions in multiple systems may play an important role in the formation of close binary systems. Further surveys to determine the true frequency of wide companions to close binary systems would be invaluable.

### 5.6 Protoplanetary discs

#### 5.6.1 The formation and evolution of discs

The calculation resolves gaseous discs with radii greater than $`10`$ AU around the young stars and brown dwarfs (see Section 3.2). We find that discs with typical radii of $`50`$ AU form around many of the stars and brown dwarfs soon after their creation due to the infall of gas with high specific angular momentum. Early on, these discs contain a large fraction of the mass of the star/disc system and are gravitationally unstable \[Lin & Pringle 1990\]. Spiral density waves form and lead to the rapid transport of angular momentum outwards and mass inwards via gravitational torques (e.g. Figures 5 and 6; Bate et al. 2002b). This serves to increase the mass of the central object(s) much faster than for a low-mass disc whose evolution is expected to be driven by the magneto-rotational instability (Balbus & Hawley 1991; Hawley & Balbus 1991). Even so, in many cases, the angular momentum transport via gravitational torques is insufficient to cope with the rapid infall of gas onto the disc \[Bonnell 1994\] and the disc fragments to form additional stars and brown dwarfs (Figures 5 and 6; Sections 5.4 and 5.5; Bate et al. 2002a).

Although many objects form with large discs initially, star-disc encounters truncate the majority of these discs. In Figure 14, we plot the closest encounter distance for each object during the calculation as a function of its final mass. Most of the stars and brown dwarfs have had encounters closer than 4 AU (recall that sink particle interactions are softened within 4 AU so minimum separations less than this are unreliable). During such an encounter, any large disc is truncated so that its radius is $`1/3`$ of the minimum separation during the encounter \[Hall 1997\]. However, even if an object has a close encounter, it may still end up with a large disc. New discs can form from gas accreted from the cloud after the encounter if the object has not been ejected from the cloud. In some cases, objects repeatedly have their discs stripped away and then replenished by new infall. Evidence for replenished discs can be seen in Figure 14, where objects are circled if they have resolved discs at the time of their ejection from the cloud, or at the end of the calculation. Many of these objects have had encounters closer than the resolution limit for discs.

The formation, truncation and replenishment of discs makes the issue of whether a star or brown dwarf ultimately has a resolved disc very complicated. In Table 3, we consider those 18 brown dwarfs that have stopped accreting by the end of the calculation and all of the stars (regardless of whether or not they are still evolving, since most stars are still evolving when the calculation is stopped). We find that the median closest encounter distance for a brown dwarf is $`10`$ AU and only 1 of the 18 brown dwarfs is ejected with a resolved disc (radius $`60`$ AU). The other brown dwarfs will possess smaller discs, but in this calculation we are unable to resolve discs with radii $`\text{ }<10`$ AU (Section 2.2). By contrast, most of the stars have had encounters closer than 1 AU, but more than half of these are surrounded by resolved discs at the end of the calculation (many of these are discs around binary and multiple systems).

The different results for the stars and brown dwarfs are due to the fact that the brown dwarfs must be ejected from the cloud before they accrete enough material to become stars. Three effects are involved. First, for the brown dwarf to be ejected before it becomes a star, it must be ejected soon after it forms. Three of the 18 definite brown dwarfs are ejected so quickly that they do not have time to accrete the high angular momentum gas required to form a large disc. Second, the encounters that result in ejection are usually close. Of the 18 definite brown dwarfs, 14 have undergone dynamical encounters with separations less than 23 AU, destroying any previously resolved discs. Finally, by definition, the last encounter a brown dwarf has is the one that ejects it from the cloud. After this, the disc cannot be replenished by new accretion from the cloud.

The stars, on the other hand, can remain in the dense gas and accrete for much longer. During this time, they may undergo close encounters and have their discs stripped away but if they remain in the cloud there is plenty of time for them to accumulate new discs. Thus, as seen in Table 3, the stars tend to have even closer encounters than the brown dwarfs (because they tend to remain in stellar groups longer), but most still have resolved discs when the calculation is stopped.

In summary, the issue of disc formation is complex. Brown dwarfs must quickly be ejected from the cloud, in order to remain as brown dwarfs. Thus, they tend to have their original discs, which are often severely truncated during the ejection process. Stars, by definition, are those objects that remain in the cloud long enough to accrete sufficient gas. Therefore, they have more time in which to undergo close encounters. They typically have smaller minimum encounter distances, but they also have plenty of time to replenish their discs and to accumulate large discs. The chaotic nature of disc formation, truncation, and replenishment naturally leads to a broad range of disc radii and masses. Thus, it should not be surprising that the discs around some objects in young clusters are undetectable or have already dispersed, while other objects with the same apparent age still have discs (e.g. Strom et al. 1989; Beckwith et al. 1990; Hillenbrand et al. 1998; Haisch et al. 2001a,b; Rebull et al. 2002; Armitage, Clarke & Palla 2002).

#### 5.6.2 Disc frequencies and sizes

At the end of the calculation, there are 11 resolved discs. Their radii and the objects they encircle are listed in Table 4. Six of the discs surround multiple systems, two of the discs encircle single objects that formed shortly before the calculation was stopped, and the remaining three surround single objects: one star and one brown dwarf, each of which were ejected during the calculation, and another substellar object that appears to be in the process of being ejected.

Because many of these discs reside in unstable multiple systems when the calculation is stopped, we are limited in the conclusions we can draw about final disc properties. However, as discussed in Section 5.4.3, of the 18 brown dwarfs that have been ejected and have finished accreting, only one has a resolved disc. Another object seems to be in the process of being ejected with a resolved disc (Table 4, object 47) at the end of the calculation. Recall from Section 5.4.3 that discs around brown dwarfs which are truncated to $`10`$ AU in radius would be expected to evolve viscously to a radius $`10`$ by the time they are observed. Thus, the vast majority of discs should have radii $`10`$ AU and only a small fraction ($`510`$ %) of brown dwarfs should have discs significantly larger than this (i.e. $`>20`$ AU).

Nine stellar systems (8 single stars and one close binary) have also dispersed from the cloud. Only one of the single stars has a resolved disc. Thus, these objects also have a low frequency of large discs $`10`$%. However, of the remaining 8 discs, 5 surround systems with stellar primaries and another surrounds a forming object so isolated that it will almost certainly become a star. The two remaining discs surround the close binary brown dwarf system (39,41 in Table 2) and a newly formed substellar object whose fate is unclear. Neglecting any further evolution of the unstable systems, we can place upper limits on the frequencies of discs with radii $`>20`$ AU around brown dwarf systems of $`3/2015`$% and around stars of $`7/1547`$%. Thus, given even these highly optimistic estimates, most of the stars and brown dwarfs do not retain discs large enough to form our solar system.

Is this size distribution of discs realistic? In a recent study of disc lifetimes in the Taurus star-forming region, Armitage et al. (2002) find that around 30% of stars lose their discs within 1 Myr while the remainder have disc lifetimes in the $`110`$ Myr range. The latter range is consistent with the lifetimes of discs that have an initial dispersion of half an order of magnitude in mass, but the very short disc lifetimes require another explanation. The large fraction of objects that suffer severe disc truncation in the calculation presented here may explain these very short disc lifetimes.

The only star-forming region for which we currently have direct information on the size distribution of circumstellar discs is the Orion Trapezium Cluster. This is because, unlike most star forming regions, a reflection nebula behind many of the young stars enables us to observe discs in silhouette (O’Dell & Wen 1994; McCaughrean & O’Dell 1996; O’Dell 2001). As discussed in Section 4.2, although our simulation only forms a small stellar association rather than a cluster on the scale of the Trapezium Cluster, the resulting stellar densities are comparable. Thus, the types of dynamical encounters that occur in our simulation may also have occurred during the formation of the Trapezium Cluster.

We know from infrared excesses that most of the stars and brown dwarfs in the Trapezium Cluster have discs (Hillenbrand et al. 1998; Lada et al. 2000). Lada et al. find that $`80`$% of the stars in the cluster have excesses indicative of discs, while for brown dwarfs the figure is $`65`$% \[Muench et al. 2001\]. Furthermore, of $`350`$ objects observed by the HST, $`150`$ are observed to have discs or proplyds indicating the presence of discs (O’Dell & Wen 1994; Stauffer et al. 1994; O’Dell & Wong 1996; O’Dell 2001). However, only $`40`$ of these exhibit resolved silhouette or ionised discs. HST can resolve discs and proplyds in Orion down to radii of $`40`$ AU. Thus, although the vast majority of stars in Orion have discs, the frequency of discs with radii greater than $`40`$ AU is only $`10`$% (Rodmann 2002; McCaughrean & Rodmann, in preparation). These numbers are broadly consistent with our results. The implication is that most stars in dense star-forming environments are incapable of forming large planetary systems like our own.

Of course, the Orion Trapezium is a high-mass star forming region and, as demonstrated by the proplyds, disc destruction by the O stars may influence the disc sizes along with dynamical encounters. Furthermore, the fact that the stellar densities in Orion are similar to those obtained in the calculation presented here does not guarantee that the encounters will be similar. Currently in Orion, there is a very low probability of an encounter at $`<100`$ AU (Scally & Clarke 2001). For encounters such as those discussed here to have occured in the Trapezium cluster, the stars must have formed in small groups which dispersed to form the Trapezium cluster as it appears today (i.e. initially it must have contained extensive substructure). Small-scale substructure would have disappeared by the current age of the cluster (Scally & Clarke 2002). Therefore, to test the prediction that most discs are small due to encounters properly, it is necessary not only to increase the resolution to observe smaller silhouette discs in Orion, but also to determine the disc size distribution in lower mass star-forming regions such as $`\rho `$ Ophiuchus. This should be possible with the Atacama Large Millimeter Array (ALMA) or the Submillimetre Array (SMA).

If the severe truncation of discs that we find is ruled out by future observations, we will be forced to find ways to reduce the number of dynamical truncations. One possibility is that the equation of state that we use results in discs that fragment too easily. Not only would this overestimate the number of brown dwarfs (Section 5.4.1), but it would decrease the number of large discs for two reasons. First, some discs that should survive would break up into fragments. Second, the increased number of objects would lead to more star-disc encounters and disc truncations. This possible dependence of the statistical properties on the equation of state should be investigated in future calculations.

## 6 Conclusions

We have presented the results from one of the most complex hydrodynamical star formation calculations to date. The calculation follows the collapse and fragmentation of a large-scale turbulent molecular cloud to form a stellar cluster consisting of 50 stars and brown dwarfs. The opacity limit for fragmentation is mimicked by the use of a non-isothermal equation of state and the resolution is sufficient to resolve all fragmentation before this physical limit is reached. Binary stars with separations as small as 1 AU and circumstellar discs with radii down to $`10`$ AU are resolved. The calculation allows us to study the formation mechanisms of stars and brown dwarfs and to determine a wide range of statistical properties for comparison with observation.

We find that star formation is a highly dynamic and chaotic process. A true appreciation of the process can only be obtained from examining an animation of the calculation. These can be downloaded from http://www.astro.ex.ac.uk/people/mbate or http://www.ukaff.ac.uk/starcluster. Fragmentation occurs both in dense molecular cloud cores and in massive circumstellar discs. Star-disc encounters form binaries and truncate discs. Stellar encounters disrupt bound multiple systems. The star formation occurs on the dynamical timescale of the cloud and the star formation efficiency across the cloud is variable with a low global efficiency ($`12`$% when we stop the calculation) but local efficiencies as high as $`50`$%. The high local efficiencies in dense molecular cores result in bursts of star formation because the rapid conversion of gas into stars depletes the high-density gas to such an extent that the star formation essentially comes to a halt until more gas has fallen into the core. When enough new gas has accumulated, another burst of star formation occurs.

We find that the opacity limit for fragmentation sets an initial mass for all fragments of $`0.005`$ M. Subsequently, the fragments accrete from the surrounding gas. Those that manage to accrete enough mass become stars ($`M\text{ }>0.075`$ M), while the rest are left as brown dwarfs. We propose that the initial mass function results from this accretion process, with each fragment accreting according to the conditions in which it is formed. The calculation produces a mass function that is consistent with a Salpeter slope ($`\mathrm{\Gamma }=1.35`$) above 0.5 M, a roughly flat distribution ($`\mathrm{\Gamma }=0`$) in the range $`0.0060.5`$ M, and a sharp cutoff below $`0.005`$ M. This is consistent with observational surveys.

Those objects that end up as brown dwarfs stop accreting before they reach stellar masses because they are ejected from the dense gas soon after their formation by dynamical interactions in unstable multiple systems. Thus, they can be viewed as ‘failed stars’. This ejection mechanism is very efficient, producing roughly equal numbers of stars and brown dwarfs (see also Bate et al. 2002a). However, the close interactions that occur during these dynamical ejections results in a low frequencies ($`5`$%) of binary brown dwarf systems. Similarly, the fraction of brown dwarfs with large (radii $`\text{ }>20`$ AU) circumstellar discs is $`5`$%. The accuracy of these frequencies is limited by our small number statistics (for example, we can only exclude a binary brown dwarf frequency of 20% at the 94% confidence level). However, further simulations will increase the significance of the predictions. Therefore, observational surveys to determine accurately the frequencies of binary brown dwarfs and the sizes of discs around brown dwarfs should be performed now so that we can test the models.

The calculation produces several binary and higher-order multiple systems. The opacity limit for fragmentation results in an initial minimum binary separation of $`10`$ AU. Despite this, we find that 7 close binary systems (separations $`<10`$ AU) exist when the calculation is stopped. These systems are produced by the hardening of initially wider multiple systems through a combination of dynamical encounters, gas accretion, and/or the interaction with circumbinary and circumtriple discs (see also Bate et al. 2002b). These mechanisms lead to close binaries having a bias towards equal-mass systems and a higher frequency of close binaries for higher-mass stars. Many of the close binaries also have wider companions. The resulting frequency of close binary systems is $`16`$%, consistent with observations. Thus, close binary systems need not be formed by fragmentation in situ.

Perhaps the most surprising result of this calculation is that most of the circumstellar discs in the calculation are severely truncated by dynamical encounters. Most young brown dwarfs should have discs with radii of $`10`$ AU, with large discs ($`\text{ }>20`$ AU) occurring around only $`5`$% of brown dwarfs. The discs around many stars are also severely truncated with the majority having radii $`\text{ }<20`$ AU (i.e. too small to form our solar system). Such severe disc truncation, and the associated low masses, may explain the observation that approximately 1/3 of the young stars in Taurus have lost their discs when they are only $`1`$ Myr old (Armitage et al. 2002). Currently, the only star-forming region in which we have information on the size distribution of circumstellar discs is the Orion Trapezium Cluster thanks to the silhouette discs. Our results are consistent with the sizes of discs in the Trapezium Cluster. However, massive stars are known to be evaporating discs in the Trapezium Cluster and the cluster is much larger than the system we are able to model. Thus, we strongly encourage observations to determine the size distribution of discs in low-mass star-forming regions such as $`\rho `$ Ophiuchus.

This is the first of a new generation of star-formation calculations that resolves all fragmentation and allows us to compare a wide range of statistical properties of stars and brown dwarfs with observations. Future calculations will determine the dependence of these properties on the initial Jeans mass in the cloud, the properties of the turbulence, and will improve the statistical significance of the results. In this way, we hope to understand better the origin of stars and brown dwarfs.

## Acknowledgments

We thank Mark McCaughrean, Cathie Clarke, and Eduardo Delgado-Donate for helpful discussions, and the referee, Ant Whitworth, for useful comments. The computations reported here were performed using the U.K. Astrophysical Fluids Facility (UKAFF).

## Appendix A Resolving the Jeans mass

Bate & Burkert \[Bate & Burkert 1997\] performed calculations of the isothermal collapse of a Jeans-unstable spherical molecular cloud core in solid-body rotation with a density distribution

$$\rho =\rho _0[1+0.1\mathrm{cos}(m\varphi )]$$
(6)
with $`m=2`$ and where $`\varphi `$ is the azimuthal angle about the rotation ($`z`$) axis. The ratios of the thermal and rotational energies to the magnitude of the gravitational potential energy of the cloud were $`\alpha =0.26`$ and $`\beta =0.16`$, respectively. The calculations were performed with a second-order finite difference hydrodynamics code (Burkert & Bodenheimer 1993, 1996) and a smoothed particle hydrodynamics (SPH) code (Bate et al. 1995).

They performed a high-resolution calculation with the grid code and a series of calculations with different resolutions using the SPH code. The initial $`m=2`$ density perturbation provides two overdense regions. As the cloud collapses, these overdense regions merge into an elongated structure. The two ends of this elongated structure each contain more than a Jeans mass so that, as the structure collapses, it fragments into a binary. Gas falls into the region between the binary, forming a dense gaseous filament. Thus, the result of the collapse is that the cloud fragments into a binary separated by a filament of gas. The subsequent evolution of the filament depends on the thermal behaviour of the gas \[Bate & Burkert 1997\].

Bate and Burkert obtained good agreement between the grid code and the SPH code provided sufficient resolution was used for the SPH calculations ($`2\times 10^4`$ SPH particles). With too few particles ($`1\times 10^4`$), the collapse of the elongated structure was delayed and when it finally did collapse it formed a filament without fragments at each end. The fragmentation was incorrectly modelled because the Jeans mass at each end of the elongated structure was not sufficiently resolved. Bate and Burkert derived an empirical resolution criterion from these calculations that the minimum Jeans mass during a calculation should be resolved by no fewer than twice the number of neighbouring particles over which SPH quantities are smoothed (i.e. $`2N_{\mathrm{neigh}}`$).

In the calculation presented in this paper, we wish to satistfy this Jeans mass criterion, but also to model the most massive molecular cloud possible. Bate and Burkert’s calculations only pin the resolution requirement down to between $`N_{\mathrm{neigh}}`$ and $`2N_{\mathrm{neigh}}`$ particles to resolve the local Jeans mass. A factor of two increase in the number of particles results in an increase of $`N^{4/3}=2.5`$ in computational time. For a calculation that requires $`10^5`$ CPU hours, this can make the difference between being able to perform the calculation and its being impractical.

To obtain a better estimate of the minimum number of particles required to resolve a Jeans mass sufficiently, we performed an identical SPH calculation to those presented in Bate & Burkert \[Bate & Burkert 1997\], but using 15000 particles. The results of this calculation are presented in Figure 15 in an identical manner to the results in Bate & Burkert \[Bate & Burkert 1997\]. It can be seen that 15000 particles is also sufficient to resolve the fragmentation and, thus, the resolution criterion can be further refined to $`1.5N_{\mathrm{neigh}}`$ particles being sufficient to resolve the minimum Jeans mass. Following Bate and Burkert, the maximum density for which the local Jeans mass can be resolved in the test calculation is then given by

$$\rho _{\mathrm{res}}\left(\frac{3}{4\pi }\right)\left(\frac{5R_\mathrm{g}T}{2G\mu }\right)^3\left(\frac{N_{\mathrm{tot}}}{1.5N_{\mathrm{neigh}}}\frac{1}{M_{\mathrm{tot}}}\right)^2$$
(7)
where $`R_\mathrm{g}`$ is the gas constant, $`T`$ is the temperature of the gas, $`G`$ is the gravitational constant, $`\mu `$ is the mean molecular weight, $`N_{\mathrm{tot}}`$ is the total number of particles, and $`M_{\mathrm{tot}}`$ is the total mass of the cloud. The density contour given by the thick line in Figure 15 gives this critical density. As expected, the fragmentation to form the binary occurs just as this critical density is surpassed.

Finally, we note that the above test calculation is purely isothermal, so the Jeans mass decreases monotonically with increasing density. Thus, the resolution criterion is derived from a calculation where a Jeans-unstable clump at the resolution limit of the calculation must continue to collapse and become more and more pronounced. For the calculation discussed in the main text of this paper, the minimum Jeans mass occurs at the density where the equation of state becomes non-isothermal (equation 1). A clump that is marginally Jeans unstable at this density cannot collapse because, as soon as the gas is compressed, it heats up and no longer contains a Jeans mass. Therefore, any clump that fragments with this equation of state must contain more mass than twice the minimum Jeans mass of 0.0011 M (1.1 MJ) and, thus, more than $`3N_{\mathrm{neigh}}`$ SPH particles.