On the evolution of a star cluster and its multiple stellar systems following gas dispersal

Abstract

We investigate the evolution, following gas dispersal, of a star cluster produced from a hydrodynamical calculation of the collapse and fragmentation of a turbulent molecular cloud. We find that when the gas, initially comprising ≈60 per cent of the mass, is removed, the system settles into a bound cluster containing ≈30–40 per cent of the stellar mass surrounding by an expanding halo of ejected stars. The bound cluster expands from an initial radius of <0.05 to 1–2 pc over ≈4–10 Myr, depending on how quickly the gas is removed, implying that stellar clusters may begin with far higher stellar densities than usually assumed. With rapid gas dispersal, the most massive stars are found to be mass segregated for the first ∼1 Myr of evolution, but classical mass segregation only develops for cases with long gas removal time-scales. Eventually, many of the most massive stars are expelled from the bound cluster. Despite the high initial stellar density and the extensive dynamical evolution of the system, we find that the stellar multiplicity is almost constant during the 10 Myr of evolution. This is because the primordial multiple systems are formed in a clustered environment and, thus, by their nature are already resistant to further evolution. The majority of multiple system evolution is confined to the decay of high-order systems (particularly quadruple systems) and the formation of a significant population of very wide (10⁴–10⁵ au) multiple systems in the expanding halo. This formation mechanism for wide binaries potentially solves the problem of how most stars apparently form in clusters and yet a substantial population of wide binaries exist in the field. We also find that many of these wide binaries and the binaries produced by the decay of high-order multiple systems have unequal mass components, potentially solving the problem that hydrodynamical simulations of star formation are found to underproduce unequal mass solar-type binaries.

methods: N-body simulations, binaries: general, stars: formation, stars: low-mass, brown dwarfs

1 INTRODUCTION

For many decades, it has been possible to perform N-body simulations of the evolution of stellar clusters, examining their dynamical evolution and the evolution of their stellar populations (von Hoerner 1960; van Albada 1968; Aarseth & Woolf 1972). Some early work was even done on the dynamical evolution of young clusters. For example, Aarseth & Hills (1972) studied the dynamical evolution of a stellar cluster with initial subclustering and found that subclustering increased both the fraction of stars ejected during the relaxation of the cluster and the formation of binary stars. However, in the last two decades, during which observational studies of star-forming regions have dramatically increased in number, much more attention has been paid to the dynamical evolution of young stellar clusters.

Kroupa (1995a) hypothesised that all stars were formed in binary systems and studied whether such initial conditions could lead to the observed field population of binaries. With an initial period, P, distribution that was flat in log P ranging from 10³ to 10^7.5 d, Kroupa found that subsequent dynamical interactions in dense stellar clusters could reduce the binary fraction from 100 per cent and modified the period distribution to match the properties found for field stars. On the other hand, simulations that begin only with binaries cannot reproduce the observed frequencies or properties of triples (Portegies Zwart et al. 2004) and simulations that begin with primordial binaries and triples cannot reproduce the observed frequencies of quadruples (van den Berk, Portegies Zwart & McMillan 2007). From his simulations of modest-sized clusters containing 200 binaries, Kroupa (1995b) also found that the clusters would be surrounded by a halo of ejected stars with a lower binary fraction. Kroupa (1998) found that ejected binary systems tended to be massive with a preference for equal-mass components. He also found that the mean system mass was approximately independent of ejection velocity from 2–30 km s⁻¹, with dynamical ejection providing relatively massive stars as far as 300 pc from the cluster within 10 Myr (see also Sterzik & Durisen 1995). The binary fraction decreased monotonically with ejection velocity for low-density clusters, but for dense clusters the dependence was more complex. However, long-period systems (periods >10⁴ d) could not be ejected with large velocities (>30 km s⁻¹).

More recently, attention has turned to the effects of gas dispersal on the evolution of young clusters and, following the pioneering work of Aarseth & Hills (1972), to the effects of substructure. Simple analytical models indicate that a young cluster requires a star formation efficiency of η > 50 per cent to survive the rapid dispersal of gas as a bound cluster (Hills 1980). However, more sophisticated semi-analytic and numerical models show that the survival of a bound cluster is more complex. If the gas is removed on a time-scale significantly longer than the crossing time, then η as low as 20 per cent can still result in a bound cluster (Tutukov 1978; Elmegreen 1983; Mathieu 1983; Verschueren & David 1989; Geyer & Burkert 2001; Baumgardt & Kroupa 2007) and even if the gas is removed quickly many stars can be lost, but a bound core can remain (Lada, Margulis & Dearborn 1984; Goodwin 1997; Adams 2000; Geyer & Burkert 2001; Kroupa, Aarseth & Hurley 2001; Boily & Kroupa 2003a,b; Baumgardt & Kroupa 2007). For example, Geyer & Burkert (2001) and Boily & Kroupa (2003b) find that with η= 0.4 and rapid expulsion around 30–40 per cent of the initial cluster could remain bound. Baumgardt & Kroupa (2007) find that for instantaneous gas removal η > 1/3 is required for such a remnant cluster to be formed, but that if the gas is removed over several crossing times, around 50 per cent of the cluster can remain bound even for η as low as 15–20 per cent. Finally, Goodwin (2009) notes that the survivability of a cluster depends not only on the star formation efficiency, but on the cluster's dynamical state. In particular, if the gas removal occurs when the cluster is still dynamically cool (subvirial), then cluster survivability can be greatly increased.

The dynamical state of a cluster is also important for mass segregation. Bonnell & Davies (1998) showed that the Orion Nebula Cluster (ONC) is not old enough for the observed segregation of its massive stars to be due to dynamical segregation alone. Instead, they argued that primordial segregation was required. Recently, however, McMillan, Vesperini & Portegies Zwart (2007), Fellhauer, Wilkinson & Kroupa (2009) and Moeckel & Bonnell (2009) have shown that if the cluster is formed via the merging of subclusters, mass segregation can be accelerated. Similarly, Allison et al. (2009a) investigated the evolution of dynamically cool clusters with fractal initial structure and found that dynamical mass segregation of the most massive stars can occur quickly.

However, the above studies have begun with artificial initial conditions. Assumptions must be made about the initial structure of the stellar clusters (e.g. Plummer spheres, fractal or subclustered stellar distributions), the virial ratio (e.g. in virial equilibrium, subvirial so they collapse, or supervirial so they expand), the initial populations of binary (e.g. 100 per cent binaries) and higher order multiple systems (typically, no triples or higher order systems, with the exception of van den Berk et al. 2007), and the properties of these multiples such as their distributions of separations and mass ratios.

Recently, it has become possible to perform hydrodynamical simulations of the collapse of a molecular cloud to form a star cluster containing in excess of 1000 of stars and brown dwarfs while resolving the opacity limit for fragmentation (Bate 2009a). Such calculations allow detailed comparisons to be made with the observed properties of young clusters, such as the stellar initial mass function (IMF), multiplicity, the properties of binary and multiple systems, and the global properties of stellar clusters. Bate's calculation included gravity and hydrodynamics using a barotropic equation of state and used sink particles to model the stars and brown dwarfs. The sink particles had accretion radii of only 5 au, allowing circumstellar discs to be resolved down to ≈10 au in radius, while binaries and multiple systems with separations as close as 1 au could be followed. Despite the limited amount of physics included in the calculations, many of the stellar properties obtained from the calculation were in good agreement with the observed properties of stellar systems. For example, the observed stellar multiplicity as a function of primary mass was found to be well reproduced, as were the separations and mass ratios of low-mass stars and brown dwarfs. Even the distribution of the orientation of the orbital planes of triple stellar systems was found to be in good agreement with observations. The two main deficiencies of the calculation were that it produced a ratio of brown dwarfs to stars well in excess of that which is observed, and solar-type binaries were found to have a deficit of low-mass ratio systems. The first of these is likely to be due to the neglect of radiative feedback from the protostars (Bate 2009b; Offner et al. 2009), while the solution to the latter problem is currently unclear. However, the implication of Bate's calculation is that many stellar properties result primarily from dissipative gravitational dynamics and that additional physical processes such as radiative transfer and magnetic fields may not be very important for determining these properties.

One objection that could be raised to Bate (2009a) is that, due to computational limitations, the formation of the cluster could only be followed for 2.85 × 10⁵ yr (1.5 global free-fall times) which was only 1.5 × 10⁵ yr after the first star formed. However, in many cases, the stellar properties were compared to star-forming regions that were typically ∼1 Myr old, or even to field stars. When Bate's calculation was stopped, there were many triple, quadruple and higher order stellar systems, and many of there systems were dynamically unstable. Furthermore, of the 500 M_⊙ of gas originally contained in the molecular cloud only 191 M_⊙ (38 per cent) had been converted to stars by the end of the calculation. Star formation would have continued if the calculation were continued. But even if the gas was dispersed at this point and star formation ceased, there is the question of how the stellar properties would evolve over million-year time-scales and longer.

In this paper, we take the endpoint of the hydrodynamical simulation of Bate (2009a) as the starting point for N-body simulations. The advantage of these starting conditions is that the stellar cluster has been formed self-consistently from the evolution of a gravitationally unstable turbulent molecular cloud. Thus, the structure of the cluster and the properties of binary and multiple systems have all been consistently set up, albeit with a restricted number of physical processes (e.g. no radiative feedback or magnetic fields). Our aims are to determine the long-term (10 Myr) evolution of the cluster and how its properties and those of the stellar multiple systems evolve during the dispersal of the gas in the cluster. The only previous similar work we are aware of is that of Hurley & Bekki (2008) who also performed smoothed particle hydrodynamics (sph) simulations of star formation to generate the initial conditions for N-body simulations. However, their work studied star cluster formation on giant molecular cloud scales, assumed equal-mass stars, and did not resolve primordial multiple systems.

