Skip to main content Accessibility help


  • Access
  • Cited by 4


      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Pre-Mainsequence Stellar Evolution in N-Body Models
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Pre-Mainsequence Stellar Evolution in N-Body Models
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Pre-Mainsequence Stellar Evolution in N-Body Models
        Available formats
Export citation


We provide a set of analytic fits to the radii of pre-mainsequence stars in the mass range 0.1 < M/M < 8.0. We incorporate the formulae in N-body cluster models for evolution from the beginning of pre-main sequence. In models with 1 000 stars and high initial cluster densities, pre-mainsequence evolution causes roughly twice the number of collisions between stars than in similar models with evolution begun only from the zero-age main sequence. The collisions are often all part of a runaway sequence that creates one relatively massive star.


Existing N-body simulations of stellar clusters usually begin with all stars on the zero-age main sequence (ZAMS). Our motivation to include pre-mainsequence (preMS) evolution in such models goes beyond the desire for completeness. Given the upper mass limit (M>10 M ) attributed to conventional star formation (Stahler, Palla, & Ho 2000) and because massive stars are mainly found in the denser central regions of star clusters (e.g. Hillenbrand & Hartmann 1998), Bonnell, Bate, & Zinnecker (1998) first suggested collisions between preMS stars as a means to create massive stars. PreMS stars are larger than their main-sequence counterparts and, although the phase is short-lived, there is a possibility that the increased collision likelihood has an effect.

The journey from molecular gas cloud to protostar to the ZAMS is complicated. If random motions owing to MHD turbulence are sustained inside a molecular cloud, density inhomogeneities form and quickly collapse under their own gravity to become protostars (Stahler & Palla 2005). Typically, when modelling preMS evolution, this self-gravitating fragment of a protostellar cloud is where the models begin (Tout, Livio, & Bonnell 1999) and is where we define the zero age of the pre-main sequence.

In the 1950s, Henyey constructed the first detailed numerical models of young stars, assuming that preMS stars are radiatively stable (Henyey, Lelevier, & Levée 1955). This led to nearly horizontal tracks in the Hertzsprung–Russell (H–R) diagram because the stars contract at almost constant luminosity. Hayashi (1961) realised that the H opacity (and Kramer’s opacity, when the temperature increases) forces the young stars to radiate at nearly constant effective temperature T eff and so to follow nearly vertical tracks in an H–R diagram. These young stars are in a quasi-hydrostatic equilibrium. They contract on a time-scale much greater than the free-fall time-scale (Hayashi 1961). Although protostars are initially homogeneous and isothermal, they collapse non-homogeneously, creating high central density and temperature. It is this that eventually gives rise to the suitable conditions for hydrogen fusion in stars with M>0.08 M . Below this mass, the central conditions are never sufficient for hydrogen fusion and the stars collapse to become brown dwarfs (Kumar 1963; Hayashi & Nakano 1963).

Stars with M ≤ 0.7 M behave as though they are fully convective and reach the ZAMS at the end of their Hayashi track. When M ≥ 0.7 M , the star develops a radiative core which is large enough that it evolves more like Henyey’s initial models, turning on to the near-horizontal Henyey track where contraction continues at approximately constant luminosity. This behaviour can be seen in Figure 1, which shows the calculated preMS tracks modelled by Tout et al. (1999). However Stahler et al. (2000) argue that stars more massive than 10 M have their contraction disrupted by both strong winds and radiation pressure and thus the conventional theory of protostellar infall fails in this case. Zinnecker & Yorke (2007) discussed three competing mechanisms for massive star formation. These are the turbulent core model of McKee & Tan (2003), the competitive accretion process of Bonnell et al. (1997) and stellar collisions and merging of stars during their preMS evolution or later. It is this third case that we investigate here.

Figure 1. An H–R diagram showing the ZAMS in bold and the preMS tracks for a range of masses. The distinction between the vertical Hayashi track and horizontal Henyey track is apparent here. Data is taken from Tout et al. (1996) and Tout, Livio, & Bonnell (1999).