The paper is structured as follows. Section 2 briefly describes the initial conditions and the numerical method for the calculations. The results are discussed in Section 3. In Section 4, we discuss the implications of the results for our understanding of star cluster evolution. Our conclusions are given in Section 5.

2 COMPUTATIONAL METHOD

The calculations discussed in this paper are continuations of a hydrodynamical calculation that began with a turbulent molecular cloud and produced a star cluster. The hydrodynamical calculation was performed using a sph code (Bate, Bonnell & Price 1995) and the cloud was evolved for 1.5 initial cloud free-fall times (2.85 × 10⁵ yr, where the free-fall time was t_ff= 6.0 × 10¹² s or 1.90 × 10⁵ yr). The stars and brown dwarfs were modelled using sink particles. In this paper, we take the sink particles at the endpoint of the hydrodynamical calculation as the starting point for a series of N-body simulations and evolve the stellar cluster to an age of 10 Myr. Note that for all the times we quote, the starting point is taken to be that of the hydrodynamical calculation, not the N-body simulations, unless stated otherwise.

2.1 Initial conditions

The initial conditions for the N-body simulations are taken from the main star cluster formation calculation of Bate (2009a). The initial conditions for the hydrodynamical calculation consisted of a 500 M_⊙, uniform, spherical molecular cloud represented by 3.5 × 10⁷sph particles. The cloud's radius was 0.404 pc (83 300 au). At the initial temperature of 10 K, the mean thermal Jeans mass was 1 M_⊙ (i.e. the cloud contained 500 thermal Jeans masses). Although the cloud was uniform in density, an initial supersonic turbulent velocity field was imposed in the same manner as Ostriker, Stone & Gammie (2001). A divergence-free random Gaussian velocity field was generated with a power spectrum P(k) ∝k⁻⁴, where k is the wavenumber. In three dimensions, this results in a velocity dispersion that varies with distance, λ, as σ(λ) ∝λ^1/2, in agreement with the observed Larson scaling relations for molecular clouds (Larson 1981). The velocity field was normalized so that the kinetic energy of the turbulence equalled the magnitude of the gravitational potential energy of the cloud. Thus, the initial rms Mach number of the turbulence was M= 13.7. The ‘turbulence’ decayed as the calculation was evolved; there was no ‘turbulent driving’. The cloud was allowed to evolve freely into the vacuum surrounding it; there were no boundary conditions applied to the simulation.

2.1.1 Sink particles

To model the thermal behaviour of the gas without performing radiative transfer, the sph calculation used a barotropic equation of state (see Bate, Bonnell & Bromm 2003 for further details). This equation of state kept the temperature at 10 K for densities <10⁻¹³ g cm⁻³, while above this density the temperature increased with density in order to mimic the opacity limit for fragmentation (Low & Lynden-Bell 1976; Rees 1976; Silk 1977a,b; Boyd & Whitworth 2005). When the central density of a fragment exceeded ρ_s= 10⁻¹¹ g cm⁻³, it was replaced by a sink particle (Bate et al. 1995). In the main calculation of Bate (2009a), sink particles were formed by replacing the sph gas particles contained within r_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 fell within this radius was accreted by the point mass if it was bound and its specific angular momentum was less than that required to form a circular orbit at radius r_acc from the sink particle. Sink particles interacted with the gas only via gravity and accretion. The gravitational acceleration between two sink particles was Newtonian for r≥ 4 au, but was 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 had in the hydrodynamical simulation even if, in reality, a binary's orbit would have been hardened. Sink particles were permitted to merge if they passed within 2 R_⊙ of each other; only one merger occurred during the simulation.

2.1.2 The end state of the hydrodynamical calculation

The end state of the hydrodynamical calculation consisted of a cluster of 1253 stars and brown dwarfs with a combined mass of 191 M_⊙ while 309 M_⊙ of gas remained. Note that Bate (2009a) states that there were 1254 objects, but, in fact, it was not noticed when the paper was written that two of these had merged and had been replaced by a single object. This has no significant effect on the statistics that were reported in the original paper. The final cluster was created from the merger of around half a dozen subclusters that formed in gaseous filaments that were produced in the turbulent gas. The subclusters formed in a dynamically cool configuration and fell together to produce a single cluster surrounded by an extended halo of unbound, dynamically ejected stars and brown dwarfs. This cluster was very dense, with a half-mass radius of only 10⁴ au (0.05 pc) and a velocity dispersion of around 4 km s⁻¹, with the velocities of the ejected halo objects extending beyond 20 km s⁻¹. Many of the stars were in binary or higher order multiple systems and the stellar multiplicity was found to be a strongly increasing function of primary mass, with values in good agreement with observed values. Many of the trends for the separation and mass ratio distributions of observed binaries were also reproduced [e.g. the trend for very-low-mass (VLM) binaries to have smaller separations and mass ratios nearer unity than for low-mass and solar-type stars], with the exception that the unequal mass solar-type binaries were underproduced.

Further details regarding the initial conditions, the sph code, the evolution of the hydrodynamical star cluster formation calculation and the statistical properties of the stars and brown dwarfs can be found in Bate (2009a). We now turn to the question of how this initial cluster evolves on 10 Myr time-scales.

2.2 The N-body calculations

The final positions, masses and velocities of the sink particles in the hydrodynamic simulation are taken as the initial conditions for a gravitational N-body simulation. The integration was performed with the nbody6 code (Aarseth 2003), which uses a fourth-order Hermite integrator. We advanced the system for 10 Myr from the end of the hydrodynamic simulation, and the relative energy error was of the order of 10⁻⁵. It is usual practice when performing N-body calculations with a modest number of stars to generate many sets of initial conditions and consider ensemble quantities, so that general trends can be separated from chaotic variations in the integrations. In this case, with initial conditions generated not by random realizations of assumed bulk cluster properties but by a self-consistent and expensive hydrodynamical calculation, we are forced to take a slightly different approach to disentangling general results from those particular to the exact initial conditions.

10 sets of perturbed initial conditions were generated by taking the initial conditions and randomizing each body's position by 10⁻⁴ times its radius from the cluster's centre of mass. That is, each position vector r was adjusted so that r→r+ 10⁻⁴rw, with w a randomly oriented unit vector. The initial specific potential energy of each body was thus perturbed by an amount greater than the integration accuracy. Stars in binaries were treated as single bodies for the perturbation, so that their initial binary properties remained the same. By comparing results of the main integration to those of the perturbed versions of the initial conditions, we can ascertain which long-term results are robust and general to the type of cluster generated by the hydrodynamic simulation, and which are sensitive to the exact state of the cluster.

The main calculation and the 10 randomly perturbed calculations are all run assuming that the gas left at the end of the hydrodynamical calculation is instantaneously removed. This is the most extreme gas removal history possible and, based on previous work, should result in the least-bound cluster remnant possible (Lada et al. 1984; Goodwin 1997; Adams 2000; Geyer & Burkert 2001; Kroupa et al. 2001; Boily & Kroupa 2003a,b; Baumgardt & Kroupa 2007). If the gas is removed over a finite period of time, we would expect more of the cluster to remain bound and it might also affect the frequencies and properties of the multiple stellar systems. To investigate how much of a difference slow removal of the gas makes, we have performed three additional calculations that include an additional gravitational potential to model the contribution of the gas.

The gas mass at the end of the hydrodynamic calculation is 309 M_⊙. While the actual gas distribution is not spherical (see fig. 1 of Bate 2009a), we approximate it as a Plummer sphere with scale radius 0.25 pc so that the half-mass radius is similar to that of the actual gas distribution (see Fig. 1). The mass of the Plummer potential is decreased with time, t, from the beginning of the N-body calculations so that the gas mass is given by

1

where M₀ is the original gas mass (309 M_⊙) and τ_gas is a gas removal time-scale. We ran three cases: moderately fast gas dispersal, where the mass was reduced to 10 per cent of its initial value after 0.3 Myr of N-body evolution (total age 0.6 Myr); slow dispersal, where the mass was reduced to 10 per cent at 1 Myr of N-body evolution (total age 1.3 Myr); and an (unphysical) test case where the gas mass remained constant throughout the calculation.

Figure 1

The solid line gives the cumulative radial distribution of the 309 M_⊙ of gas that remained at the end of the hydrodynamical calculation of Bate (2009a), although the gas distribution is not spherically symmetric. Note that the half-mass radius of the stellar cluster is only 0.05 pc. The dotted line gives the cumulative radial distribution of a Plummer sphere with the same total mass and a characteristic radius of 0.25 pc. The two distributions have similar half-mass radii.

A summary of the main properties of the calculations discussed in the next section is presented in Table 1. The plots and discussion that follow are based on the unperturbed initial conditions with instantaneous gas removal, with comparison to the randomized integrations and those with gas potentials made where appropriate.

Table 1

A summary of the calculations discussed in this paper.

Calculation	End time (Myr)	Remnant cluster		Mass segregation	Final numbers of systems
		Mass (per cent)	Radius (pc)		Single	Binary	Triple	Quadruple