An alternative model was investigated by Baumgardt & Klessen (2011) and it is interesting to compare our results with theirs. They base their pre-mainsequence evolution on models of accreting stars constructed by Bernasconi & Maeder (1996). The physics used by Bernasconi & Maeder (1996) is not very different from that used by Tout et al. (1999) and indeed the evolutionary tracks in the H–R diagram are very similar once accretion has ceased and the stars contract down Hayashi tracks. The major difference between the models we present here and those of Baumgardt & Klessen (2011) is that all their stars begin as protostellar cores of 0.1M and accrete at a constant rate until they arrive at a suitable mass to populate an IMF up to 15M. Their stars then undergo a preMS phase in which the logarithm of their radii shrinks linearly with time until they reach the ZAMS, whereupon they begin main-sequence evolution. Instead we start with a set of pre-mainsequence stars with masses populating a similar IMF up to 4M that begin their evolution already fully grown on Hayashi tracks. Our fitting formulae are a little more complicated because we attempt to to fit both the Hayashi and Henyey phases of pre-mainsequence evolution together. Because of this our stars remain larger for longer (see Section 4) after any accretion has ceased. We include this behaviour in cluster models with a variety of initial conditions to test its effect.


We use the nbody6 1 (Aarseth 1999). This code incorporates stellar evolution and binary interaction by empirical formulae fitted to detailed stellar models. Two major additions to these standard packages are required. First we must have a reasonable empirical representation of how preMS stars, and particularly their radii, evolve. Secondly we must include a model of what happens when they collide.

2.1 Parameterization of the preMS evolution

We construct a set of fitting functions for the preMS evolution of stars with masses M in the range 0.1 < M/M < 8.0. We choose to find analytical fits rather than to tabulate data, as did Wiersma, Sills, & Portegeis Zwart (2006), in our models because this approach is in line with the parameterized treatment of stellar evolution in the other parts of the N-body code. Such fits also have the advantage that they are continuous and differentiable and this makes modelling preMS evolution with accretion simpler because the fits are smooth functions of mass. We began with detailed stellar models originally constructed by Tout et al. (1999) with the Cambridge stars code (Eggleton 1971) described by Pols et al. (1995). We used data for nine stars with solar metallicity Z = 0.02 and masses M∈{0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 8.0} M . These masses were chosen because stars with M < 0.08 M never meet the ZAMS and for M>8 M the time-scale of preMS evolution is short enough that it is safe to neglect it. The time taken for our 8 M star to contract from about 90R to the main sequence is about 3×105 yr. It takes only the first 1.25×104 yr to contract to 40R, a radius it does not exceed again until it becomes a red giant. It then spends about 5×104 yr around 36R while it burns its supply of deuterium fuel and then loses its convective envelope, so moving from the Hayashi to the Henyey phase. In Section 3.2 we compare this with the expected collision time-scale in the densest clusters that we model. We find that neglect of preMS evolution for stars above 8M is reasonable for these clusters but would not be for much denser systems.

Given the general difficulty of deciding the birth time of a star and hence its age, we measure the age of preMS stars backwards from the ZAMS. We therefore devise a slightly artificial preMS time-scale τpreMS by taking the time from which the Tout et al. (1999) models begin until their radius meets the ZAMS value. To do this we use the analytic ZAMS radius formula of Tout et al. (1996). It is a function of mass and metallicity but all our fits are made for solar metallicity Z = 0.02, for which differences from the detailed ZAMS models reported by Pols et al. (1995) are at most 1.2%. It is imperative that these particular errors are small because all other formulae calculate properties relative to these ZAMS values. Using a least squares method, we fit these times as a function of mass to get the preMS time-scale in the form

(1) \begin{eqnarray} &&\log _{10}\left(\frac{\tau _{\rm {preMS}}}{\rm {yr}}\right)=43.6-35.8\left(\frac{M}{M_{\odot }}\right)^{0.015}\nonumber\\ &&\quad\times \exp {\left[3.96 \times 10^{-3}\left(\frac{M}{M_{\odot }}\right)\right]} \end{eqnarray}
for our mass range. Next we define a scaled time
(2) \begin{equation} \tau =1-t/\tau _{\rm {preMS}} \end{equation}
such that time runs from τ = 1 at the beginning of the preMS track to the ZAMS at τ = 0 and hence τ is strictly in the range τ∈[0, 1]. These functions take the form
(3) \begin{equation} R=R_{\rm {ZAMS}}\left(M\right)10^{f(\tau )}, \end{equation}
(4) \begin{equation} f(\tau )=\frac{\alpha \tau ^3+\beta \tau ^4+\gamma \tau ^5}{1.05-\tau }. \end{equation}
For each of the nine models we obtain a best fit for α, β and γ and then fit these three coefficients as functions of mass. Several iterations were necessary to converge on good formulae. The coefficients α, β and γ are all well represented as piecewise cubics in mass of the forms
(5) \begin{eqnarray} \alpha^{(j)}\left(\frac{M}{M_{\odot }}\right)&=&\sum _{i=0}^3\alpha _i^{(j)}\left(\frac{M}{M_{\odot }}\right)^i, \nonumber \nonumber\\ \beta ^{(j)}\left(\frac{M}{M_{\odot }}\right)&=&\sum _{i=0}^3\beta _i^{(j)}\left(\frac{M}{M_{\odot }}\right)^i\nonumber\\ \text{and}\quad \gamma ^{(j)}\left(\frac{M}{M_{\odot }}\right)&=&\sum _{i=0}^3\gamma _i^{(j)}\left(\frac{M}{M_{\odot }}\right)^i, \end{eqnarray}
(6) \begin{equation} j= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 & M \le 1\,M_{\odot },\\ 2 & 1\,M_{\odot } < M < 2\,M_{\odot }, \\ 3 & M \ge 2\,M_{\odot } \end{array}\right. \end{equation}
and α i (j), β i (j) and γ i (j) are listed in Table 1. Our fits (illustrated in Figure 2) are all physically reasonable. They incorporate both the Hayashi and Henyey tracks. Sadly the fits are not as accurate as we may have liked. However we do not wish to make the formulae excessively complicated for what is a rather short phase of evolution. Our biggest deviations are around 5M, shown in the figure, for which we overestimate the radius by a factor of between about 1.25 and 1.5 for around 20% of its pre-mainsequence lifetime. Averaged over all masses we tend to underestimate as much as overestimate.

Figure 2. Fitted preMS radii for 1.0 M , 5.0 M and 8.0 M (bottom to top) at τ = 0.8. The solid lines are the detailed models and the dashed lines the fits. Note the reversed scale in τ, with τ = 1 at the beginning of the preMS track and τ = 0 when the star reaches the ZAMS.

Table 1. Coefficients for Equation (3).

2.2 Treatment of collisions

In addition to the treatment of stellar evolution we must also consider collisions of preMS stars with themselves and stars of other types. We model a collision between two preMS stars as if both are n = 3/2 polytropes because stars on Hayashi tracks are fully convective and most collisions occur whilst stars are largest. The gravitational energy of an n = 3/2 polytrope of mass M and radius R is