Hydrodynamical (Bate 2009a)	0.3	–	–	≈10 most masssive stars	905	90	23	25
Instant gas removal	10	30	2	–	956	98	23	8
T0.1= 0.3 Myr	10	30	1.8	≳1 M⊙	941	100	20	13
T0.1= 1 Myr	10	40	1.3	≳0.3 M⊙	962	97	19	10
No gas removal	10	50	0.8	≳0.3 M⊙	966	93	11	17

Note. The four N-body calculations with different gas removal time-scales take their initial conditions as the stars (sink particles) at the end of the hydrodynamical calculation of Bate (2009a). T_0.1 is the time from the start of the N-body evolution at which the gas mass has fallen to 10 per cent of its initial value. In addition to the calculations listed above, the instantaneous gas removal case was run 10 times with randomly perturbed initial conditions. In the table, we summarize the properties of the remnant clusters that are left at the end (10 Myr) of the N-body calculations (their radii and the mass fractions of the stars they contain). We also comment on the mass segregation that is observed (e.g. stars with masses >0.3 M_⊙ at the end of the case with no gas removal). Finally, we give the numbers of single and multiple systems at the end of the hydrodynamical calculation of Bate (2009a) and each of the unperturbed N-body calculations. The N-body evolution always results in a significant decrease in the number of quadruple systems and small increases in the numbers of single and binary systems. Slower gas removal time-scales generally produce more compact and populous remnant clusters that show more mass segregation.

3 RESULTS

3.1 Evolution of cluster properties

3.1.1 Substructure

In the hydrodynamic simulation of Bate (2009a), the stars form in a structured fashion, with smaller subclusters merging to form a final cluster consisting of a tightly bound core with radius ≈0.05 pc surrounded by an expanding halo of ejected stars. The reader interested in the hydrodynamical evolution of the cluster is encouraged to consult Fig. 1 of Bate (2009a). In Fig. 2 of this paper, we display the global evolution of the stellar cluster during the 10 Myr of N-body evolution. Each star is represented by an open circle with its area proportional to the star's mass. On large scales, the stars in the halo are seen to expand to distances exceeding 100 pc from the cluster in less than 10 Myr. On small scales, a bound cluster remnant containing approximately 30 per cent of the stars expands from a radius of ≈0.05 to ≈2 pc over several Myr.

Figure 2

The N-body evolution of the unperturbed star cluster with instantaneous gas removal. Each star is plotted using an open circle with the area of the circle proportional to the star's mass. Stars with masses ≥3 M_⊙ are plotted in red, while those with masses 1 ≤M < 3 M_⊙ are plotted in blue. Binaries and multiples can be seen as almost concentric circles when the components have different masses. We show the structure on 200, 20 and 5 pc scales at four different times. On the largest scales, it can be seen that some ejected stars in the halo reach distances >100 pc in less than 10 Myr and that these stars are a mixture of massive and low-mass stars. On the smallest scales, the remnant bound cluster containing approximately 1/3 of the stars expands from its initial radius of ≈0.05 to ≈2 pc. The most massive star and a small cluster of ≈20 companion stars is ejected from the remnant cluster at ≈1 Myr and can be seen to the lower left of the main cluster remnant in the lower two 5 pc panels.

The degree of subclustering can be quantified by the parameter Q (Cartwright & Whitworth 2004). Values of Q < 0.8 indicate substructure and can be associated with a fractal dimension, while values greater than 0.8 are associated with radial density variation. Calculating Q requires both an area and radius to be assigned to a cluster; as the cluster at early times is decidedly non-spherical, we calculate the area as the convex hull of the stars' positions projected on to a plane, and take the radius to be that of a circle having the same area. The evolution of Q over the first 1.5 Myr of the cluster evolution is shown in Fig. 3.

Figure 3

The value of the parameter Q over the first 1.5 Myr. At each time the quantity is determined for three orthogonal projections of the cluster, shown as the points, and the means of those values are connected by the lines. The black line gives the evolution of the entire stellar cluster and halo, while the grey line gives the evolution of the inner 3 pc region. The inner region evolves from a subclustered stellar distribution (Q < 0.8) to a core and halo structure (Q > 1) during the hydrodynamical calculation. When the gas is removed, the remnant core expands to fill most of the inner 3 pc and Q decreases to approximately unity.

At the earliest measured times the cluster consists of 154 stars, and the value of Q≈ 0.4 corresponds to an estimated fractal dimension of 1.5, i.e. highly substructured. Over the next ∼4 × 10⁴ yr, the early subclusters merge, the number of stars increases to 579 and Q rises above 0.8. This indicates a fractal dimension of 3, i.e. a distribution lacking any significant substructure. By the end of the hydrodynamical calculation, the number of stars has reached its final value of 1253, and Q has risen in excess of 1.5. In Cartwright & Whitworth (2004), larger values of Q are identified as corresponding to an increasingly steep radial density profile of a spherical cluster. In our case, the cluster is not well-fit by a single power-law density profile, but rather has a core–halo structure, which complicates a straightforward interpretation of Q. Because the cluster is characterized by a bound remnant surrounded by an expanding halo, the exact value we calculate is quite sensitive to the radial extent of stars that we consider. In Fig. 3, following gas removal, two lines are plotted. The black line includes all the stars, arranged with a core–halo structure, while for the grey line the analysis has been restricted to the central 3-pc of the cluster. In this central region, the cluster core expands upon gas removal to occupy most of the inner region and the value of Q drops to unity, reflecting an absence of significant structure and the fact that the halo is only present on larger scales. For our purposes, it suffices to note that Q remains quite high for the entire 10 Myr run, showing the persistence of the core–halo structure.

We emphasize that the evolution of the cluster from a highly clustered to a core–halo structure occurred while the cluster was still embedded and forming; the number of stars increases by nearly an order of magnitude during the transition from a fractal dimension ≈1.5 to ≈3.0. This highlights the importance of coupled gas and gravitational dynamics in young cluster evolution, and the difficulty in assigning realistic initial conditions to purely N-body studies.

3.1.2 Lagrangian radii and stellar velocities

The radial mass structure of the cluster can be seen in the evolution of the Lagrangian radii, i.e. the radii encompassing a constant percentage of the cluster's mass. This calculation is dependent on what one considers the cluster centre. At each time, we identify the star which is located in the region of the greatest stellar density, which we define as the star with the smallest sphere encompassing 50 neighbours, and set this star as the origin. The behaviour of these diagnostics depends sensitively on the details of the integration, especially the behaviour of the most massive stars. We first discuss the main features of the Lagrangian radii for the unperturbed initial conditions, shown in Fig. 4, and then note how these vary among the randomly perturbed clusters and those with a gas potential.

Figure 4

Lagrangian radii during the span of the main N-body calculation. From bottom to top, solid curves enclose 0.025, 0.05, 0.1. 0.2, 0.3, 0.4, 0.5, 0.7 and 0.9 of the total stellar mass. The dotted line shows the median radius of the five most massive stars.

The first Myr of N-body evolution is marked by general expansion of the outer 70 per cent of the cluster's mass, with the Lagrangian radii increasing by an order of magnitude. In contrast, the inner 10 per cent of the mass expands by only a factor of 2. Between 1 and 2 Myr, these inner radii expand by an order of magnitude, catching up to the expansion of the outer radii. At later times, there is a divergence as the outer radii continue to expand, while the inner radii (certainly the inner 20 per cent, and after 8 Myr the inner 30 per cent) halt their growth, leaving a remnant cluster core at the centre of an expanding halo. The evolution of this remnant cluster core is further explored in Fig. 5 where we plot the velocity dispersion of the stars in the remnant as a function of both time and of the remnant cluster's radius. As expected for a system near virial equilibrium, the velocities of the stars decrease as the remnant expands roughly as R^−1/2.

Figure 5

The evolution of the stellar velocity dispersion and the radius of the bound remnant cluster in the main N-body calculation. In the top panel, we give the time evolution of the radius containing 30 per cent of the stars (dashed line), centred on the star at the location of the highest stellar density. It can be seen that the remnant cluster expands by a factor of ≈20 in size. The velocity dispersion of the stars in the remnant (taking the centre of mass velocity for all binaries and counting them only once) is given by the solid lines in both panels. The velocity dispersion decreases from ≈4 to ≈1 km⁻¹. The spike at ≈1 Myr is associated with the ejection of the most massive binary from the cluster. In the lower panel, the velocity dispersion is plotted as a function of the remnant cluster's radius. The dotted line has a slope of −1/2, as expected if the remnant cluster remains close to virial equilibrium as it expands.

The 1 Myr delay in the expansion of the inner 10 per cent of the mass is explained by the behaviour of the most massive stars. The five most massive stars in the cluster are by themselves approximately 10 per cent of the total stellar mass, and if they are concentrated in the cluster centre, then they will dominate the inner Lagrangian radii. To give an idea of these stars' radii from the cluster centre, we also plot in Fig. 4 the median radius of the five most massive stars as a dotted line. At about 1 Myr, a massive binary, comprised of a 5.35 M_⊙ primary (also the most massive star in the cluster) and a 3.54 M_⊙ secondary, is ejected from the cluster centre (see the lower two right-hand panels of Fig. 2). This event also appears in the velocity dispersion evolution depicted in Fig. 5. The massive binary is accompanied by a small group of ≈20 companion stars. The inner 10 per cent Lagrangian radii track this binary and its companions for a few ×10⁴ yr, at which point they begin to reflect the overall cluster structure rather than the positions of a few massive stars.

This feature, a delayed expansion of the inner Lagrangian radii due to a clustering of the most massive stars, is seen in nine of the 10 randomly perturbed simulations, though the time of the ejection ranges from ∼4 × 10⁴ yr to 2 Myr after the start of the N-body evolution. In one case the massive stars remain in the cluster centre throughout the 10 Myr. In other words, in only one of our 11 runs with instantaneous gas removal do the massive stars remain continually associated with the centre of the cluster for more than 2 Myr. In two cases, the massive stars return to the centre of the cluster, one of which sees a second ejection. In that perturbed realization, after the second ejection of massive stars at about 3 Myr, the inner Lagrange radii expand freely. While a cluster remnant is still identifiable, it is globally dispersing unlike any of the other realizations. Another feature is the behaviour of the 20 per cent radius in the unperturbed case, which is fairly flat after its initial expansion. In six of the perturbed runs, this radius increases by a factor of 3 or less between 2 and 10 Myr; in two cases, this radius is flat for most of the evolution but grows over the final 2 Myr; and, in one case, it expands freely after the ejection of the massive stars. The global evolution of the cluster's mass profile thus appears to be surprisingly sensitive to the initial conditions. We can say with some confidence that the massive stars are unlikely to remain associated with the main cluster remnant, but the final characteristics of that remnant (i.e. its radius, mass and its stability or expansion) display significant variation.

As expected, the N-body calculations that include a gas potential and slower gas removal result in smaller final Lagrangian radii, and larger membership of the remnant cluster. In Fig. 6, we show the evolution of the Lagrangian radii in the slow gas dispersal case (where the gas mass is reduced to 10 per cent of its initial value after 1 Myr of N-body evolution). Comparing this with Fig. 4, we find that somewhat more of the stellar mass (about 40 per cent compared to 30 per cent) remains in a remnant cluster and this remnant cluster is somewhat smaller in the slow gas dispersal case (radius approximation 1.3 pc rather than 2 pc). The effect of the gas potential is more clearly demonstrated in Fig. 7, where we show the Lagrangian radii at 10 Myr for each of the unperturbed cases, which we label by the time-scale, T_0.1, for the gas mass to be reduced to 10 per cent of its initial value. The inner Lagrangian radii all display the expected trend, with longer gas-removal times yielding a more compact cluster and populous cluster. At 10 Myr, the instantaneous, 0.3 Myr, 1.0 Myr and no gas removal cases have 552, 595, 756 and 944 stars within 3 pc of the densest star, respectively.

Figure 6

Lagrangian radii during the span of the N-body calculation that includes a gas potential which decreases to 10 per cent of its initial value over 1 Myr of evolution (total age 1.3 Myr). From bottom to top, solid curves enclose 0.025, 0.05, 0.1. 0.2, 0.3, 0.4, 0.5, 0.7 and 0.9 of the total stellar mass. The dotted line shows the median radius of the five most massive stars. In comparison with Fig. 4, the remnant cluster is smaller and more populous, while more of the most massive stars have been ejected.

Figure 7

Lagrangian radii at the end of the N-body calculations (10 Myr) as a function of the gas removal time-scale, T_0.1, which gives the time from the beginning of the N-body evolution at which the gas potential has been reduced to 10 per cent of its initial value (0, 0.3, 1 Myr and no gas removal). Longer gas removal times yield more compact and populous clusters.

3.2 Mass segregation

Mass segregation is seen in several young clusters, typically appearing as a concentration of the most massive stars (e.g. the ONC or Mon R2; Hillenbrand & Hartmann 1998; Carpenter et al. 1997) rather than a general trend affecting all mass ranges, as would be expected in older, dynamically mature clusters. Recent works (McMillan et al. 2007; Fellhauer et al. 2009; Allison et al. 2009a; Moeckel & Bonnell 2009) show that a structured mode of star formation, with merging substructures forming the final quasi-spherical cluster, can lead to segregation affecting only the massive stars, mirroring observations. Moeckel & Bonnell (2009) noted that the final state of the Bate (2009a) hydrodynamic simulation showed a similar mass segregation signal to the ONC. Here, we examine how this rapid mass segregation evolves on longer time-scales.

With the realization that mass segregation in very young clusters is qualitatively different from the segregation that develops as a spherical cluster evolves over many relaxation times, new methods of quantifying segregation have been introduced. We use two diagnostics here: first, a standard technique, the cumulative radial distribution of stars in different mass bins; secondly, a variation of the minimum spanning tree (MST) method introduced in Allison et al. (2009b), in which we measure the effective surface density of the most massive stars and compare it to the expected surface density. This is the method used in Moeckel & Bonnell (2009), and full details are found there. Briefly, the expected surface density of n stars in the cluster is found by calculating the area of the convex hull, H_c(n), for many samples of n random stars. The median, ±1σ, ± 2σ and ±3σ of the quantity H_c(n)/n define the ‘normal’ surface density of n random stars. This distribution is then compared to the surface density of the n most massive stars. With both of these methods, we restrict our analysis to stars within 3 pc of the star at the location of the greatest stellar density. This radius encompasses the final cluster remnant at 10 Myr, and prevents very distant escaping stars from interfering with a sensible interpretation of the results.

In Fig. 8, we show the state of the cluster's segregation at three times; the beginning of the N-body simulation, just after 1 Myr when the most massive stars are about to be ejected from the cluster centre, and at around 2 Myr, which is representative of later times. At the earliest time, the cumulative distribution shows very little evidence of mass segregation (as noted by Bate 2009a), while the surface density analysis shows that the 10 most massive stars have a surface density higher than the 3σ value expected for the cluster. At around 1 Myr, the cumulative distribution shows a clear difference in the radial distribution of stars above 1 M_⊙. The surface density of the four most massive stars is still very high, and both plots reflect the lowered surface density of the cluster after its initial expansion. It is interesting to note that the surface density of the cluster between 1 and 2 Myr has dropped by only a factor of ∼2, while the inner Lagrangian radii (Fig. 4) have increased by over an order of magnitude, emphasizing that the most massive stars were dominating the inner Lagrangian radii at early times.

Figure 8

The mass segregation in the unperturbed cluster with instantaneous gas removal at three times; the end state of the hydrodynamic simulation, just before the ejection of the most massive stars from the cluster centre and at ∼2 Myr. Top row: segregation measured by the cumulative radial distribution of stars in different mass bins, measured from the position of the star in the region of highest stellar density. The lines, from light grey to black, show stars in the mass bins M/M_⊙ < 0.1, 0.1 ≤M/M_⊙ < 0.3, 0.3 ≤M/M_⊙ < 1.0, 1.0 ≤M/M_⊙. Bottom row: segregation measured by the effective surface density of the most massive n stars (see the text for further explanation). The grey line is the median value for the surface density of n random stars, and the shaded regions show the ±1σ, ± 2σ and ±3σ values. The blue line is the surface density of the n most massive stars.

By 3 Myr, mass segregation among the four most massive stars remaining in the cluster has re-established itself in the surface density diagnostic, and remains for the remainder of the simulation. However, in contrast to the segregation at 1 Myr, this is due to the presence of two massive and well-separated binaries, and not a real clustering of four stars. The disappearance and return of the mass segregation signal is dependent on what one considers part of the cluster versus the field. In this case, when dispersing massive stars cease to be considered part of the cluster (by our 3 pc selection criterion), the remaining massive stars are left in a configuration that appears to be clustered. In addition to the cluster membership determination, the amount of segregation seen in a cluster depends on the viewing angle. A projection that minimizes the separation between two massive binaries will overemphasize the clustering. This is important to bear in mind when considering clustering involving only a handful of massive stars.

In Table 2, we examine the radial properties of the remnant cluster and halo of ejected stars using logarithmically spaced bins. In agreement with the diagnostics discussed above, there is no evidence of mass segregation. Some of the most massive stars in the system are found in each radial bin. The velocities of the ejected stars in the halo increase with radius outside of 3 pc from the cluster centre (i.e. faster moving stars have travelled further from the cluster centre).

Table 2

Radial properties of the remnant cluster (radius ≈2 pc) and stellar halo (extending beyond 200 pc) at the end of the unperturbed N-body evolution with instantaneous gas removal (age 10 Myr).

Quantity/Distance range	<1 pc	1–3 pc	3–10 pc	10–30 pc	30–100 pc	>100 pc
Median mass [M⊙]	0.046	0.033	0.035	0.047	0.046	0.050
Upper quartile mass [M⊙]	0.14	0.09	0.11	0.14	0.15	0.21
Maximum mass [M⊙]	2.9	5.3	2.3	1.7	3.7	3.7
Velocity dispersion [km s−1]	0.8	1.3	0.7	2.1	5.7	17
Number objects	276	236	223	276	169	73
Number binaries	39	20	27	32	10	4
Binary fraction	0.16	0.09	0.14	0.13	0.06	0.06

Note. This table can be compared with table 2 of Bate (2009a) for the properties of the cluster at the end of the hydrodynamical evolution. There is no evidence for radial mass segregation of the cluster in terms of the median, upper quartile or maximum masses. Indeed, we note that some of the most massive stars in the cluster are present in each of the radial bins. The binary fraction is more or less uniform, except beyond 30 pc where it drops by about a factor of 2. The velocities of the ejected stars in the halo increase roughly linearly with distance beyond 3 pc (i.e. faster moving stars have travelled further).