(7) \begin{equation} \Omega =-\frac{6}{7}\frac{GM^2}{R}. \end{equation}
The internal energy of a star in virial equilibrium is
(8) \begin{equation} U=-\frac{1}{3(\gamma -1)}\Omega . \end{equation}
This is a reasonable approximation to make because the preMS stars are in quasi-hydrostatic equilibrium. For a perfect gas, γ = 5/3, so
(9) \begin{equation} U=\frac{3}{7}\frac{GM^2}{R}. \end{equation}
We consider a collision between two preMS stars with masses M 1 and M 2, radii R 1 and R 2 and total mass M = M 1 + M 2. We assume that colliding stars merge while a small fraction of the total mass ξM is lost from the cluster (Davies et al. 1993), and let energy be conserved during the collision. The two stars are moving with high velocity when they collide but the kinetic energy at infinity, or at the apogee of an eccentric binary orbit, in their relative orbits can be neglected. The initial velocity of the colliding stars would need to be of order $10^3\,\rm {km}\,\rm {s}^{-1}$ to have kinetic energy comparable to the energy lost with ξM, whereas in globular clusters and star-forming regions the velocities are typically $10\,\rm {km}\,\rm {s}^{-1}$ (e.g. Portegies Zwart, McMillan, & Gieles 2010). Conserving the internal energy and ignoring this orbital contribution, we find that the radius R * of the new coalesced star is
(10) \begin{equation} R_{\ast }=\frac{7M^2}{3}\left(\frac{M_1^2}{R_1}+\frac{M_2^2}{R_2}\right)^{-1}(3-4\xi )(1-\xi ) \end{equation}
and its mass is
(11) \begin{equation} M_{\ast } = M_1 + M_2 - \xi M. \end{equation}
It is difficult to estimate mass loss in any stellar collision but, for MS stars, SPH calculations by Benz & Hills (1987) indicate about 10%. Thus ξ = 0.1. Laycock & Sills (2005) also found that a mass loss of a few per cent of the total mass per collision with a preMS star is a reasonable approximation. We invert the R = R(M *, τ) function (3) to find the rejuvenated age at which to restart evolution so that the star continues to contract from its new radius.

The outcome of a preMS/MS collision depends on the mass of the MS star. Below 0.7 M the MS star is still mostly convective and, to a good approximation, the collision remnant can be modelled as above. At higher masses, the MS star has very little convective envelope and a dense convective core appears above 1.1 M . Mass is therefore most likely to be accreted on to the surface of the star but not mixed as above and evolution would not restart on the preMS. All other collisions with different stellar types are treated simply as accretion, so the non-preMS stellar type is kept upon collision and the mass of the star is simply increased. Such prescription is again supported by Laycock & Sills (2005) who showed that the collision product, when the collision involves a preMS star, does not depend strongly on the impact parameter nor the initial velocity.


Our initial conditions are all Plummer (1911) models in virial equilibrium with no primordial binary stars. There is no interstellar gas nor mass segregation. Our stellar evolution is as described in Section 2.1 and by Hurley, Pols, & Tout (2000) and Hurley, Tout, & Pols (2002). We use the code nbody6 (Aarseth 1999). To explore a didactic set of N-body models, we keep the initial number of stars N, the initial mass function for the stars and the time-scale of evolution constant in our series of models. We vary only the half-mass radius R 0.5, the radius within which half the mass of the cluster is contained. We fix N to be 1 000. For the masses of the stars we use a modified Kroupa initial mass function (Kroupa 2001), with M min = 0.1M , M max = 4.0M and $\overline{M}=0.4\,M_{\odot }$ . We chose the length-scale $\overline{R}\in \lbrace 0.02,0.05,0.1,0.2,0.4\rbrace \rm {pc}$ . Our half-mass densities $\overline{\rho }_{0.5}$ range from 103 to 107 M pc − 3. The half-mass radius $R_{0.5}\approx 0.8\overline{R}$ .

Because our preMS fitting is restricted to the mass range 0.1 < M/M < 8.0 (Equation (3)), if (a) different IMF parameters were chosen with M max ⩾ 8.0M or (b) several stars collide to form a coalesced star with M * ≥ 8.0 M then these high-mass stars evolve from the ZAMS. This approximation is justified because the preMS evolution of stars of these high masses is so rapid that neglecting it loses little information (see Section 3.2). We evolved each cluster for t = 10Myr, approximately τpreMS for a 0.5 M star, with the logic that most interesting preMS behaviour would have occurred by this time. We made ten models for each of the length-scales, both with, all the stars started at the top of their Hayashi tracks at τ = 1, and without, all stars started on the ZAMS at τ = 0, preMS evolution. We shall hereinafter call these preMS and ZAMS runs respectively.

3.1 Densities