All four mass-removal scenarios begin with the same initial stellar positions and, therefore, all begin with the 10 most massive stars being centrally concentrated. While this primordial segregation is transient and only to be observed in the earliest stages of cluster evolution, the more compact remnant clusters that result from slower gas removal have a direct effect on the later mass segregation. Calculated at 10 Myr, the half-mass relaxation times of the calculations with instantaneous gas removal, T_0.1= 0.3 Myr, T_0.1= 1 Myr and no gas removal are 20.8, 12.3, 8.2 and 3.3 Myr, respectively. Thus, whereas the main calculation is dynamically unrelaxed, the two calculations with slow or no gas removal have evolved for more than one relaxation time, which is the time-scale on which mass segregation manifests itself. This is the classical form of mass segregation, due to two-body effects in an almost spherical cluster, in contrast to the primordial segregation that is an imprint of star formation processes and clustered star formation. As such, it should appear in cumulative distribution plots as a general trend, with more massive populations progressively more centrally concentrated.

In Fig. 9, we show the cumulative radial distributions for the main calculation and the three cases with gas potentials at 10 Myr. The case with T_0.1= 0.3 Myr shows the expected signal, with the distribution of the most massive stars beginning to be separated from the lower mass stars. The calculations with T_0.1= 1 Myr and without gas removal show segregation of the second-most massive bin as well. We stress that while the earliest (earlier than 1 Myr for these clusters) observations of segregation are a result of a forming cluster's morphology, this later type of mass segregation is expected due to the dynamical maturity of these clusters, and does not contain information about the star formation process.

Figure 9

The dependence of the mass segregation at 10 Myr on the gas removal time-scale. From left to right, the four calculations are: instantaneous gas removal, T_0.1= 0.3, 1 Myr and no gas removal. The segregation is measured by the cumulative radial distribution of stars in different mass bins, measured from the position of the star in the region of highest stellar density. The lines, from light grey to black, show stars in the mass bins M/M_⊙ < 0.1, 0.1 ≤M/M_⊙ < 0.3, 0.3 ≤M/M_⊙ < 1.0, 1.0 ≤M/M_⊙. The slow and no gas removal cases display classical mass segregation with more massive stars being progressively more centrally condensed (particularly in the two highest mass bins).

3.3 The evolution of multiple systems

Bate (2009a) quantified the fraction of stars and brown dwarfs that were in multiple systems using the multiplicity fraction, defined as

2

where S is the number of single stars within a given mass range and B, T and Q are the numbers of binary, triple and quadruple systems, respectively, for which the primary has a mass in the same mass range. This differs from the companion star fraction (csf) that is also often used and where the numerator has the form B+ 2T+ 3Q. Bate chose the multiplicity fraction following Hubber & Whitworth (2005) who pointed out that this measure is more robust observationally in the sense that, if a new member of a multiple system is found (e.g. a binary is found to be a triple), the quantity remains unchanged. It is also more robust for simulations in the sense that, if a high-order system decays because it is unstable, the numerator only changes if a quadruple decays into two binaries (which is quite rare). Furthermore, if the denominator is much larger than the numerator (e.g. for brown dwarfs where the multiplicity fraction is low), the production of a few single objects does not result in a large change to the value of MF. This was useful for Bate (2009a) because many of the high-order systems in existence at the end of the calculations were likely to undergo further dynamical evolution. Indeed, the long-term evolution of the multiple systems produced by the main calculation of Bate (2009a) is the topic of this section.

In Fig. 10, we plot the multiplicity fraction of the stars and brown dwarfs as a function of stellar mass at the end of the hydrodynamical calculation (i.e. the beginning of the N-body evolution) and at the end of the main N-body calculation at an age of 10 Myr (top left and top right panels, respectively). The top left panel is a reproduction of fig. 15 of Bate (2009a). The mass bins were chosen for comparison with the results of various observational surveys. The filled blue squares give the multiplicity fractions from the calculation, while the surrounding blue hatched regions give the range in stellar masses over which the fraction is calculated and the 1σ (68 per cent) uncertainty on the multiplicity fraction. The black open boxes and their associated error bars and/or upper/lower limits give the results from a variety of observational surveys (see the figure caption).

Figure 10

Multiplicity fraction as a function of primary mass. The top left-hand panel gives the results at the end of the hydrodynamical calculation (0.3 Myr), while the other panels give the results at the end of the N-body evolutions (10 Myr) for three cases: instantaneous gas removal (top right-hand panel), T_0.1= 0.3 Myr (lower left-hand panel) and with no gas removal (lower right-hand panel). The blue filled squares surrounded by shaded regions give the results from the calculations with their statistical uncertainties. The open black squares with error bars and/or upper/lower limits give the observed multiplicity fractions from the surveys of Close et al. (2003), Basri & Reiners (2006), Fischer & Marcy (1992), Duquennoy & Mayor (1991), Preibisch et al. (1999) and Mason et al. (1998), from left to right. There is very little evolution of the multiplicity fraction over the 10 Myr of N-body evolution, even for the extreme case of no gas removal. All sets of results are in reasonable agreement with the observed multiplicity as a function of primary mass.

There is very little evolution of the multiplicity fraction during the main N-body evolution of the cluster with instantaneous gas removal. Both the results at the end of the hydrodynamical calculation and those at the end of the N-body evolution are in good quantitative and qualitative agreement with the observed multiplicity fraction (see Bate 2009a for a full discussion), which is found to be a strongly increasing function of primary mass. We also note that the decrease in binary fraction in the outer parts of the ejected halo of stars found by Bate (2009a) is maintained to ages of 10 Myr (Table 2), but that the transition radius at which the binary fraction begins to decrease with distance has moved from 0.15 to 30 pc as the halo has expanded.

For the randomly perturbed sets of initial conditions, although each calculation differs in detail with respect to which binaries and multiple systems survive or evolve dynamically, any variations in the values of the multiplicity as a function of primary mass are smaller than the statistical uncertainties already plotted in Fig. 10. Typically, for each mass range, one or two more or fewer multiple systems survive in each random realization for primary masses <0.50 M_⊙. For the more massive systems, there tends to be somewhat more variation. For example, for the most massive systems (primary masses >1.2 M_⊙), the multiplicity ranges from 6/22 = 27 per cent to 10/20 = 50 per cent. But the statistical uncertainties are also largest for these multiplicities due to the small numbers of massive stars.

The lower panels of Fig. 10 give the multiplicity fractions for the two of the calculations aimed at investigating the dependence of the results on the gas removal time. The results from the T_0.1= 1 Myr (lower left-hand panel) and no gas removal (lower right-hand panel) calculations are given. It can be seen that the multiplicity fractions decrease slightly overall with increased gas removal time-scales, but only marginally despite the fact that the remnant clusters become significantly more condensed and populous. Even in the extreme case of no gas removal, the decreases in the multiplicity fractions for different primary masses between the instantaneous gas removal and no gas removal cases are smaller than the statistical uncertainties. As for the randomly perturbed sets of initial conditions, only the most massive systems (primary masses >1.2 M_⊙) have a large change in the multiplicity fraction, but, again, the statistical uncertainties are also largest in this mass range due to the small numbers of massive stars.

3.3.1 Triple and quadruple systems

Although the multiplicity fraction does not evolve significantly, this is partially due to the fact that this measurement of the fraction of binaries and multiples was deliberately chosen to be relatively insensitive to the break-up of dynamically unstable multiple systems. In fact, during the N-body evolution of the unperturbed case with instantaneous gas removal, the number of quadruple systems drops from 25 to 8, while the number of triple systems remains the same (though the specific make-up of the 23 systems changes). The number of binary systems increases slightly from 90 to 98 and the number of single objects also increases from 904 to 956 (Table 1). Much of this evolution occurs quite quickly, with the number of quadruple systems having dropped to 10 after 0.3 Myr of N-body evolution (i.e. when the cluster age is 0.6 Myr). Thus, high-order systems (particularly quadruple systems) are converted into lower order systems through dynamical disruption, but the multiplicity fractions themselves are not significantly changed.

The final overall frequencies of triple and quadruple systems are 2.1 and 0.7 per cent, respectively. However, as discussed by Bate (2009a), these frequencies are a strong function of primary mass. At the end of the N-body evolution, there are no quadruple systems with primary masses less than 0.2 M_⊙. There are two surviving VLM (0.03–0.10 M_⊙) triple systems out of a total of 331 systems (giving a frequency of 0.6 per cent). At the end of the hydrodynamical calculation, there were three triples and two quadruple systems with primaries in this mass range. Thus, the fraction of high-order VLM multiples has significantly decreased. For low-mass stars, there are four triple systems out of the 133 systems with primaries in the range of 0.1–0.2 M_⊙ (i.e. 3 per cent) at the end of the N-body evolution, whereas at the end of the hydrodynamical evolution there were four triples and two quadruples. The frequencies of triples or quadruples for M-dwarf primaries with masses 0.2–0.5 M_⊙ is nine out of 96 systems (five triples and four quadruples), very similar to the seven triples and one quadruple found from the 99 M-dwarf primaries surveyed by Fischer & Marcy (1992). At the end of the hydrodynamical evolution, there were five triples and 10 quadruples with primaries in this mass range. For higher mass primaries, the frequencies of triple and quadruple systems are 18 per cent for 0.5–0.8 M_⊙ and 29 per cent for primary masses >0.8 M_⊙, but with large statistical uncertainties. The number of high-order systems with these higher mass primaries is similar after the N-body evolution to that at the end of the hydrodynamical evolution, but whereas the number of triples has increased from 10 to 12, the number of quadruples has dropped from 11 to 4.