Our models have initial half-mass densities in the range $10^3 < \overline{\rho }_{0.5}/M_{\odot }\,{\rm pc}^{-3} < 10^7$ . The initial densities for the smallest clusters seem a little extreme but there is increasing evidence that the initial densities of open clusters are higher than previously thought (Parker et al. 2009), that bound clusters can expand quickly (Bastian et al. 2008; Moeckel & Bate 2010) and that rapid expansion can occur in the core (Kroupa, Aarseth, & Hurley 2001). Indeed Table 2 shows that the densest clusters expand the most because of the longer dynamical time required and, interestingly, end up with very similar half-mass radii to clusters that were initially somewhat sparser. The final densities vary because the initially denser clusters lose more mass throughout the model. There is thus some uncertainty in extrapolating back in time to estimate an initial cluster scaling.

Table 2. Evolution of the half-mass density $\overline{\rho }_{0.5}=M_{\rm {total}}/ \frac{4}{3}\pi R_{0.5}^3$ for $\overline{R}\in \lbrace {0.05,0.10,0.20}\,\rm {pc}\rbrace$ .

3.2 Collision frequency

First we looked at the number of collisions between stars for each $\overline{R}$ with and without preMS evolution. Averaged over ten simulations, the results are listed in Table 3, while Figure 3 shows the data and the averages over the 100 models. As expected, more stars collided in the initially denser clusters, as did more of the stars evolved from the top of their preMS tracks. The stars which collided in the ZAMS models had R≈2R whereas in the preMS runs the radii were at least twice this when they collided.

Figure 3. Number of collisions and their averages for 10 runs with and without preMS evolution for $\overline{R}=0.8\, R_{0.5}\in \lbrace 0.02,0.05,0.1,0.2,0.4\rbrace \,\text{pc}$ .

Table 3. Collision statistics for the models.

The gravitational focusing cross-section for stars with masses M 1 and M 2, and radii R 1 and R 2 is

(12) \begin{equation} \sigma (M_1,M_2)=\pi (R_1+R_2)^2\left(1+\frac{2G(M_1+M_2)}{v_{\rm {rel}}^2(R_1+R_2)}\right), \end{equation}
where v rel is the relative velocity of the two stars at infinity or at apocentre. Because v rel is usually small, this reduces to
\begin{equation*} \sigma (M_1,M_2)\propto (R_1+R_2). \end{equation*}
To determine whether we are justified in ignoring the pre-mainsequence evolution of stars over 8M we check the collision rate for such stars in our models. For our densest clusters with a half-mass radius of 0.02 pc the stellar number density n ≈ 9 × 10 − 16R − 3 and the typical relative velocity v rel ≈ 11.5kms − 1 ≈ 520Ryr − 1. We estimate the collision cross-section by setting R 1 + R 2 = 100R and M 1 + M 2 = 10M and so find σ ≈ 9 × 106R2 . From this we estimate a collision time-scale
(13) \begin{equation} t_{\rm coll} \approx \frac{1}{n\sigma v_{\rm rel}} \approx 2.4\times 10^5\,\rm yr. \end{equation}
This is just a little shorter than the entire pre-mainsequence lifetime of 3×105 yr for our 8M. We recall that it actually spends only 1.25×104 yr above 40R. So we can expect stars more massive than 8M to contract to the main sequence before they collide again, but only just in our densest cluster models. If we were to model higher densities we would need to model the pre-mainsequence evolution of more massive stars.

3.3 Runaways

In nearly every case, these collisions involved one or two stars colliding many times. We call a star that collides with another more than once a runaway. These were usually the most massive stars in the initial configuration, with initial masses 3 ⩽ M initial/M ⩽ 4. Some gained over ten times their initial mass within 10Myr through multiple collisions. Table 4 is a list of the ten most massive runaways in both preMS and ZAMS models. This shows that not only do more collisions occur in the preMS models but also that the runaways end up somewhat more massive. The mean mass gained by runaways in the preMS models was 12.0 M (averaged over 34 runaways) and 8.8 M (24 runaways) in the ZAMS models.

Table 4. The ten most massive runaways.

1 Time is the age of the cluster when the runaway star collided for the last time.

2 This star later evolved to a 6.8 M black hole by 8.5 Myr and absorbed another preMS star.