The destruction of high-order multiple systems is qualitatively the same for the randomly pertubed realizations and even for the calculations that include a gas potential. In Table 1, we give the total numbers of single, binary, triple and quadruple systems at the end of the hydrodynamical calculation and at the end of the four N-body calculations with different gas removal time-scales. For gas removal time-scales of 1 Myr or less, the numbers are statistically indistinguishable. For the calculation without gas removal, the number of triples is lower, but on the other hand the number of quadruples is somewhat higher than in the cases with gas removal. Closer inspection (see the next section) shows that many of these quadruples are very wide systems in the halo.

3.3.2 Separation distributions

Along with the evolution of the frequencies of multiple systems, we wish to understand how the characteristics of the multiple systems evolve. In Fig. 11, we present the separation (semimajor axis) distributions of the stellar (primary masses greater than 0.10 M_⊙; left-hand panels) and VLM multiples (right-hand panels). These distributions are compared with the survey of Duquennoy & Mayor (1991) and the listing of VLM multiples maintained by N. Siegler, C. Gelino and A. Burgasser at http://vlmbinaries.org/(last updated 2009 July 28), respectively. The filled histograms give the separations of binary systems, while the double hashed region adds the separations from triple systems (two separations for each triple, determined by subdividing the triple into a binary with a wider companion), and the single hashed region includes the separations of quadruple systems (three separations for each quadruple which may be composed of two binary components or a triple with a wider companion). The top row gives the distributions at the end of the hydrodynamical calculation, while the second row gives the distributions at the end of the N-body calculation (10 Myr) with instantaneous gas removal.

Figure 11

The distributions of separations (semimajor axes) of multiple systems with stellar (left-hand panels) and VLM (right-hand panels) primaries. From top to bottom, the four rows give the results: at the end of the hydrodynamical calculation (top row; 0.3 Myr; reproduced from Bate 2009a) and at the end of the N-body evolution (10 Myr) with instantaneous gas removal (second row), gas removal time-scale T_0.1= 1 Myr (third row) and no gas removal (bottom row). The T_0.1= 0.3 Myr case has been omitted as the results lie between those of the instantaneous gas removal and T_0.1= 1 Myr cases. The solid, double hashed and single hashed histograms give the orbital separations of binaries, triples and quadruples, respectively (each triple contributes two separations, each quadruple contributes three separations). In the stellar graph, the curve gives the G-dwarf separation distribution (scaled to match the area under the histogram in each case) from Duquennoy & Mayor (1991). In the VLM systems graph, the open black histogram gives the separation distribution (scaled to match the number in the 10–100 au range of the calculated histogram in each case) of the known VLM multiple systems maintained by N. Siegler at http://vlmbinaries.org/(last updated on 2009 July 28). The vertical dotted line gives the resolution limit of the hydrodynamical calculation as determined by the gravitational softening and accretion radii of the sink particles.

We note that the gravitational softening and finite sink particle accretion radii used in the hydrodynamical calculation of Bate (2009a) result in a pile-up of stellar binaries with separations of a few au. This was investigated by Bate who reran part of the hydrodynamical calculation without gravitational softening and using smaller accretion radii and, as expected, the pile-up disappeared and a bell-shaped separation distribution was recovered (but with a much smaller number of objects because the rerun calculation could not be followed for as long).

Comparing the upper two left-hand panels in Fig. 11 for the stellar multiple systems, we find that, as mentioned earlier, many of the quadruple systems are broken up, resulting in a marked decrease in the contribution of the quadruple systems to the separation distributions. There is also an increase in the number of triple systems with component separations <10 au. The other significant change is in the number of very wide systems (separations 10⁴–10⁵ au). While at the end of the hydrodynamical calculation, there were only two separations in this range (one binary, and one from the wide component of a quadruple), at 10 Myr there are seven binaries, five triples and two quadruples with separations in this range. These wide systems are found in the halo of ejected objects and are formed when two ejected objects have similar velocities and as they depart from the cluster, they happen to be weakly mutually bound. The overall effect of the break-up of stellar quadruples with intermediate separations (1–1000 au) and the formation of very wide systems is to significantly broaden the separation distribution. In fact, the separation distribution from 10–10⁵ au appears to be even flatter than that observed by Duquennoy & Mayor (1991), and more consistent with the flat separation distribution of wide (500–5000 au) binaries found for pre-main-sequence stars by Kraus & Hillenbrand (2009).

For the VLM systems (primary masses <0.1 M_⊙), the evolution is less dramatic (upper two right-hand panels of Fig. 11). Many of the triples and all of the quadruple systems present at the end of the hydrodynamical calculation have been broken up by the end of the N-body evolution. However, these comprise a smaller fraction of the initial separation distribution than for the stellar systems. There are two main effects of the long-term evolution. First, all the systems with separations greater than 60 au at the end of the hydrodynamical calculation have been destroyed. Secondly, two extremely wide VLM binaries have been formed in the dispersing halo population: a 26 000 au system consisting of 40 and 14 M_J brown dwarfs and a 34 000 au system consisting of 39 and 10 M_J brown dwarfs. This is in good agreement with the observed separation distribution of VLM binaries: field VLM binaries are almost always found to have projected separations ≲40 au (Martín et al. 2000; Reid et al. 2001; Close et al. 2003; Burgasser et al. 2007) with only a handful of wider field systems found till date (220, 1800, 5100, 6700 au; Billères et al. 2005; Artigau et al. 2007; Caballero 2007; Radigan et al. 2009), while four out of five VLM binaries with separations in the 100–1000 au range have been found in the star-forming regions (Luhman 2004; Jayawardhana & Ivanov 2006; Close et al. 2007; Béjar et al. 2008).

The evolution of the separation distributions discussed above is qualitatively similar for all the randomly perturbed realizations, though quantitatively there is some variation. The number of wide (10⁴–10⁵ au) systems displays significant variation. In the unperturbed case, seven binary, five triple and two quadruple separations fell within this separation range. The mean values for all 11 realizations are: 6.8 ± 3.5, 3.6 ± 2.7 and 4.5 ± 2.5 with the total number of separations in this range being 15 ± 7.2 with the minimum being five and the maximum being 32. Thus, the unperturbed case is similar to the mean case, but sometimes there can be a lot more or a lot fewer wide multiple systems formed. This is not unexpected, since these systems require stars or systems to be ejected with similar velocities, something that is very sensitive to small perturbations.

In the lower two rows of Fig. 11, we present results from the N-body calculations with gas potentials at 10 Myr. The third row presents the results from the calculation with a gas removal time-scale of T_0.1= 1 Myr, while the bottom row presents the results from the calculation with no gas removal. The T_0.1= 0.3 Myr case has been omitted because, as might be expected, the results are nicely bracketed by the cases of instantaneous gas removal and T_0.1= 1 Myr. As expected, longer gas removal time-scales destroy more wide systems (100 − 10⁴ au), with triple and quadruple systems preferentially disrupted. However, there is very little difference between the calculations with different gas removal time-scales at separations less than 100 au and the formation of significant numbers of very wide systems (separations >10⁴ au) discussed above still occurs even without gas removal. For the VLM systems, there is also very little dependence on the gas removal time-scale. In all cases, the main evolution is that the number of systems with separations of 100–10⁴ au decreases substantially and, in all but the no gas removal case, a few very wide (>10⁴ au) systems are formed in the halo. Even in the extreme case of no gas removal at 10 Myr, the separation distribution is in good agreement with observations.

3.3.3 Mass ratio distributions

As we have seen, subsequent dynamical evolution of the Bate (2009a) cluster results in the formation of very wide systems and evolution of the separation of VLM binaries that is in good agreement with that observed. How do the mass ratio distributions evolve?

In Fig. 12, we present the mass ratio distributions of binaries with primaries that have masses ≥0.5 M_⊙ (left-hand panels), are M-dwarfs with masses 0.1 ≤M < 0.5 M_⊙ (centre panels) or are VLM objects (right-hand panels). Again, the top panels are at the end of the hydrodynamical calculation (see also fig. 19 in Bate 2009a), while the other panels give the results at the end of the N-body evolutions at an age of 10 Myr. The case with instantaneous gas removal is given by the middle panels, while the lower panels give the results with no gas removal. We compare the M-dwarf mass ratio distribution to that of Fischer & Marcy (1992), and the higher mass stars to the mass ratio distribution of solar-type stars obtained by Duquennoy & Mayor (1991). The VLM mass ratio distribution is compared with the listing of VLM multiples maintained by N. Siegler, C. Gelino and A. Burgasser at http://vlmbinaries.org/(last updated 2009 July 28). The figure only includes the mass ratios of binaries, but we include binaries that are components of triple and quadruple systems. A triple system composed of a binary with a wider companion contributes the mass ratio from the binary, as does a quadruple composed of a triple with a wider companion. A quadruple composed of two binaries orbiting each other contributes two mass ratios – one from each of the binaries.

Figure 12

The mass ratio distributions of binary systems with stellar primaries in the mass ranges M₁ > 0.5 M_⊙ (left-hand panels) and M₁= 0.1 − 0.5 M_⊙ (centre panels) and VLM primaries (right-hand panels; M₁ < 0.1 M_⊙). The top panels give the results at the end of the hydrodynamical calculation (0.3 Myr; Bate 2009a), while the results at the end of the N-body evolutions (10 Myr) with instantaneous gas removal and no gas removal are given in the middle and lower panels, respectively. The solid black lines give the observed mass ratio distributions of Duquennoy & Mayor (1991) for G-dwarfs (left-hand panels), Fischer & Marcy (1992) for M₁= 0.3–0.57 M_⊙ (centre, solid line) and M₁= 0.2–0.57 M_⊙ (centre, dashed line), and of the known VLM binary systems maintained by N. Siegler, C. Gelino and A. Burgasser at http://vlmbinaries.org/(right-hand panels). The observed mass ratio distributions have been scaled so that the areas under the distributions (M₂/M₁= 0.4–1.0 only for the centre panels) match those from the simulations. In all cases, VLM binaries are biased towards equal masses when compared with M-dwarf binaries (primary masses in the range M₁= 0.1–0.5 M_⊙). At 10 Myr with instantaneous gas removal, 75 per cent of the VLM binaries have M₂/M₁ > 0.6 while for the M-dwarf binaries the fraction is only 55 per cent. There is only weak evolution of the VLM and M-dwarf mass ratio distributions during the N-body evolution, even in the most extreme case of no gas removal. However, with fast gas removal, the N-body evolution results in a significant increase in the fraction of unequal mass solar-type binaries (the fraction of systems with mass ratios M₂/M₁ < 0.4 increases from 26 to 42 per cent in the case with instantaneous gas removal).

As noted by Bate, at the end of the hydrodynamical calculation, the ratio of near-equal mass systems to systems with dissimilar masses decreases going from VLM objects to M-dwarfs in a similar way to the observed mass ratio distributions, but that the trend is not as strong as in the observed systems. Specifically, 71 per cent of the VLM binaries had M₂/M₁ > 0.6 while for primary masses 0.1–0.5 M_⊙ the fraction was only 51 per cent. The M-dwarf mass ratio distribution is consistent with Fischer & Marcy's distribution. The VLM binaries, although biased towards equal-mass systems, are not as strongly biased as is observed. The main problem with the hydrodynamical binary mass ratio distributions, however, is that amongst ‘solar-type’ binaries with primary masses greater than 0.5 M_⊙ the proportion with mass ratios M₂/M₁ < 0.5 is much lower than that found by Duquennoy & Mayor.

At the end of the N-body evolution with instantaneous gas removal, the mass ratio distributions of the VLM and M-dwarf binaries have not changed significantly. 75 per cent of the VLM binaries have M₂/M₁ > 0.6 and 55 per cent of the M-dwarfs. However, there is a substantial change in the mass ratio distribution of solar-type binaries. The number of near-equal mass binaries has decreased somewhat, while the number of unequal mass binaries has substantially increased. Whereas only 9/34 = 26 per cent had M₂/M₁ < 0.4 at the end of the hydrodynamical calculation, the fraction is 15/36 = 42 per cent at the end of the N-body evolution.

As with the evolution of the separation distributions, the evolution of the mass ratio distributions discussed above is qualitatively similar for all the randomly perturbed realizations. The number of unequal mass solar-type binaries displays some variation. In the unperturbed case there were 15 with M₂/M₁ < 0.4 (out of 36), while the mean value is slightly lower at 10.7 ± 2.3, and the minimum and maximum values are 7 and 15, respectively. There is also some variation in the number of unequal mass VLM binaries. The unperturbed case left 9 with M₂/M₁ < 0.6 (out of 36), while the mean value is 8.7 ± 2.0 and the minimum and maximum values are 6 and 12, respectively.

The mass ratio distributions resulting from the N-body calculations with gas potentials also differ very little from those obtained with instantaneous gas removal. Even the extreme case of no gas removal differs only slightly from the instantaneous case (compare the middle and lower rows of Fig. 12).

To illustrate what has caused the increase in unequal mass solar-type binaries in the case with instantaneous gas removal, we discuss the formation of the three binaries with mass ratios 0.2 < M₂/M₁ < 0.4. One of these systems was formed when a triple consisting of a 2.5 au binary containing 0.96 and 0.73 M_⊙ stars with an 0.21 M_⊙ companion at 12 au broke up leaving the 0.73 and the 0.21 M_⊙ stars in a tight (0.1 au separation) binary with a low-mass ratio. Another was formed when a quadruple system consisting of two binaries with 21 au and 5 au semimajor axes orbiting each other with a 95 au semimajor axis broke up forming a new binary with a mass ratio of M₂/M₁= 0.31 consisting of one component from each of the original binaries in a 12 au semimajor axis. The last was formed when 1.74 and 0.40 M_⊙ stars that had been ejected from the cluster with similar velocities happened to find themselves weakly mutually bound in a 4100 au binary. As was seen in the previous section, a significant number of wide systems are formed in this manner and they tend to contain (although not exclusively) one of the more massive stars (>0.5 M_⊙), as might be expected for a wide pair of objects to be bound. Another low-mass ratio example consists of 0.70 and 0.10 M_⊙ stars in a 7500 au binary.