Portegies Zwart & McMillan (2002) found a similar phenomenon in their rather extreme models of colliding MS stars in dense star clusters where MS/MS collisions typically form one runaway star. In Monte Carlo stellar dynamical models Freitag, Gürkan, & Rasio (2006), while investigating a runaway mechanism to create an intermediate-mass black hole in compact stellar clusters, also found that only one runaway object was formed. Our results are in line with these findings. Thus multiple collisions between preMS stars early on in cluster evolution can be a viable mechanism to create a few massive stars in otherwise low-mass clusters but is unlikely to be the means to populate the top of the IMF (Zinnecker & Yorke 2007). This conclusion was reached by Baumgardt & Klessen (2011) who also found typically only one runaway massive star in their simulated clusters.

We found that the characteristic time to the first runaway is about the same as that to form a small core, 0.05pc in radius. The collisions, or merges, then occur within this core. They also usually take place at the pericentre of highly eccentric binary orbits. In a few cases a fly-by was seen to induce collisions between two stars in an eccentric binary system.

Usually the most massive runaway stars formed in the very dense $\overline{R}=0.02\,\rm {pc}$ model. For $\overline{R}=0.4$ no runaway stars formed at all and only three formed in the preMS model when $\overline{R}=0.2$ . However there is significant variation between the models. In the two least dense models the effect of the preMS evolution is insignificant because the clusters are sparse enough that by the time two stars come sufficiently close to collide the preMS phase is over.


We have not attempted to create a completely realistic model of a cluster. The background gas that would exist in a protostellar cloud has been neglected, as have the effects of accretion of this gas (see the treatment by Baumgardt & Klessen 2011). Nor have we modelled the effects of stellar accretion discs which would enhance the likelihood of collisions and also change their nature. Our aim here has been to show that preMS evolution can have important consequences and to identify when the preMS phase should be considered in more detail in the future.

It is not straightforward to make a direct comparison with Baumgardt & Klessen (2011) because they let their stars form by accretion on to 0.1 M cores, at constant accretion rates to populate the IMF, whereas we start with an IMF and evolve our protostars as if they were coeval. However, examining the time-integrated cross-section of our analytic fits compared to theirs (see Figure 4) shows that their models systematically underestimate the radii of the preMS stars below 6 M compared to ours. This is important when looking at collision frequencies. Calculating the time and mass-integrated cross-section

(14) \begin{equation} I_{\rm {IMF}}=\int _{0.1\,M_\odot }^{4.0\,M_\odot } \int _{\tau _{\rm preMS}} \eta (M) R(M,\tau )\,\mathrm{d}\tau \,\mathrm{d}M, \end{equation}
where η(M) is the Kroupa IMF normalised so that
(15) \begin{equation} \int _{0.1\,M_\odot }^{4.0\,M_\odot }\eta (M)\,\mathrm{d}M = 1, \end{equation}
we find that for the original models I IMF = 11.5, for our analytic fits I IMF = 12.9 and for the fits of Baumgardt & Klessen (2011), I IMF = 7.9. Though we should therefore expect to find more collisions in the preMS regime than they do the differences in the models make it almost impossible to identify whether this is actually the case.

Figure 4. The integrated cross-section I = ∫τpreMS Rdτ, taken over the preMS lifetime given by Equation (1). The models here are those of Tout et al. (1999). The fits of Baumgardt & Klessen (2011) are integrated from the time their accretion phase stops to when their models reached the ZAMS. Our cross-sections are relatively larger for the longer lived lower-mass stars but smaller between about 5.5 and 7.5M.

The inclusion of gas is a natural next step in these models. We would expect an interplay between gas expulsion, which decreases the star formation efficiency, accretion rates and the effect of gas accretion on to stars. Gas accretion not only leads to larger stars and hence larger collision cross-sections but, in clusters, forces contraction of the whole system and thence leads to more merging stars (Bonnell et al. 1998; Bonnell & Bate 2002), an increase in the number of binary stars and early mass segregation (Moeckel & Clarke 2011).

We would also expect our preMS stars to form discs while accreting or indeed as a result of tidal disruption during a collision (Davies et al. 2006). The effect of such protoplanetary discs on the number of collisions, their outcome and the time-scale on which they take place should also be modelled in some way and included in a more detailed N-body model. Star–disc interactions can lead to the disc being stripped away (Hall, Clarke, & Pringle 1996; Scally & Clarke 2001) or the formation of binary stars (Clarke & Pringle 1993) and, because many massive stars live in tight binary systems (Zinnecker & Bate 2002), it will be important to know whether the preMS phase contributes to this.

In conclusion, although there is much physics that we have not included, our models show that preMS evolution increases both the number of collisions, when the density is sufficiently large, and the amount of mass gained by the final collision remnants. Most collisions are part of a runaway process so preMS collisions alone are probably insufficient to populate the upper IMF. We identify an initial half-mass cluster density of 104 M pc − 3 below which preMS evolution can be neglected because collisions are sufficiently rare. However caution should be exercised owing to large uncertainty in the early length scales of clusters.


CAT is very grateful to Churchill College for his Fellowship. We thank the referee for his helpful suggestions for improvement.

1 This code is available to download at


Aarseth, S. J. 1999, PASP, 111, 1333
Bastian, N., Gieles, M., Goodwin, S. P., Trancho, G., Smith, L. J., Konstantopoulos, I., & Efremov, Y. 2008, MNRAS, 389, 223
Baumgardt, H., & Klessen, R. S. 2011, MNRAS, 413, 1810
Benz, W., & Hills, J. G. 1987, ApJ, 323, 614
Bernasconi, P. A., & Maeder, A. 1996, A&A, 307, 829
Bonnell, I. A., & Bate, M. R. 2002, MNRAS, 336, 659
Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 1997, MNRAS, 285, 201
Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS, 298, 93
Clarke, C. J., & Pringle, J. E. 1993, MNRAS, 261, 190
Davies, M. B., Bate, M. R., Bonnell, I. A., Bailey, V. C., & Tout, C. A. 2006, MNRAS, 370, 2038
Davies, M. B., Ruffert, M., Benz, W., & Muller, E. 1993, A&A, 272, 430
Eggleton, P. P. 1971, MNRAS, 151, 351
Freitag, M., Gürkan, M. A., & Rasio, F. A. 2006, MNRAS, 368, 141
Hall, S. M., Clarke, C. J., & Pringle, J. E. 1996, MNRAS, 278, 303
Hayashi, C. 1961 PASJ, 13, 450
Hayashi, C., & Nakano, T. 1963, Progress of Theoretical Physics, 30, 460
Henyey, L. G., Lelevier, R., & Levée, R. D. 1955, PASP, 67, 154
Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540
Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
Kroupa, P. 2001, MNRAS, 322, 231
Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699
Kumar, S. S. 1963, ApJ, 137, 1121
Laycock, D., & Sills, A. 2005, ApJ, 627, 277
McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850
Moeckel, N., & Bate, M. R. 2010, MNRAS, 404, 721
Moeckel, N., & Clarke, C. J. 2011, MNRAS, 410, 2799
Parker, R. J., Goodwin, S. P., Kroupa, P., & Kouwenhoven, M. B. N. 2009, MNRAS, 397, 1577
Plummer, H. C. 1911, MNRAS, 71, 460
Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964
Portegies Zwart, S. F., & McMillan, S. L. W. 2002, ApJ, 576, 899
Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431
Scally, A., & Clarke, C. 2001, MNRAS, 325, 449
Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Weinheim: Wiley-VCH)
Stahler, S. W., Palla, F., & Ho, P. T. P. 2000, in Protostars and Planets IV, ed. Mannings, V., Boss, A. P., & Russell, S. S. (Tucson: Univ. Arizona Press), 327
Tout, C. A., Livio, M., & Bonnell, I. A. 1999, MNRAS, 310, 360
Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257
Wiersma, R., Sills, A., & Portegies Zwart, S. 2006, ApJ, 637, 838
Zinnecker, H., & Bate, M. R. 2002, in ASP Conf. Ser. 267, Hot Star Workshop III: The Earliest Phases of Massive Star Birth, ed. Crowther, P. (San Francisco: ASP), 209
Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481