The formation of low-mass ratio solar-type binaries during the N-body evolution with instantaneous gas removal goes some way to addressing the problem that hydrodynamical simulations of star formation appear to underproduce unequal mass solar-type binaries (Delgado-Donate et al. 2004; Bate 2009a; Clarke 2009). Bate (2009a) suggested that if the mass ratios of triples and quadruples were included, this might also help to address the disagreement (i.e. essentially arguing that some of Duquennoy & Mayor's unequal mass wide binaries might have in fact been unidentified triples or quadruples). In Fig. 13, we plot the mass ratio distribution at the end of the N-body evolution of systems with primary masses >0.5 M_⊙ including the mass ratios between the wide components of triples and quadruple systems, but excluding VLM companions. For example, for a triple, the mass ratio is that between the close binary and the wider companion. For a quadruple consisting of two binaries, we include the mass ratio of the two binaries. We exclude VLM companions because the survey of Duquennoy & Mayor (1991) was not sensitive to VLM objects. The inclusion of triple and quadruple systems, in addition to the formation of unequal mass binaries during the N-body evolution, does not completely remove the discrepancy between the simulation and the survey of Duquennoy & Mayor (1991), but it does go a long way towards addressing it. We also note that the recent surveys of binaries with separations 5–5000 au in star-forming regions (Kraus et al. 2008; Kraus & Hillenbrand 2009) find uniform mass ratio distributions for primary masses ≳0.75 M_⊙, unlike that of Duquennoy & Mayor (1991). However, Kraus & Hillenbrand (2009) also found an indication of a variation between different star-forming regions with 9/12 of the wide binaries (500–5000 au) in Taurus having M₂/M₁ > 0.75 while in Upper Sco 6/11 of the binaries had M₂/M₁ < 0.25. In summary, the mass ratio distribution of binaries warrants further investigation, both observationally and theoretically.

Figure 13

The mass ratio distributions of binary, triple and quadruple systems with stellar primaries with masses M₁ > 0.5 M_⊙ at the end (10 Myr) of the unperturbed N-body evolution with instantaneous gas removal, excluding VLM companions (masses <0.1 M_⊙). The double hashed region gives the mass ratio distribution for binaries only. The single hashed region includes the mass ratios of triple and quadruple systems taking the mass ratio of the two most widely separated components. The solid black line gives the fit to the observed mass ratio distribution of Duquennoy & Mayor (1991) for G-dwarfs. The observed mass ratio distribution has been scaled so that the areas under the distribution matches that from the simulation.

4 DISCUSSION

For the first time, we have studied the evolution following gas dispersal of a stellar cluster produced from a hydrodynamical calculation of the collapse and fragmentation of a turbulent molecular cloud. Since only one hydrodynamical star formation calculation that forms a substantial cluster (≫100 stars) and resolves binary and multiple stellar systems has ever been performed (Bate 2009a) our results are, by necessity, limited to one particular type of cluster in terms of total mass, stellar density, etc. However, while it is difficult to extrapolate our results more widely (e.g. to more or less massive clusters), the calculations presented here at least illustrate the types of evolution that are likely to occur in other types of clusters.

4.1 Cluster evolution during gas dispersal

The cluster studied here results from the hydrodynamical collapse and fragmentation of a dense turbulent molecular cloud and, just before gas dispersal, it is very compact (half-mass radius 0.05 pc), with a very high central stellar density (∼10⁶ pc³), surrounded by a halo of ejected objects. Despite this initial state, upon removal of the gas (which makes up 62 per cent of the total mass), a bound remnant cluster, that contains about 30–40 per cent of the mass and number of stars, expands to a radius of ≈1–2 pc over a period of ≈4–10 Myr (depending on the gas removal time-scale). This evolution is consistent with the many past N-body studies of the effects gas expulsion on cluster evolution which have started from idealized initial conditions. For example, Baumgardt & Kroupa (2007) performed a wide variety of N-body simulations of the cluster dissolution following gas dispersal with varying star formation efficiencies. They found that their surviving clusters expanded by up to a factor of 10 in size. The obvious question raised by this sort of evolution is whether or not real stellar clusters may originate from such compact configurations? The amount of gas removed from the cluster we study here is large and only 30–40 per cent of the total number of stars remain in the remnant cluster. However, these numbers are not unreasonable. A star formation efficiency of ≈40 per cent has been inferred for the ONC (Hillenbrand & Hartmann 1998), while Weidner et al. (2007) find observationally that low-mass clusters keep at most 50 per cent of their stars during gas loss.

Recently, the observational evidence for significant expansion of young clusters has also been mounting. Scally, Clarke & McCaughrean (2005) showed that if the current estimates of the virial ratio of the ONC (which show that the cluster is unbound) are correct, the size and age of the ONC imply either that it became unbound only very recently, or else that it has expanded quasi-statically. In the latter case, they infer that its initial central density may have exceeded its current value of ∼10⁴ pc⁻³ (McCaughrean & Stauffer 1994; Hillenbrand & Hartmann 1998) by one to two orders of magnitude which would make it similar to the cluster studied here. Parker et al. (2009) also argue that the central region of the ONC was two orders of magnitude denser in the past with a half-mass radius of only 0.1–0.2 pc based on its current population of wide binaries. Moving to clusters other than the ONC, Bastian et al. (2008) studied several dozen Galactic and extragalactic stellar clusters and found evidence for rapid expansion of the cluster cores during the first 20 Myr of their evolution. Examining the evolution of cluster sizes with age, they found that clusters may begin with core radii ≲0.1 pc in size.

How else might such young cluster evolution be tested? From their N-body simulations, Baumgardt & Kroupa (2007) found that if gas expulsion is rapid, the stars acquire strongly radially anisotropic velocity dispersions outside their half-mass radii. Similarly, in the cluster evolution presented here, the ejected halo population which comprises >60 per cent of the stars and stellar mass are on unbound essentially radial trajectories with velocities that increase with distance from the cluster remnant (faster moving stars have travelled further).

Finally, as noted in Section 3.1.2, a common feature of the N-body calculations is for the most massive binary to be ejected from the remnant cluster, sometimes accompanied by a small group of ≈20 companion stars. Many Herbig Ae/Be stars are associated with small groups of stars (Hillenbrand 1994, 1995; Hillenbrand et al. 1995). Since the ejection of such a group occurs in the several of our randomly perturbed realizations, we suggest that this might be a significant formation mechanism of Herbig Ae/Be stars in small groups as opposed to an intermediate-mass star forming with a few other stars in an isolated molecular cloud.

4.2 The evolution of binary and multiple stellar systems

As noted in Section 1, many previous N-body simulations have investigated the evolution of binary properties, particularly binary frequency and separation distributions, during the early evolution of a young cluster (e.g. Kroupa 1995a,b, 1998; Parker et al. 2009). Many of the conclusions reached from these studies are also applicable to the cluster we have studied here. In particular, we find a cluster remnant surrounded by a halo of ejected stars with a lower binary fraction (Kroupa 1995b). We also find that some of the most massive stars can be ejected and reach large distances from the cluster within 10 Myr (Kroupa 1998). However, unlike Kroupa (1995a) and Parker et al. (2009), we do not find that many primordial binary systems are disrupted. Rather, the vast majority of binary systems formed during the star formation process (i.e. the hydrodynamical evolution) survive the gas dispersal phase with little evolution of binary frequency or their separation distribution. This is because stars do not form in the hydrodynamical simulations with 100 per cent of them being in binaries, nor do they form with a flat distribution in log-separation. Rather, star formation and dynamical interaction cannot be separated – they proceed simultaneously, resulting in a ‘primordial’ binary population that is already somewhat stable against further evolution. To try to separate star formation and the properties of primordial binary properties is likely to severely overestimate the importance of the dynamical evolution of the binary population early in the life of a stellar cluster.

By contrast, we note that higher order multiple systems (which are usually not included in purely N-body studies at all) do undergo significant evolution during the early evolution of a stellar cluster. In particular, many quadruple systems that result from the star formation process are found to decay into binaries or triples. Thus, we assert that observational studies of higher order multiple systems in very young stellar clusters are likely to find a much higher proportion of high-order multiples than is found in open clusters or the field population (e.g. Reipurth 2000).

Finally, as presented in Section 3.3, we find that a significant fraction of very wide (separations 10⁴–10⁵ au) multiple systems (around half of which are actually triples or quadruples) are produced in the halo of objects ejected from the cluster. This has also been noted in recent pure N-body simulations of dispersing clusters (Kouwenhoven et al. 2010). Thus, even though the hydrodynamical star formation process in a very dense cluster such as that discussed here (half-mass radius only 0.05 pc, or 10⁴ au at the end of the star formation) cannot produce very wide systems by direct fragmentation, such wide systems can be produced during cluster dissolution in numbers that are large enough to explain the observed field population of wide systems (Duquennoy & Mayor 1991). Such formation of wide systems during cluster dissolution potentially solves the problem (Parker et al. 2009) of how a broad log-separation distribution (e.g. Duquennoy & Mayor 1991) can be obtained if most stars form in clusters (Adams & Myers 2001; Lada & Lada 2003) which cannot contain wide binaries. Moreover, many of these wide binary systems and the binaries formed from the decay of higher order multiple systems have unequal mass components. This also helps to fix the deficit of unequal mass solar-type binaries (relative to the observations of Duquennoy & Mayor 1991) that is usually present at the end of the hydrodynamical star formation phase (Delgado-Donate et al. 2004; Bate 2009a; Clarke 2009). Thus, in summary, the formation of significant numbers of wide multiple systems during cluster dissolution potentially solves several problems simultaneously.

5 CONCLUSIONS

We have presented the results of N-body calculations that take the endpoint of the hydrodynamical star cluster formation calculation of Bate (2009a) and evolve the stellar cluster to an age of 10 Myr. The calculations are unique in the sense that they begin from a star cluster that has self-consistently formed from the hydrodynamical collapse and fragmentation of a molecular cloud to form a stellar cluster containing binary and higher order multiple systems rather than beginning from arbitrary initial conditions. Although we are limited to studying one particular type of cluster in terms of the number of stars, stellar density, etc, we find many of the same evolutionary processes that have been studied in past N-body calculations that begin from arbitrary initial conditions. For example, we find that a bound remnant cluster containing ≈30–40 per cent of the stars by number or mass remains despite the low star formation efficiency (38 per cent) and this cluster is surrounded by a halo of ejected stars. The remnant cluster expands from an initial radius of less than 0.05 pc to ≈1–2 pc over a period of ≈4–10 Myr (depending on the gas removal time-scale). This result, along with other recent numerical and observational studies, suggests that young clusters may begin from much higher stellar densities (10⁵–10⁷ pc⁻³) than usually assumed.

When investigating mass segregation in our young clusters, we differentiate between primordial mass segregation, which is a result of the star formation process, and classical mass segregation, which is due to dynamical relaxation. We find that primordial mass segregation is only evident for the few most massive stars (in agreement with observations of the ONC), but only lasts a short time (≈1 Myr). Later, classical mass segregation may develop with higher mass stars being increasingly centrally condensed, but this depends on the gas removal time-scale. Slower gas removal results in shorter relaxation times and, thus, more classical mass segregation.

Although many of the global cluster evolutionary phases have been discussed before in the literature, we also find some new evolutionary processes at work. We find that in the majority of our calculations, some of the most massive stars in the cluster are ejected from the bound cluster and are accompanied by a small group of companion stars. Thus, we propose that this mechanism might be responsible for the formation of Herbig Ae/Be stars that are accompanied by small stellar groups.

Turning to the evolution of binary and multiple systems, we note that past N-body studies of the evolution of young clusters that begin with high binary frequencies and wide separation distributions have emphasized the rapid destruction of binary systems. However, in our calculations, the binary frequency remains almost unchanged throughout the N-body evolution. We emphasize that the processes of star formation and the formation of multiple stellar systems cannot be separated during the formation of a cluster. Instead, they occur simultaneously, resulting in a ‘primordial’ binary population that is already somewhat stable against further evolution. That said, we do find that the higher order multiple systems produced by the star formation process, particularly quadruples, evolve significantly with most decaying into lower order systems. This is potentially observable by comparing the frequencies of triple and quadruple systems in young star-forming regions with those found in open clusters. We also find very wide (10⁴–10⁵ au) binary and higher order multiple systems are formed in the dispersing halo of ejected objects when objects happen to be expelled with similar velocities and find themselves weakly bound. These wide binaries are formed in sufficient numbers to explain the observed separation distribution of solar-type stars despite the fact that the original star-forming cloud is too small to form such systems by direct fragmentation. Finally, whereas the cluster at the end of the hydrodynamical evolution was deficient in unequal mass solar-type binaries compared with observations, we find that many of the wide binaries formed in the halo and those binaries resulting from the decay of higher order multiple systems have unequal mass components. Thus, these processes may help to bring the mass ratio distributions of binaries produced by hydrodynamical simulations into agreement with the observed mass ratio distributions of solar-type stars.

MRB is grateful for the support of a European Young Investigator (EURYI) Award and for the hospitality and financial support provided by the Isaac Newton Institute for Mathematical Sciences as part of its Dynamics of Discs and Planets programme where part of this paper was written up. This paper has made use of the Very Low Mass Binaries archive (http://vlmbinaries.org/) maintained by N. Siegler, C. Gelino and A. Burgasser. This work, conducted as part of the award ‘The Formation of Stars and Planets: Radiation Hydrodynamical and Magnetohydrodynamical Simulations’ made under the European Heads of Research Councils and European Science Foundation EURYI Awards scheme, was supported by funds from the Participating Organisations of EURYI and the EC Sixth Framework Programme.