## Abstract

We present a new study of the collisional runaway scenario to form an intermediate-mass black hole (IMBH, MBH≳ 100 M) at the centre of a young, compact stellar cluster. The first phase is the formation of a very dense central core of massive stars (M*≃ 30–120 M) through mass segregation and gravothermal collapse. Previous work established the conditions for this to happen before the massive stars evolve off the main sequence (MS). In this and a companion paper, we investigate the next stage by implementing direct collisions between stars. Using a Monte Carlo stellar dynamics code, we follow the core collapse and subsequent collisional phase in more than 100 models with varying cluster mass, size, and initial concentration. Collisions are treated either as ideal, ‘sticky-sphere’ mergers or using realistic prescriptions derived from 3D hydrodynamics computations. In all cases for which the core collapse happens in less than the MS lifetime of massive stars (≃3 Myr), we obtain the growth of a single very massive star (VMS, M*≃ 400–4000 M) through a runaway sequence of mergers. Mass loss from collisions, even for velocity dispersions as high as σv∼ 1000 km s−1, does not prevent the runaway. The region of cluster parameter space leading to runaway is even more extended than predicted in previous work because, in clusters with σv > 300 km s−1, collisions accelerate (and, in extreme cases, drive) core collapse. Although the VMS grows rapidly to ≳1000 M in models exhibiting runaway, we cannot predict accurately its final mass. This is because the termination of the runaway process must eventually be determined by a complex interplay between stellar dynamics, hydrodynamics, and the stellar evolution of the VMS. In the vast majority of cases, we find that the time between successive collisions becomes much shorter than the thermal time-scale of the VMS. Therefore, our assumption that all stars return quickly to the MS after a collision must eventually break down for the runaway product, and the stellar evolution of the VMS becomes very uncertain. For the same reason, the final fate of the VMS, including its possible collapse to an IMBH, remains unclear.

## 1 INTRODUCTION

In our first paper (Freitag, Rasio & Baumgardt 2006; hereafter Paper I), we have presented our numerical approach to the study of stellar collisions in young, dense star clusters with a broad stellar mass spectrum. It is based on the use of me(ssy)**2 (for ‘Monte Carlo Experiments with Spherically SYmmetric Stellar SYstems’), a Monte Carlo (MC) code to simulate the long-term evolution of spherical clusters subject to relaxation, collisions, stellar evolution and taking into account the possible presence of a central massive object (Freitag & Benz 2001, 2002). Our main motivation is to investigate in detail the runaway growth of a massive object (‘very massive star’, hereafter VMS) during core collapse (Ebisuzaki et al. 2001; Portegies Zwart & McMillan 2002; Gürkan, Freitag & Rasio 2004, hereafter GFR04). This provides a natural path to the formation of an intermediate-mass black hole (IMBH) in a dense star cluster, or a seed black hole (BH) in a proto-galactic nucleus In Paper I, we presented a number of test calculations to validate our MC code and compare its results for simple idealized systems to those from N-body, Fokker–Planck (FP) and gaseous-model codes. These comparisons established that we can reliably follow the relaxation-driven evolution of clusters with a broad mass function, all the way to core collapse. We also found good agreement with models of dense clusters in which mergers between stars create a mass spectrum starting from a single-mass population, thus accelerating collapse and leading to collisional runaway (Quinlan & Shapiro 1990).

Here we use me(ssy)**2 to perform a large number of calculations covering broadly the parameter space of young star clusters. Specifically, we vary systematically the number of stars (in the range 105–108, using up to 9 × 106 particles) and the size of the cluster (with half-mass radius values from 0.02 to 5 pc), and we consider systems with low and high initial concentration (King parameter W0= 3 and W0= 8). Our work is guided by the main finding of GFR04, namely that for clusters with realistic initial mass functions (IMFs) (for which the ratio of the maximum stellar mass to the average mass is larger than ∼100), mass segregation drives the cluster to core collapse in a time not longer than 15 per cent of the initial central relaxation time. It is therefore expected that in any cluster for which this relaxational core-collapse time is shorter that the main-sequence (MS) lifetime of the most massive stars (M*≈ 100 M), that is, ∼3 Myr, these objects will eventually collide with each other in the high-density core, likely in a runaway fashion. As we show in Section 2, our results confirm this expectation.

Our work on the subject was inspired by the investigations of Portegies Zwart and collaborators who revived the study of the runaway scenario thanks to N-body simulations (Portegies Zwart et al. 1999; Portegies Zwart & McMillan 2002). Although many key aspects of the process were already known from older works (in particular Spitzer & Saslaw 1966; Colgate 1967; Sanders 1970; Begelman & Rees 1978; Vishniac 1978; Lee 1987; Quinlan & Shapiro 1990; see Paper I, for a review of the field), these recent N-body studies were seminal in considering the effects of collisions in a cluster with a realistically broad IMF. This demonstrated explicitly for the first time how mass segregation can lead to a collisional phase by concentrating massive stars in a small central volume and made it clear that the post-MS evolution of the most massive stars initially present can prevent the runaway phase by driving cluster expansion. Furthermore, Portegies Zwart and collaborators showed that, contrary to the somewhat unrealistic situation in single-mass clusters (Lee 1987; Quinlan & Shapiro 1990), dynamically formed binaries, far from preventing the runaway phase by halting the central concentration increase, promote it, as many binary interactions lead to stellar mergers. As noted in Paper I and suggested by N-body simulations with higher particle numbers (Portegies Zwart et al. 2004), in a broad-IMF cluster with N*≳ 106, binaries probably cannot form through three-body interactions before the collisional phase is reached.

It is important to stress, however, that our MC simulations probe a different regime from that explored by the direct N-body approach. We consider systems with higher number of stars in the central region (inside a few core radii). For this reason, some of our most important results are only superficially similar to the findings of Portegies Zwart et al. In particular, in such large-N clusters with a long initial central collision time, a genuine core collapse, itself driven by mass segregation, appears to be the condition for collisional runaway. The situation for systems with a smaller number of stars is more complex because there may not be enough massive stars in the central regions to drive such a core collapse (Portegies Zwart et al. 1999; Portegies Zwart & McMillan 2002; McMillan & Portegies Zwart 2004; Portegies Zwart et al. 2004).

Our paper is organized as follows. In Section 2, we present the simulations we have performed and explain their results. In Section 3, we summarize our findings and discuss avenues for future research in the field. Detailed presentations of the runaway scenario, the physics at play, our numerical methods and test computations are to be found in Paper I.

## 2 SIMULATIONS

### 2.1 Initial conditions and units

In this work, we consider the evolution of isolated spherical stellar cluster with a broad IMF, subject to two-body relaxation and stellar collisions. When not stated otherwise, we use a Salpeter IMF, for which the number, dN*, of stars with masses between M* and M*+ dM* is given by

1

with α= 2.35, Mmin= 0.2 M and Mmax= 120 M. For this stellar population, the average mass is 〈M*〉≃ 0.69 M. There is no initial mass segregation. To investigate the role of the initial concentration of the cluster, we consider King models with W0= 3 or 8 (Binney & Tremaine 1987, Section 4.4). We do not enforce tidal truncation because it was shown in GFR04 that this does not affect the central regions which completely dominate the evolution.

We refer to GFR04 (Section 3 and Table 1) for detailed explanations about the important physical parameters of such clusters and units. In keep with the tradition, when not stated otherwise, we are using the N-body unit system (Hénon 1971) defined by G= 1, Mcl(0) = 1 (initial total cluster mass) and Ucl (0) =−1/2 (initial cluster potential energy). As time-unit, we prefer the ‘FP’ time TFP to the N-body unit TNB because the former is a relaxation time while the latter is a dynamical time; they are related to each other by TFP=N*/ln (γcN*) TNB where N* is the total number of stars in the cluster. As explained in Paper I, we conservatively use γc= 0.01 to determine relaxation times in this work. For King models, the N-body length unit is close to the half-mass radius, RNB≃ 1.19 Rh for W0= 3 and RNB≃ 1.15 Rh for W0= 8. The core radius (Spitzer 1987, equations 1–34 and Paper I) is Rcore≃ 0.543 RNB and encompasses a fraction 0.238 of the total mass for W0= 3; for W0= 8, these values are 0.121 RNB and 0.0531, respectively.

Table 1

Properties of Simulated Clusters.

Name W0 N* Np RNB (pc) Coll. Peculiarities tra (TFPtra (Myr)
K3-1 105 105 0.3 SS 1–120 M 4.38 × 10−2 2.77
K3-2 105 105 0.4 SS … 1.08 × 10−2 2.25
K3-3 105 105 0.4 SS 1–120 M … …
K3-4 105 105 0.5 SS … … …
K3-5 105 105 0.5 SS … 1.05 × 10−2 3.03
K3-6 105 105 0.5 SS … 1.07 × 10−2 3.08
K3-7 105 105 0.5 SS … 1.02 × 10−2 2.94
K3-8 105 105 0.5 SS 1–120 M … …
K3-9 3 × 105 3 × 105 0.03 SS … 9.89 × 10−3 0.0634
K3-10 3 × 105 3 × 105 0.1 SS … 1.06 × 10−2 0.414
K3-11 3 × 105 3 × 105 0.1 SPH … 1.04 × 10−2 0.406
K3-12 3 × 105 3 × 105 0.2 SS … 1.05 × 10−2 1.16
K3-13 3 × 105 3 × 105 0.3 SS … 9.60 × 10−3 1.95
K3-14 3 × 105 3 × 105 0.3 SS Constant 5 per cent mass loss 1.03 × 10−2 2.09
K3-15 3 × 105 3 × 105 0.3 SS Constant 10 per cent mass loss 1.03 × 10−2 2.09
K3-16 3 × 105 3 × 105 0.3 SS Constant 30 per cent mass loss … …
K3-17 3 × 105 3 × 105 0.3 SS Fixed central object 1.01 × 10−2 2.05
K3-18 3 × 105 3 × 105 0.3 5 per cent mass loss 9.35 × 10−3 1.9
K3-19 3 × 105 3 × 105 0.3 SS Constant R* for VMS 9.66 × 10−3 1.96
K3-20 3 × 105 3 × 105 0.3 SS 0.2–10 M 9.04 × 10−2 19.9
K3-21 3 × 105 3 × 105 0.3 SS Kroupa 0.08–120 M 1.02 × 10−2 2.54
K3-22 3 × 105 3 × 105 0.3 SPH … 8.88 × 10−3 1.8
K3-23 3 × 105 3 × 105 0.4 SS … 9.50 × 10−3 2.97
K3-24 3 × 105 3 × 105 0.4 SPH … 9.28 × 10−3 2.9
K3-25 3 × 105 3 × 105 0.4 SS … 9.69 × 10−3 3.02
K3-26 3 × 105 3 × 105 0.4 SS … 1.02 × 10−2 3.18
K3-27 3 × 105 3 × 105 0.4 SS Check for orbit overlap 9.22 × 10−3 2.85
K3-28 3 × 105 3 × 105 0.5 SS … … …
K3-29 3 × 105 3 × 105 0.5 SS 0.2–10 M … …
K3-30 106 106 0.1 SS Kroupa 0.01–120 M 3.50 × 10−3 0.327
K3-31 106 106 0.2 SS … 7.07 × 10−3 1.87
K3-32 106 106 0.2 SPH … 8.00 × 10−3 1.39
K3-33 106 106 0.3 SS … 7.94 × 10−3 2.54
K3-34 106 106 0.4 SS … … …
K3-35 3 × 106 3 × 105 0.3 SS … … …
K3-36 3 × 106 3 × 106 0.1 SS … 6.44 × 10−3 0.617
K3-37 3 × 106 3 × 106 0.2 SS … 6.84 × 10−3 1.85
K3-38 3 × 106 3 × 106 0.2 SPH … 7.25 × 10−3 1.97
K3-39 3 × 106 3 × 106 0.3 SS … … …
K3-40 107 3 × 105 0.1 SS … 7.75 × 10−3 1.21
K3-41 107 3 × 105 0.1 SS … 7.53 × 10−3 1.18
K3-42 107 3 × 105 0.2 SS … 8.79 × 10−3 3.89
K3-43 107 3 × 105 0.25 SS … … …
K3-44 107 3 × 105 0.3 SS … … …
K3-45 107 3 × 105 0.4 SS … … …
K3-46 107 106 0.1 SPH … 6.41 × 10−3 0.998
K3-47 107 106 0.2 SS … 6.27 × 10−3 2.76
K3-48 107 106 0.25 SS … … …
K3-49 3 × 107 3 × 105 0.1 SS … 5.24 × 10−3 1.3
K3-50 108 3 × 105 0.03 SS … 3.49 × 10−4 0.0237
K3-51 108 3 × 105 0.04 SPH … 2.99 × 10−3 0.312
K3-52 108 3 × 105 0.05 SS … 7.36 × 10−4 0.107
K3-53 108 3 × 105 0.05 SS No relaxation 8.90 × 10−4 0.13
K3-54 108 3 × 105 0.1 SS … 1.94 × 10−3 0.801
K3-55 108 3 × 105 0.2 SS … 3.41 × 10−3 3.98
K3-56 108 3 × 105 0.2 SPH … … …
K3-57 108 3 × 105 0.25 SS … … …
K3-58 108 3 × 105 0.3 SS … … …
K3-59 108 3 × 105 0.4 SS … … …
K3-60 108 3 × 105 0.5 SS … … …
K3-61 108 106 0.04 SS … 5.05 × 10−4 0.0524
K3-62 108 106 0.1 SPH … 3.19 × 10−3 1.31
K3-63 108 106 0.1 SPH … 3.20 × 10−3 1.31
K3-64 108 106 0.1 SS … 1.90 × 10−3 0.78
K3-65 108 106 0.2 SPH … … …
K8-1 3 × 105 3 × 105 0.9 SS … 3.97 × 10−4 0.418
K8-2 3 × 105 3 × 105 SS … 4.63 × 10−4 1.62
K8-3 3 × 105 3 × 105 SS … … …
K8-4 3 × 105 3 × 105 SS … … …
K8-5 3 × 105 3 × 106 1.2 SPH … 3.89 × 10−4 0.618
K8-6 3 × 105 3 × 106 4.7 SPH … … …
K8-7 106 106 0.1 SS … 4.00 × 10−4 0.0246
K8-8 106 106 SS … 4.32 × 10−4 0.841
K8-9 106 106 SS … 4.82 × 10−4 2.65
K8-10 106 106 SS … … …
K8-11 1.4 × 106 1.4 × 106 0.96 SS Fixed central object 4.02 × 10−4 0.843
K8-12 1.4 × 106 1.4 × 106 0.96 SS Fixed central object 3.91 × 10−4 0.819
K8-13 1.4 × 106 1.4 × 106 0.96 SS Fixed central object, tidal truncation 3.75 × 10−4 0.786
K8-14 1.4 × 106 1.4 × 106 0.96 SS … 4.66 × 10−4 0.976
K8-15 1.4 × 106 1.4 × 106 0.96 SS … 3.98 × 10−4 0.834
K8-16 1.4 × 106 1.4 × 106 0.96 SS Fixed central object, stellar evolution 4.19 × 10−4 0.878
K8-17 1.4 × 106 1.4 × 106 0.96 SS Fixed central object, stellar evolution 3.91 × 10−4 0.819
K8-18 1.4 × 106 1.4 × 106 0.96 SS … 4.32 × 10−4 0.904
K8-19 3 × 106 3 × 106 0.2 SPH … 3.31 × 10−4 0.0892
K8-20 3 × 106 3 × 106 SS … 3.43 × 10−4 1.03
K8-21 3 × 106 3 × 106 SPH … 3.43 × 10−4 1.03
K8-22 3 × 106 3 × 106 SS … 3.70 × 10−4 1.12
K8-23 3 × 106 3 × 106 SS … 3.74 × 10−4 1.13
K8-24 3 × 106 3 × 106 SS Fixed central object 2.99 × 10−4 0.9
K8-25 3 × 106 3 × 106 1.2 SPH … 4.22 × 10−4 1.65
K8-26 3 × 106 3 × 106 1.4 SS … 4.27 × 10−4 2.13
K8-27 3 × 106 3 × 106 SS … … …
K8-28 3 × 106 3 × 106 4.7 SPH … … …
K8-29 9 × 106 9 × 106 SS … 3.78 × 10−4 1.78
K8-30 107 3 × 106 SS … 3.61 × 10−4 1.78
K8-31 107 3 × 106 1.4 SS … 3.12 × 10−4 2.55
K8-32 107 3 × 106 1.4 SS … 3.86 × 10−4 3.15
K8-33 107 3 × 106 SS … … …
K8-34 3 × 107 3 × 106 SS … 3.14 × 10−4 2.45
K8-35 3 × 107 3 × 106 1.4 SS … … …
K8-36 108 3 × 106 0.1 SS … 7.49 × 10−5 0.0308
K8-37 108 3 × 106 0.6 SS … 2.40 × 10−4 1.45
K8-38 108 3 × 106 0.8 SS … 2.74 × 10−4 2.55
K8-39 108 3 × 106 SS … … …
K8-40 108 3 × 106 1.4 SS … … …
K8-41 108 3 × 106 SS … … …
MGG9-K8 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 … …
MGG9-K9 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 1.55 × 10−4 1.95
MGG9-K12a 12 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 5.10 × 10−6 0.0644
MGG9-K12b 12 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 5.98 × 10−6 0.0755
Name W0 N* Np RNB (pc) Coll. Peculiarities tra (TFPtra (Myr)
K3-1 105 105 0.3 SS 1–120 M 4.38 × 10−2 2.77
K3-2 105 105 0.4 SS … 1.08 × 10−2 2.25
K3-3 105 105 0.4 SS 1–120 M … …
K3-4 105 105 0.5 SS … … …
K3-5 105 105 0.5 SS … 1.05 × 10−2 3.03
K3-6 105 105 0.5 SS … 1.07 × 10−2 3.08
K3-7 105 105 0.5 SS … 1.02 × 10−2 2.94
K3-8 105 105 0.5 SS 1–120 M … …
K3-9 3 × 105 3 × 105 0.03 SS … 9.89 × 10−3 0.0634
K3-10 3 × 105 3 × 105 0.1 SS … 1.06 × 10−2 0.414
K3-11 3 × 105 3 × 105 0.1 SPH … 1.04 × 10−2 0.406
K3-12 3 × 105 3 × 105 0.2 SS … 1.05 × 10−2 1.16
K3-13 3 × 105 3 × 105 0.3 SS … 9.60 × 10−3 1.95
K3-14 3 × 105 3 × 105 0.3 SS Constant 5 per cent mass loss 1.03 × 10−2 2.09
K3-15 3 × 105 3 × 105 0.3 SS Constant 10 per cent mass loss 1.03 × 10−2 2.09
K3-16 3 × 105 3 × 105 0.3 SS Constant 30 per cent mass loss … …
K3-17 3 × 105 3 × 105 0.3 SS Fixed central object 1.01 × 10−2 2.05
K3-18 3 × 105 3 × 105 0.3 5 per cent mass loss 9.35 × 10−3 1.9
K3-19 3 × 105 3 × 105 0.3 SS Constant R* for VMS 9.66 × 10−3 1.96
K3-20 3 × 105 3 × 105 0.3 SS 0.2–10 M 9.04 × 10−2 19.9
K3-21 3 × 105 3 × 105 0.3 SS Kroupa 0.08–120 M 1.02 × 10−2 2.54
K3-22 3 × 105 3 × 105 0.3 SPH … 8.88 × 10−3 1.8
K3-23 3 × 105 3 × 105 0.4 SS … 9.50 × 10−3 2.97
K3-24 3 × 105 3 × 105 0.4 SPH … 9.28 × 10−3 2.9
K3-25 3 × 105 3 × 105 0.4 SS … 9.69 × 10−3 3.02
K3-26 3 × 105 3 × 105 0.4 SS … 1.02 × 10−2 3.18
K3-27 3 × 105 3 × 105 0.4 SS Check for orbit overlap 9.22 × 10−3 2.85
K3-28 3 × 105 3 × 105 0.5 SS … … …
K3-29 3 × 105 3 × 105 0.5 SS 0.2–10 M … …
K3-30 106 106 0.1 SS Kroupa 0.01–120 M 3.50 × 10−3 0.327
K3-31 106 106 0.2 SS … 7.07 × 10−3 1.87
K3-32 106 106 0.2 SPH … 8.00 × 10−3 1.39
K3-33 106 106 0.3 SS … 7.94 × 10−3 2.54
K3-34 106 106 0.4 SS … … …
K3-35 3 × 106 3 × 105 0.3 SS … … …
K3-36 3 × 106 3 × 106 0.1 SS … 6.44 × 10−3 0.617
K3-37 3 × 106 3 × 106 0.2 SS … 6.84 × 10−3 1.85
K3-38 3 × 106 3 × 106 0.2 SPH … 7.25 × 10−3 1.97
K3-39 3 × 106 3 × 106 0.3 SS … … …
K3-40 107 3 × 105 0.1 SS … 7.75 × 10−3 1.21
K3-41 107 3 × 105 0.1 SS … 7.53 × 10−3 1.18
K3-42 107 3 × 105 0.2 SS … 8.79 × 10−3 3.89
K3-43 107 3 × 105 0.25 SS … … …
K3-44 107 3 × 105 0.3 SS … … …
K3-45 107 3 × 105 0.4 SS … … …
K3-46 107 106 0.1 SPH … 6.41 × 10−3 0.998
K3-47 107 106 0.2 SS … 6.27 × 10−3 2.76
K3-48 107 106 0.25 SS … … …
K3-49 3 × 107 3 × 105 0.1 SS … 5.24 × 10−3 1.3
K3-50 108 3 × 105 0.03 SS … 3.49 × 10−4 0.0237
K3-51 108 3 × 105 0.04 SPH … 2.99 × 10−3 0.312
K3-52 108 3 × 105 0.05 SS … 7.36 × 10−4 0.107
K3-53 108 3 × 105 0.05 SS No relaxation 8.90 × 10−4 0.13
K3-54 108 3 × 105 0.1 SS … 1.94 × 10−3 0.801
K3-55 108 3 × 105 0.2 SS … 3.41 × 10−3 3.98
K3-56 108 3 × 105 0.2 SPH … … …
K3-57 108 3 × 105 0.25 SS … … …
K3-58 108 3 × 105 0.3 SS … … …
K3-59 108 3 × 105 0.4 SS … … …
K3-60 108 3 × 105 0.5 SS … … …
K3-61 108 106 0.04 SS … 5.05 × 10−4 0.0524
K3-62 108 106 0.1 SPH … 3.19 × 10−3 1.31
K3-63 108 106 0.1 SPH … 3.20 × 10−3 1.31
K3-64 108 106 0.1 SS … 1.90 × 10−3 0.78
K3-65 108 106 0.2 SPH … … …
K8-1 3 × 105 3 × 105 0.9 SS … 3.97 × 10−4 0.418
K8-2 3 × 105 3 × 105 SS … 4.63 × 10−4 1.62
K8-3 3 × 105 3 × 105 SS … … …
K8-4 3 × 105 3 × 105 SS … … …
K8-5 3 × 105 3 × 106 1.2 SPH … 3.89 × 10−4 0.618
K8-6 3 × 105 3 × 106 4.7 SPH … … …
K8-7 106 106 0.1 SS … 4.00 × 10−4 0.0246
K8-8 106 106 SS … 4.32 × 10−4 0.841
K8-9 106 106 SS … 4.82 × 10−4 2.65
K8-10 106 106 SS … … …
K8-11 1.4 × 106 1.4 × 106 0.96 SS Fixed central object 4.02 × 10−4 0.843
K8-12 1.4 × 106 1.4 × 106 0.96 SS Fixed central object 3.91 × 10−4 0.819
K8-13 1.4 × 106 1.4 × 106 0.96 SS Fixed central object, tidal truncation 3.75 × 10−4 0.786
K8-14 1.4 × 106 1.4 × 106 0.96 SS … 4.66 × 10−4 0.976
K8-15 1.4 × 106 1.4 × 106 0.96 SS … 3.98 × 10−4 0.834
K8-16 1.4 × 106 1.4 × 106 0.96 SS Fixed central object, stellar evolution 4.19 × 10−4 0.878
K8-17 1.4 × 106 1.4 × 106 0.96 SS Fixed central object, stellar evolution 3.91 × 10−4 0.819
K8-18 1.4 × 106 1.4 × 106 0.96 SS … 4.32 × 10−4 0.904
K8-19 3 × 106 3 × 106 0.2 SPH … 3.31 × 10−4 0.0892
K8-20 3 × 106 3 × 106 SS … 3.43 × 10−4 1.03
K8-21 3 × 106 3 × 106 SPH … 3.43 × 10−4 1.03
K8-22 3 × 106 3 × 106 SS … 3.70 × 10−4 1.12
K8-23 3 × 106 3 × 106 SS … 3.74 × 10−4 1.13
K8-24 3 × 106 3 × 106 SS Fixed central object 2.99 × 10−4 0.9
K8-25 3 × 106 3 × 106 1.2 SPH … 4.22 × 10−4 1.65
K8-26 3 × 106 3 × 106 1.4 SS … 4.27 × 10−4 2.13
K8-27 3 × 106 3 × 106 SS … … …
K8-28 3 × 106 3 × 106 4.7 SPH … … …
K8-29 9 × 106 9 × 106 SS … 3.78 × 10−4 1.78
K8-30 107 3 × 106 SS … 3.61 × 10−4 1.78
K8-31 107 3 × 106 1.4 SS … 3.12 × 10−4 2.55
K8-32 107 3 × 106 1.4 SS … 3.86 × 10−4 3.15
K8-33 107 3 × 106 SS … … …
K8-34 3 × 107 3 × 106 SS … 3.14 × 10−4 2.45
K8-35 3 × 107 3 × 106 1.4 SS … … …
K8-36 108 3 × 106 0.1 SS … 7.49 × 10−5 0.0308
K8-37 108 3 × 106 0.6 SS … 2.40 × 10−4 1.45
K8-38 108 3 × 106 0.8 SS … 2.74 × 10−4 2.55
K8-39 108 3 × 106 SS … … …
K8-40 108 3 × 106 1.4 SS … … …
K8-41 108 3 × 106 SS … … …
MGG9-K8 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 … …
MGG9-K9 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 1.55 × 10−4 1.95
MGG9-K12a 12 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 5.10 × 10−6 0.0644
MGG9-K12b 12 2.7 × 106 2.7 × 106 2.6 SPH Model for MGG-9 5.98 × 10−6 0.0755

When not otherwise mentioned in the ‘Peculiarities’ column, clusters were modelled with a Salpeter (α= 2.35) IMF extending from 0.2 to 120 M, our standard MR relation and prescription for collisional rejuvenation and MS lifetime. The data listed in the column entitled ‘Coll.’ indicates the type of collision treatment used. ‘SS’ stands for sticky-sphere approximation; ‘SPH’ for the prescriptions for mass loss and merger derived from the SPH simulations of Freitag & Benz (2005);‘X’ is a special case for which a simple merger criterion was used [λmerg=−1.14–0.5(lv− 3.0); see Section 3.2.5 of Paper I] and a flat mass-loss rate of 5 per cent was assumed in all cases. tra is the time at which runaway growth of a VMS started. It is given in FP time-units (TFP) and Myr. No value is given for clusters which didn't experience collisional runaway.

Expressions for half-mass and local relaxation times are given in GFR04 and Paper I. We denote the initial values of the half-mass and central relaxation times by trh (0) and trc (0), respectively.

Table 1 lists the initial conditions for all runs performed in this study. Np denotes the number of particles used in the simulation. We also give in this table the value tra of the time when runaway started for all runs in which a VMS formed.

### 2.2 Overview of the results of the standard set of simulations

We first establish the conditions for runaway formation of a VMS using the simplest and most favourable prescription for the collisions, that is, assuming that every time two stars come closer to come closer to each other than the sum of their radii, they merge without any mass loss. With such a clean collision physics, one expects the runaway to happen provided that:

• Core collapse, driven by mass segregation due to two-body relaxation (as studied in GFR04), occurs within t*≃ 3 Myror the cluster is initially sufficiently collisional.

• The collisional cross section is a steeply increasing function of a star's mass, ScollMη* with η > 1.

What ‘sufficiently collisional’ means is not easy to establish from simple principles. One may think that a sufficient condition is that the collision time at the centre, for a star near the top of the IMF (i.e. 120 M for our standard IMF), tcoll|$M max$ (see equation 3), be shorter than the MS lifetime for such a star. But this does not account for the structural evolution of the cluster due to relaxation and collisions themselves. Another condition is suggested by the coagulation equation (Lee 1993, 2000; Malyshkin & Goodman 2001; see also the simplified mathematical analysis encompassed in equations 6 and 7 of Paper I). In the case of strong gravitational focusing, it reduces to R*Mβ* with β > 0 (R* is the stellar radius), a condition always obeyed by MS stars (but not during the pre-MS stage). If gravitational focusing is strong (i.e. for velocity dispersion σv > 500–1000 km s−1), one needs β > 0.5, which is (although marginally) not satisfied by our standard VMS mass–radius relation (MR relation) (β= 0.47, see Paper I). In any case, the relevance of the second condition is not clear because it stems from a type of analysis which disregards the stellar dynamics, most noticeably the role of mass segregation.

Fig. 1 represents the (Mcl, RNB) plane in which are plotted, for W0= 3 and W0= 8, the conditions tcc|rlx= 0.12 trc (0) =t*= 3 Myr and tcoll |$M max$= 3 Myr. tcc|rlx denotes the relaxational core-collapse time, when collisions are absent. The value tcc|rlx= 0.12 trc(0) is the average value yielded by me(ssy)**2 when collisions are absent (see Paper I and discussion below). We also indicate the initial one-dimensional velocity dispersion at the centre,

2

This is approximately correct for any W0 value. By comparing σ0 with the escape velocity from a stellar surface, 500–1000 km s−1, one can estimate whether gravitational focusing (included in all simulations) and collisional mass loss (not accounted for in the sticky-sphere approximation) may play an important role.

Figure 1

Diagram summarizing the initial conditions and outcomes of the ‘standard’ cluster simulations with collisions treated as perfect mergers or using the SPH prescription. Each point represents the cluster mass and (initial) N-body length unit (RNB) for one simulation. Triangles correspond to W0= 3 King models, round dots to W0= 8. The standard 0.2–120 M Salpeter IMF was used in all cases. The solid diagonal lines (of negative slope) show the condition for the relaxation-driven core-collapse time to be 3 Myr, tcc|rlx= 0.12 trc(0) = 3 Myr for W0= 8 (top) and W0= 3 (bottom). The long-dashed lines (of positive slope) show where the central collision time for a 120-M star (with any other star) is 3 Myr. The dot–dashed curves indicate where the time for occurrence of runaway should be 3 Myr, according to our parametrization (see equation 4 and text). Short-dashed lines indicate approximately the 1D central velocity dispersion. Solid symbols are for simulations which resulted in runaway formation of a VMS; open symbols are for cases in which stellar evolution interrupted core collapse before a VMS could grow. In some cases, several simulations with the same RNB and N* (but different W0, initial realizations of the cluster, random sequences for the MC algorithm, number of particles or collision prescription) were carried out. These cases correspond to symbols connected together with a thin line. Asterisks indicate simulations done with a number of particles smaller than the number of stars (3 × 105 or 106 particles for W0= 3 cases, 3 × 106 particles for W0= 8). In the colour online version of this diagram, simulations using the SPH prescriptions for collision outcome are indicated by red points. Among the three models with W0= 3, N*= 108 and RNB= 0.2 pc, one was computed with the sticky-sphere approximation and experienced runaway; the other two, making use of SPH prescriptions, missed the runaway phase. This is essentially the same plane as Fig. 1 of Paper I except that RNB is used instead of Rh and that the plotting domain is shifted to higher masses and smaller sizes to match the conditions that can be treated by the MC approach and may lead to runaway evolution.

Figure 1

Diagram summarizing the initial conditions and outcomes of the ‘standard’ cluster simulations with collisions treated as perfect mergers or using the SPH prescription. Each point represents the cluster mass and (initial) N-body length unit (RNB) for one simulation. Triangles correspond to W0= 3 King models, round dots to W0= 8. The standard 0.2–120 M Salpeter IMF was used in all cases. The solid diagonal lines (of negative slope) show the condition for the relaxation-driven core-collapse time to be 3 Myr, tcc|rlx= 0.12 trc(0) = 3 Myr for W0= 8 (top) and W0= 3 (bottom). The long-dashed lines (of positive slope) show where the central collision time for a 120-M star (with any other star) is 3 Myr. The dot–dashed curves indicate where the time for occurrence of runaway should be 3 Myr, according to our parametrization (see equation 4 and text). Short-dashed lines indicate approximately the 1D central velocity dispersion. Solid symbols are for simulations which resulted in runaway formation of a VMS; open symbols are for cases in which stellar evolution interrupted core collapse before a VMS could grow. In some cases, several simulations with the same RNB and N* (but different W0, initial realizations of the cluster, random sequences for the MC algorithm, number of particles or collision prescription) were carried out. These cases correspond to symbols connected together with a thin line. Asterisks indicate simulations done with a number of particles smaller than the number of stars (3 × 105 or 106 particles for W0= 3 cases, 3 × 106 particles for W0= 8). In the colour online version of this diagram, simulations using the SPH prescriptions for collision outcome are indicated by red points. Among the three models with W0= 3, N*= 108 and RNB= 0.2 pc, one was computed with the sticky-sphere approximation and experienced runaway; the other two, making use of SPH prescriptions, missed the runaway phase. This is essentially the same plane as Fig. 1 of Paper I except that RNB is used instead of Rh and that the plotting domain is shifted to higher masses and smaller sizes to match the conditions that can be treated by the MC approach and may lead to runaway evolution.

From this diagram, it is clear that only relatively small clusters may have tcc|rlx < 3 Myr without being initially collisional (tcoll|$M max$ > 3 Myr); for this one must have N* < 106 for W0= 3 and N* < 3 × 106 for W0= 8. This fact was not included in GFR04, who considered only the tcc|rlx condition. Consequently, there are in principle more clusters which may form VMSs through collisions. Even though initial conditions with collision time shorter than a few million years are questionable, they are an interesting idealization to pave the way to more realistic simulations of the collisional formation of dense clusters.

We show in Fig. 1 the initial conditions for all simulations carried out with the standard IMF and collisions treated as pure mergers or with the smoothed particle hydrodynamics (SPH)-generated prescriptions. We also indicate the outcome of each simulation, that is, whether a VMS grew or whether the core collapse was terminated by stellar evolution of the massive MS stars before this could happen. We set the lower limit for successful VMS growth at MVMS≥ 400 M. This is a relatively arbitrary value which matters only in the case of the four simulations with W0= 3, N*= 105 and RNB= 0.5 pc. These runs were computed with the same code and parameters but using different realizations of the cluster and random sequences. In three runs, a star grew to 400–500 M before it left the MS; in the fourth, only ∼300 M was reached. To determine the time at which the runway started, tra, we look for an increase of mass of a factor of 3 or more in a star more massive than 0.9 Mmax (the maximum mass in the IMF) within the last tenth of the elapsed time. This definition may seem contrived but is required to capture the onset of runaway as it would be naturally identified by inspection of collision history diagrams such as Fig. 12 or 13. In practice tra is very close to the core-collapse time one would identify by looking at the evolution of the Lagrange radii.

Figure 12

Collisional history for run K3-33. We represent the evolution of the five particles that have experienced the largest number of collisions. The top panel represents the radius (distance from the cluster centre) where the collision occurred. The bottom panel shows the evolution of the mass. A circled dot indicates that the star has merged with a larger one (usually the runaway star). Note that there are few collisions until t≃ 2.5 Myr, that is, the moment of core collapse, at which time one star starts growing very quickly.

Figure 12

Collisional history for run K3-33. We represent the evolution of the five particles that have experienced the largest number of collisions. The top panel represents the radius (distance from the cluster centre) where the collision occurred. The bottom panel shows the evolution of the mass. A circled dot indicates that the star has merged with a larger one (usually the runaway star). Note that there are few collisions until t≃ 2.5 Myr, that is, the moment of core collapse, at which time one star starts growing very quickly.

We followed the merger sequence to 1000 M at least in most cases but, except for a few exceptions, did not try to carry on the simulation until the growth was terminated by evolution off the MS. With no or little collisional mass loss, each merger brings new hydrogen to the VMS and allows it to survive until the next merger, see Figs 16 and 17. It is likely that, in real cluster, the growth will saturate through some process not included in me(ssy)**2. One possibility is the depletion of the ‘loss cone’, that is, the disappearance of the stars populating the orbits that intersect the small central volume in which the VMS move around, see Section 2.5.3. Another limiting mechanism could be that the VMS cannot radiate the energy dumped into it by collisions fast enough to keep its relatively small MS size, and instead swells and becomes ‘transparent’ to impactors, as suggested, for other reasons, by Colgate (1967). On the MS, a VMS is not only dense enough to stop impactors but also resilient to disruption despite its radiation-dominated interior. Indeed, the binding energy per unit mass actually increases from ∼3 to ∼4.5 GM/R for masses ranging from 100 to 104 M (Bond, Arnett & Carr 1984). In any case, we are here mostly concerned with detecting the onset of the runaway growth and leave the difficult prediction of its final mass for future studies.

Figure 16

Evolution of mass and various time-scales during the runaway for the same simulations as in Fig. 14. We plot the time between successive collisions, an estimate of the K–H time-scale tKH of the collision product (assuming normal MS MR relation and luminosity), the time left until exhaustion of hydrogen at the centre (from our simple ‘minimal rejuvenation’ prescription), as well as the mass of the star (right scale). No use is made of tKH during the simulations. We plot it here to show than in the vast majority of cases, the interval between collisions is (much) shorter than tKH so that the star would probably be out of thermal equilibrium, with a swollen structure.

Figure 16

Evolution of mass and various time-scales during the runaway for the same simulations as in Fig. 14. We plot the time between successive collisions, an estimate of the K–H time-scale tKH of the collision product (assuming normal MS MR relation and luminosity), the time left until exhaustion of hydrogen at the centre (from our simple ‘minimal rejuvenation’ prescription), as well as the mass of the star (right scale). No use is made of tKH during the simulations. We plot it here to show than in the vast majority of cases, the interval between collisions is (much) shorter than tKH so that the star would probably be out of thermal equilibrium, with a swollen structure.

Figure 17

Mass and time-scales evolution during the growth of the runaway star for the same simulations as in Fig. 15. See caption of Fig. 16 for explanations.

Figure 17

Mass and time-scales evolution during the growth of the runaway star for the same simulations as in Fig. 15. See caption of Fig. 16 for explanations.

Our work is limited to clusters with a relatively large number of particles: N*≥ 105, 3 × 105, for W0= 3, 8, respectively, because the MC code can only yield robust results when there are always at least a few dozens particles in the cluster core; otherwise, one cannot make a robust estimate of the central density, which is required to simulate relaxation and compute the collision probabilities.

An examination of Fig. 1 shows that, in general, the simple expectation derived from the results of GFR04 that VMS formation will occur when and only when the central relaxation time is shorter than 20 Myr is borne out by the present simulations. However, for relatively small clusters (N* < 3 × 105 for W0= 3; N* < 106 for W0= 8), a shorter central relaxation time appears to be required. This is possibly a consequence of the dearth of massive stars in the central regions: with our standard IMF, the number of stars with M* > 50 M in the core of a N*= 3 × 105, W0= 3 or a N*= 106, W0= 8 cluster is ∼30 or ∼20, respectively. Although, over the time required for core collapse ∼100 M stars may come to the centre from distances much larger than the core radius, the dynamics of core collapse is driven by a small number of particles and the applicability of the standard treatment of relaxation (and of the notion of two-body relaxation itself) becomes questionable.

For large clusters, in contrast, we see that runaway VMS formation may happen at sizes for which trc(0) > 20 Myr because collisions occur from the beginning and accelerate core collapse. Fig. 2 is a graphical attempt at quantifying this effect. On this diagram, we plot, for all simulations with our standard IMF which lead to VMS formation, tra, the time of runaway occurrence. Using equation [8–123] of Binney & Tremaine (1987), we define the initial central collision time for a star of 120 M as

3

where

Figure 2

Time of collisional runaway occurrence as a function of central initial collision time for 120-M stars (equation 3), for all simulations with standard Salpeter IMF where runaway happened. Times are given in units of the initial central relaxation time. Triangles and circles are for W0= 3, 8 models respectively. The squares are W0= 9, 12 models for the cluster MGG-9 in M 82 (see Section 2.5.2). These simulations were carried out with a Kroupa IMF extending from 0.1 to 100 M. Solid symbols correspond to the ‘sticky-sphere’ treatment of collisions, open symbols to the SPH prescriptions for collision outcomes. The dashed line is an ad-hoc parametrization of the sticky-sphere results given in equation (4). Note that for large collision times, one gets tra≃ 0.12 trc(0), in good agreement with GFR04 for the relaxational core-collapse time but with a large dispersion in the result. The star corresponds to a simulation with no two-body relaxation but the same initial conditions as those for the triangular point below it (see Fig. 6).

Figure 2

Time of collisional runaway occurrence as a function of central initial collision time for 120-M stars (equation 3), for all simulations with standard Salpeter IMF where runaway happened. Times are given in units of the initial central relaxation time. Triangles and circles are for W0= 3, 8 models respectively. The squares are W0= 9, 12 models for the cluster MGG-9 in M 82 (see Section 2.5.2). These simulations were carried out with a Kroupa IMF extending from 0.1 to 100 M. Solid symbols correspond to the ‘sticky-sphere’ treatment of collisions, open symbols to the SPH prescriptions for collision outcomes. The dashed line is an ad-hoc parametrization of the sticky-sphere results given in equation (4). Note that for large collision times, one gets tra≃ 0.12 trc(0), in good agreement with GFR04 for the relaxational core-collapse time but with a large dispersion in the result. The star corresponds to a simulation with no two-body relaxation but the same initial conditions as those for the triangular point below it (see Fig. 6).

We take M1= 120 M, M2=〈M*〉≃ 0.69 M and R1,2=R*(M1,2) according to our MS MR relation. When collisions are treated as perfect mergers (‘sticky-sphere’ approximation), a reasonable fit to our results for the runaway time is

4

Hence tra is approximately 0.12 trc(0), close to the core-collapse time found in GFR04, 0.15 trc(0), for all clusters with a relatively long initial central collision time. As tcoll,c(0) decreases, so does tra, first mildly and then steeply when tcoll,c(0) < 10−3trc(0). For very short tcoll,c(0), one has tra≃ 250 tcoll,c(0). In this regime, the core collapse appears to be driven by collisions rather than relaxation. This is made clear by comparing the results of two simulations of a cluster with W0= 3, N*= 108 (Np= 3 × 105) and RNB= 0.05 pc. The first run (K3-52) was realized with the usual physics, including two-body relaxation and collisions but, for the second, relaxation was switched off. The evolution of the Lagrange radii is plotted for both the runs in Fig. 6. With relaxation, core collapse occurs at tcc≃ 7.4 × 10−4TFP. When only collisions drive the evolution (run K3-53), it is slightly slower, with tcc≃ 8.9 × 10−4TFP, but is otherwise similar.

Figure 6

Evolution of Lagrange radii for model K3-52. This cluster is strongly collisional from the beginning. The core collapse is driven mostly by mergers and occurs on a time-scale much shorter (when measured in FP units). The dashed lines show the evolution of the same system without relaxation (K3-53). Core collapse still occurs, although on a slightly longer time-scale.

Figure 6

Evolution of Lagrange radii for model K3-52. This cluster is strongly collisional from the beginning. The core collapse is driven mostly by mergers and occurs on a time-scale much shorter (when measured in FP units). The dashed lines show the evolution of the same system without relaxation (K3-53). Core collapse still occurs, although on a slightly longer time-scale.

It is somewhat surprising that the transition to collision-dominated core collapse occurs at such a short value of tcoll,c(0) and that tra/tcoll,c(0) is so long in this regime. To some extent, this is due to our definition of tcoll,c(0), based on the most massive stars in the IMF (M1= 120 M), a choice which may not accurately capture the nature of the cluster evolution in the collisional regime. Indeed stars much less massive also experience collisions from the beginning. Using M1= 1 M, R1= 1 M to define tcoll,c(0), one finds tra≃ 2 tcoll,c(0) when tcoll,c(0) ≪trc(0). In any case this relation for tra only holds for sticky-sphere collisions. When the SPH-inspired prescription is used (with mass loss and merger requiring dmin < R1+R2) and for high velocity dispersions, collisions are less effective at driving core collapse and tra is consequently longer for a given, small tcoll,c(0) value, as is apparent on Fig. 2.

Although it is clearly only approximate, it is tempting to use equation (4) to predict which conditions will lead to collisional runaway by assuming this happens if, and only if tra < 3 Myr. Therefore we have also plotted in Fig. 1 lines indicating this ‘corrected’ runaway condition. One sees that, when the velocity dispersion is sufficiently smaller than V*, this condition reduces to the condition on the central relaxation time alone1 while at high velocities, one predicts runaway for larger cluster sizes than suggested by trc. However, the simulation results show the runaway domain for very massive clusters to reach still slightly larger RNB values. Some of this may be due to dispersion in the tra results (see below) but, by inspecting the run K3-55, another effect is discovered. In this case, the VMS star survives to at least t= 4 Myr, despite its MS lifetime being shorter than 3 Myr. This is due to the combination of (1) the prescription used to set the effective age of merger products by following their central He content (‘minimal rejuvenation’, see Paper I), (2) the relation we use for the value of the core He at the end of the MS as a function of the stellar mass (Fig. 3b of Paper I) and, (3) for this particular simulation, the assumption of perfect merger. Because the MS lifetime is nearly independent of the mass for M*≥ 100 M but more massive stars transform a larger fraction of their mass into He on the MS, each merger amounts to an effective rejuvenation if no mass loss is allowed. Clearly collisional mass loss may change the picture in this very special regime but we think the most important shortcoming of the approach is to assume the VMS to be on the MS. As we will see below, in most cases, the average time between collisions is much shorter than the thermal time-scale of a MS VMS so that it structure is likely strongly affected by collisions and our SPH mass-loss prescriptions would not be appropriate anyway.

Figure 3

Evolution of Lagrange radii for model K8-27. For this cluster, the core collapse is slightly too long to allow the runaway process to set in. The cluster undergoes two successive core collapses. The first one driven by the massive MS stars, the second by the stellar BHs (with a mass of 7 M) resulting from their evolution off the MS. For clarity, the curves have been smoothed by using a sliding average procedure with a (truncated) Gaussian kernel over 10 adjacent data points. See also Fig. 4.

Figure 3

Evolution of Lagrange radii for model K8-27. For this cluster, the core collapse is slightly too long to allow the runaway process to set in. The cluster undergoes two successive core collapses. The first one driven by the massive MS stars, the second by the stellar BHs (with a mass of 7 M) resulting from their evolution off the MS. For clarity, the curves have been smoothed by using a sliding average procedure with a (truncated) Gaussian kernel over 10 adjacent data points. See also Fig. 4.

The extension of the runaway domain to larger systems opens the possibility that proto-galactic nuclei without a massive black hole (MBH) may be subject to runaway formation of a central massive object, despite relatively long relaxation times. Present-day nuclei have a size of a few pc. However, using equation (4) as a guide, a compact RNB= 1 pc nucleus with W0= 8 needs to be at least as massive as ∼3 × 109 M to experience collision-driven collapse in less than 3 Myr. For RNB= 3 pc, the minimum mass is ∼3 × 1010 M! So it seems that this process may operate in galactic nuclei only if they are born very compact or concentrated. Indeed for W0= 10, one would expect all clusters more compact than ∼2.5 pc to complete core collapse within 3 Myr, independent of their mass.

Our results show a relatively large dispersion in the runaway time when evolution is dominated by relaxation, in which case tra should be very close to the core-collapse time found in GFR04, that is, tcc|rlx≃ 0.15 trc(0) (with little variation from one simulation to another). The present results indicate a slightly faster evolution, tcc|rlx≃ 0.12 trc(0), which is very likely the consequence of the use of a different code.2 The larger intrinsic dispersion is also a property of the present code, linked to the use of local time-steps; they impose a scheme in which the selection of the next particle pair to evolve is a random process.

Following this overview of our results for the complete set of standard simulations, we proceed, in the following subsections to a more detailed description of the runaway process based on a few illustrative cases. Presentation of non-standard cases, for which we varied the IMF, MR relation, collision prescription, or the treatment of the VMS, will be done in Section 2.5.

### 2.3 Cluster structure evolution

#### 2.3.1 Missing the collisional runaway

Figs 3 and 4 illustrate what occurs when the core-collapse time is longer than MS lifetime of massive stars. For the sake of clarity, we present a case in which the segregation-induced core collapse would have taken just slightly longer than t*= 3 Myr. The evolution of the Lagrange radii, plotted in Fig. 3, indeed indicates a first core contraction which stopped and reversed at t*. A second collapse occurs much later, at t≃ 500 Myr. What happens is made clear by Fig. 4. On the top panel, we have plotted the evolution of the number of stellar remnants; the bottom panel is the evolution of the Lagrange radii for each stellar type [MS: main-sequence stars, WD: white dwarfs (0.6 M), NS: neutron stars (1.4 M), BH: stellar black holes (7 M)]. The first collapse stops when the massive stars turn into BHs, causing strong mass loss from the central region which re-expands as its binding energy decreases. The BHs are born strongly segregated because their massive progenitors had concentrated at the centre. The BH distribution first re-expands but eventually re-collapses as a result of mass segregation. From t≃ 40 Myr onwards, they are indeed the most massive objects in the cluster. The subsequent evolution of this population of BHs will be driven by their binaries, either present before the second collapse, or dynamically formed through three-body interactions during it. The MC code used here cannot treat binaries so we leave the study of the evolution of the central dense BH subcluster out of the present work (see O'Leary et al. 2006).

Figure 4

Evolution of the various stellar types for model K8-27 (same as Fig. 3). For this cluster, the core collapse is slightly too long to allow the runaway process to set in. The stellar types accounted for are: MS stars (MS), white dwarfs (WD), neutron stars (NS) and stellar BHs (BH). For simplicity, all WDs have 0.6 M, all NS 1.4 M and all BHs 7 M. Top panel: number of compact remnants as function of time. Bottom panel: evolution of Lagrange radii. The simulation was stopped as the stellar BHs underwent core collapse on their own.

Figure 4

Evolution of the various stellar types for model K8-27 (same as Fig. 3). For this cluster, the core collapse is slightly too long to allow the runaway process to set in. The stellar types accounted for are: MS stars (MS), white dwarfs (WD), neutron stars (NS) and stellar BHs (BH). For simplicity, all WDs have 0.6 M, all NS 1.4 M and all BHs 7 M. Top panel: number of compact remnants as function of time. Bottom panel: evolution of Lagrange radii. The simulation was stopped as the stellar BHs underwent core collapse on their own.

#### 2.3.2 Core collapse with runaway

The evolution of Lagrange radii during the core collapse and runaway VMS growth are depicted, for four W0= 3 cases, in Fig. 5. The first simulation is for a cluster which not is collisional initially, that is, from equation (4), with tcoll,c(0)– as defined in equation (3) for 120-M stars – larger than ∼10−3trc(0): tcoll,c(0) = 3.0 × 106 yr = 0.056 trc(0). For this particular cluster, core collapse occurs in approximately 2.5 Myr, thus beating stellar evolution. Collisions are very rare until the runaway starts in deep core collapse, see Fig. 12. Therefore, the evolution to core collapse is very similar to what was obtained for the pure relaxation case in Paper I. Once the runaway has started, the central regions re-expand as a result of the heating effect of the accretion of deeply bound stars by the central object, as in a cluster with a central BH destroying or ejecting stars that come too close (Shapiro 1977; Duncan & Shapiro 1982; Amaro-Seoane, Freitag & Spurzem 2004; Baumgardt, Makino & Ebisuzaki 2004a,b). We cannot follow this ‘post-collapse’ evolution very far because our code slows down considerably when reaching the runaway phase, as the central relaxation and collision times plummet. Furthermore, as will be discussed in Section 2.5.3, it is likely that our way of computing collision probabilities applies correctly to the VMS only at the beginning of its growth.

Figure 5

Evolution of Lagrange radii for various clusters with W0= 3. From the least to the most collisional (initially): with solid lines (black in colour version), K3-33; with long dashes (blue) K3-37; with short dashes (magenta) K3-64; with dotted lines (red) K3-61. The top horizontal and right vertical axis indicate time and radius in physical units for model K3-33. The evolution of this cluster is driven by two-body relaxation until late into core collapse, at which point the collisions kick in. The core-collapse time is thus very similar to the value one would obtain for point-mass dynamics.

Figure 5

Evolution of Lagrange radii for various clusters with W0= 3. From the least to the most collisional (initially): with solid lines (black in colour version), K3-33; with long dashes (blue) K3-37; with short dashes (magenta) K3-64; with dotted lines (red) K3-61. The top horizontal and right vertical axis indicate time and radius in physical units for model K3-33. The evolution of this cluster is driven by two-body relaxation until late into core collapse, at which point the collisions kick in. The core-collapse time is thus very similar to the value one would obtain for point-mass dynamics.

The other simulations plotted in Fig. 5 correspond to denser clusters in which collisions play an increasingly important role in the dynamics and lead to shorter and shorter collapse time (when measured in relaxation times). The central collision times tcoll,c(0) for a 120 M star are 6.2 × 105 yr = 0.014 trc(0), 8.8 × 103 yr = 1.3 × 10−4trc(0) and 500 yr = 2.9 × 10−5trc(0). As explained in Section 2.2, in extremely dense clusters, core collapse is driven by collisions themselves, with relaxation playing only a minor role. This is illustrated in Fig. 6 for another W0= 3 cluster. For these runs, we have assumed that all collisions result in complete mergers. In contrast to gravitational two-body encounters, such mergers do not redistribute energy between stars; they continuously dissipate orbital energy, leading to shrinkage of all inner Lagrange radii with no rebound.

To get a more precise description of the structure of a cluster while it undergoes core collapse and runaway VMS formation, we look at successive radial profiles of the stellar density and velocity dispersion. To have a sufficient resolution of the central regions, which are the most affected by the dynamical evolution, we use our highest resolution simulation, model K8-29, with 9 × 106 particles (the largest number that would fit into the memory of available computers). The initial structure is a W0= 8 King model. Fig. 7 shows the evolution of the total stellar mass density. The core collapse is characterized by the formation a small central density peak, growing inside the initial constant-density core. In comparison with the core collapse of single-mass models, the evolution shown here is relatively unremarkable and does not appear to proceed in a self-similar manner. In particular, the contracting core does not leave behind a ρ∝R−γ region with γ≃ 2.23 (Freitag & Benz 2001; Baumgardt et al. 2003, and references therein). The central density profile in deep core collapse is better fitted with γ≃ 1.75. It may be tempting to interpret this value of the exponent by the formation of a Bahcall–Wolf cusp (Bahcall & Wolf 1976), similarly to what two-body relaxation creates in a stellar system dominated by a central massive object. This would not be correct, however, as the following figures make clear. When we look at the evolution of the distribution of stars with masses in the range between 3 and 120 M in Fig. 8, we see a much stronger evolution than for the overall distribution because lighter stars, which form an essentially static background, are so much more numerous. Fig. 9 yields a detailed picture of the distribution of stars in various mass ranges during the late stages of the core collapse (while the VMS is quickly growing). Clearly, only the spatial distribution of stars more massive than 10 M significantly changes and the central density peak is mostly made of stars more massive than 30 M. These objects do not form a γ≃ 1.75 cusp but rather a distribution of varying slope.

Figure 7

Evolution of the density profile for all stars with masses between 0.2 and 120 M (i.e. excluding the runaway object) in simulation K8-29. For each profile in the sequence, we indicate the time and mass of the most massive star in the cluster at this time. The FP time-unit is TFP≃ 4.72 × 109 yr. A central subcluster with a cuspy density profile forms. Its profile is compatible with a power-law of exponent −1.75 but it is not a Bahcall & Wolf (1976) cusp (see Figs 8–10 and text).

Figure 7

Evolution of the density profile for all stars with masses between 0.2 and 120 M (i.e. excluding the runaway object) in simulation K8-29. For each profile in the sequence, we indicate the time and mass of the most massive star in the cluster at this time. The FP time-unit is TFP≃ 4.72 × 109 yr. A central subcluster with a cuspy density profile forms. Its profile is compatible with a power-law of exponent −1.75 but it is not a Bahcall & Wolf (1976) cusp (see Figs 8–10 and text).

Figure 8

Similar to Fig. 7 but we only plot the density of stars with masses between 3 and 120 M. This figure confirms that the most massive stars are responsible for the central detached cusp seen in Fig. 7. However, it seems that the innermost part of the density profile is more concentrated than a ρ∝R−1.75 cusp (see Fig. 9).

Figure 8

Similar to Fig. 7 but we only plot the density of stars with masses between 3 and 120 M. This figure confirms that the most massive stars are responsible for the central detached cusp seen in Fig. 7. However, it seems that the innermost part of the density profile is more concentrated than a ρ∝R−1.75 cusp (see Fig. 9).

Figure 9

Density profiles for various stellar masses in the same simulation as in Figs 7 and 8. To investigate the structure of the cluster during the runaway phase, we have averaged the density profiles for 18 snapshots within the given time interval. The various curves correspond to different ranges in stellar mass. It is clear that only stars more massive than 10 M (with a number fraction of only ∼5 × 10−3) contribute to the formation of the cusp and that its profile is steeper than ρ∝R−1.75. The line with slope -2 is plotted as a guide. Note that the radius and density scales and ranges are different from those of Figs 7 and 8.

Figure 9

Density profiles for various stellar masses in the same simulation as in Figs 7 and 8. To investigate the structure of the cluster during the runaway phase, we have averaged the density profiles for 18 snapshots within the given time interval. The various curves correspond to different ranges in stellar mass. It is clear that only stars more massive than 10 M (with a number fraction of only ∼5 × 10−3) contribute to the formation of the cusp and that its profile is steeper than ρ∝R−1.75. The line with slope -2 is plotted as a guide. Note that the radius and density scales and ranges are different from those of Figs 7 and 8.

In Fig. 10 are plotted the velocity dispersion profiles for the same ranges of stellar masses as in Fig. 9. Most importantly, this figure shows that, inside a region roughly corresponding to the initial core [Rcore(0) = 0.12 pc], relaxation has produced partial energy equipartition with stars above 10 M having cooled to lower velocities by heating up lighter stars. This is particularly dramatic for the stars in the 30–120 M bin, which exhibit an inverse temperature gradient. Around R= 0.2 the velocity dispersion for these massive stars has in fact increased above the initial value, because the stars with higher orbital velocity are less affected by dynamical friction and have stayed in place while others were sinking. In this figure, we also plot the Keplerian velocity for a MVMS= 1000 M central object, . Clearly, we cannot resolve the radius of influence (the central region dominated by the central object's gravity) for any range of stellar mass. Only within the radius of influence would one expect a Bahcall–Wolf cusp to develop. The evolution of the central density and velocity dispersion (for another simulation) is depicted in Fig. 11. We see how the central density steadily increases while (in contrast with the single-mass core collapse) the velocity dispersion drops, a result of energy equipartition. Both effects contribute to making collisions more and more likely. Furthermore, the decrease in velocity dispersion ensures that disruptive collisions will be unlikely unless the central velocity dispersion is initially very high.

Figure 10

Similar to Fig. 9 but for the 1D velocity dispersion σv (plotted on a linear scale). As a result of (partial) energy equipartition, the most massive stars are considerably cooler than other stars in the central regions. The dotted line is the initial velocity dispersion (W0= 8 King model). The dash–dotted line indicates the Keplerian dispersion corresponding to a central object with MVMS= 1000 M. One sees that the sphere of influence of the runaway star (inside which σv≃σK) is not resolved.

Figure 10

Similar to Fig. 9 but for the 1D velocity dispersion σv (plotted on a linear scale). As a result of (partial) energy equipartition, the most massive stars are considerably cooler than other stars in the central regions. The dotted line is the initial velocity dispersion (W0= 8 King model). The dash–dotted line indicates the Keplerian dispersion corresponding to a central object with MVMS= 1000 M. One sees that the sphere of influence of the runaway star (inside which σv≃σK) is not resolved.

Figure 11

Evolution of the central density (ρc, dashed line, left scale) and 3D velocity dispersion (σ3c, solid line, right scale) for model K3-37. Note how the central parts cool down during core collapse as a result of partial equipartition between the massive stars and the lighter ones.

Figure 11

Evolution of the central density (ρc, dashed line, left scale) and 3D velocity dispersion (σ3c, solid line, right scale) for model K3-37. Note how the central parts cool down during core collapse as a result of partial equipartition between the massive stars and the lighter ones.

### 2.4 Collisional runaway sequences

We first describe results obtained in the sticky-sphere approximation for collisions. We will then present the differences introduced by the more realistic collision treatment based on SPH results.

A typical runaway history for a cluster which is not initially collisional is depicted in Fig. 12, where we follow five stars participating in collisions, including the runaway merger product. The parameters of the cluster are such that it undergoes core collapse in ∼2.5 Myr. Only a few collisions happen before deep core collapse, at which time one star detaches itself from the rest of the population by growing in a rapid succession of mergers. Two characteristics of this runaway growth are its steepness and the fact that no other star experiences significant mass increase. The initial central velocity dispersion in this cluster is σ0≃ 52 km s−1, a value for which we expect the sticky-sphere approximation to be accurate. Fig. 13 is an example of runaway in a cluster which is initially strongly collisional. In this case, many stars experience numerous collisions and start growing from the beginning of the evolution. Nevertheless, only one of these objects experiences a very fast growth at the moment of (collision-driven) core collapse. Here σ0≃ 1200 km s−1 but the escape velocity from a 1000 M VMS on the MS is Vesc≃ 3000 km s−1, so gravitational focusing is still effective, producing the steep increase of cross section with MVMS required for runaway. Because relative velocities are still significantly lower than Vesc, we do not expect strong collisional mass loss but the critical dmin for merger becomes significantly smaller than R1+R2, making mergers less frequent than assumed here.

Figure 13

Collisional history for run K3-52. Collisions occur from the beginning of the evolution and many stars grow in an orderly way until one detaches itself from the others and starts growing in a runaway fashion. See text for explanations.

Figure 13

Collisional history for run K3-52. Collisions occur from the beginning of the evolution and many stars grow in an orderly way until one detaches itself from the others and starts growing in a runaway fashion. See text for explanations.

To study which stars contribute to the build-up of the VMS, we draw merger trees. Four examples are presented in Figs 14 and 15. The first is for two clusters with W0= 3 and N*=Np= 3 × 105, the second for W0= 3 and N*=Np= 3 × 106. In each figure, panel (b) shows the case of a smaller, initially more collisional cluster. In all cases, we follow the VMS growth up to 2000 M. The cluster in Fig. 14(a) has RNB= 0.4 pc, leading to collapse in tcc= 2.9 Myr, just in time to produce collisional runaway before the massive stars would have left the MS. In this situation, typical of clusters which first become collisional in deep collapse, the merger tree has a very simple structure, with most branches lacking substructure: the stars that merge with the VMS (the trunk) are not themselves collision products. Furthermore, a large majority of these impactors have masses around the tip of the IMF, between 60 and 120 M. These stars add up to 92 per cent of the VMS mass. The distribution of the number of contributing stars (the tips of the branches, that is, stars that are not themselves collision products) is bimodal with one peak between 0.2 and 1 M and another covering the 40–120 M; it reflects the MF in the central regions (GFR04). This is to be contrasted with Fig. 14(b), obtained for a much more collisional cluster with RNB= 0.03 pc, which enters the runaway phase at t≃ 0.064 Myr. Here most stars merging with the VMS were themselves produced by previous mergers and the contribution of lighter stars is larger, although the 60–120 M range still contributes 78 per cent of the VMS mass. The cause of the difference is that, in the more compact cluster, collisions occur from the beginning, when the cluster is not segregated yet; in the larger cluster, virtually all collisions take place in the high-density core formed by the concentration of high-mass stars. The two clusters in Fig. 15, with 3 × 106 stars and sizes of RNB= 0.2 pc and 0.1 pc are both relatively collisional, as required for this W0 value and this number of stars to experience core collapse in less than 3 Myr (see Fig. 1); thus, they exhibit complex merger trees. Again, in the more compact cluster, low-mass stars play a larger but not predominant role, with some 90 per cent of the VMS mass originating in stars more massive than 60 M. Merger trees for W0= 8 clusters exhibit the same general characteristics and variety. In particular, the same kind of bimodal distribution in the masses of the contributing stars can be observed, with more stars in the 0.2–1 M range for clusters with shorter initial tcoll/trlx.

Figure 14

Merger trees for simulations of clusters with W0= 3, N*=Np= 3 × 105. Panel (a): RNB= 0.4 pc (model K3-23). Panel (b): RNB= 0.03 pc (model K3-9). We follow the growth of the runaway star to ∼2000 M. The horizontal axis indicates the mass of stars taking part in the collisions. The sequence of mergers is represented in the vertical direction from top to bottom. The right ‘trunk’ is the growing VMS. The radius of the discs is proportional to that of the corresponding star. In case (a), the cluster is not initially collisional (the collision time for a 120-M star is much larger than its MS lifetime) and most collisions occurs in deep collapse and feature stars of mass 60–120 M that have segregated to the centre. Most of them are not themselves merger products. In case (b), the cluster is initially collisional and, although the runaway also occurs in core collapse, most stars contributing to it have experienced earlier collisions.

Figure 14

Merger trees for simulations of clusters with W0= 3, N*=Np= 3 × 105. Panel (a): RNB= 0.4 pc (model K3-23). Panel (b): RNB= 0.03 pc (model K3-9). We follow the growth of the runaway star to ∼2000 M. The horizontal axis indicates the mass of stars taking part in the collisions. The sequence of mergers is represented in the vertical direction from top to bottom. The right ‘trunk’ is the growing VMS. The radius of the discs is proportional to that of the corresponding star. In case (a), the cluster is not initially collisional (the collision time for a 120-M star is much larger than its MS lifetime) and most collisions occurs in deep collapse and feature stars of mass 60–120 M that have segregated to the centre. Most of them are not themselves merger products. In case (b), the cluster is initially collisional and, although the runaway also occurs in core collapse, most stars contributing to it have experienced earlier collisions.

Figure 15

Merger trees for the simulation of a cluster with W0= 3, N*=Np= 3 × 106. Panel (a): RNB= 0.2 pc (model K3-37). Panel (b): RNB= 0.1 pc (model K3-36). We follow the growth of the runaway star to ∼2000 M. Both clusters have initial central collision times for 120-M stars shorter than 3 Myr, hence the relatively complex merger trees. The more compact cluster, case (b), is more collisional which explains a stronger contribution of low-mass stars.

Figure 15

Merger trees for the simulation of a cluster with W0= 3, N*=Np= 3 × 106. Panel (a): RNB= 0.2 pc (model K3-37). Panel (b): RNB= 0.1 pc (model K3-36). We follow the growth of the runaway star to ∼2000 M. Both clusters have initial central collision times for 120-M stars shorter than 3 Myr, hence the relatively complex merger trees. The more compact cluster, case (b), is more collisional which explains a stronger contribution of low-mass stars.

It is instructive to examine the various time-scales involved in the runaway process. This is done in Figs 16 and 17 for the same four simulations as in Figs 14 and 15. In these diagrams, we follow the VMS, merger after merger, making use of two ordinate axis, one to indicate its mass, the other to monitor important time-scales. Those are the time to the next collision Δtcoll, the time left on the MS ΔtMS (if the star were left to evolve without further collision), and the MS thermal relaxation time, that is, the Kelvin–Helmholtz (K–H) time. The latter is approximated here by tKH=GM2*(R*L*)−1 where L* is the luminosity of the star (close to the Eddington limit for VMS; see, e.g. Baraffe, Heger & Woosley 2001; Schaerer 2002). From 100 to 104 M, tKH increases slowly from ∼104 to 6 × 104 yr. Here we do not truncate at MVMS= 2000 M, but show instead the full extent of the simulations. However, in all but the first case (Fig. 16a), we stopped the simulation before the VMS left the MS, when the average interval between mergers was still much shorter than the remaining MS lifetime.

Because it corresponds to a core-collapse time just slightly shorter than 3 Myr, the first simulation of the four is one of the exceptional few for which we could afford to integrate until stellar evolution caught up with collisions, as indicated by the convergence of the Δtcoll and ΔtMS curves. In most other cases, ΔtMS still exceeded ΔtMS by orders of magnitude – suggesting many more collisions to come – when we decided to terminate the simulation. From the merger trees and the curve showing MVMS as function of the collision number, there is only little sign of a decrease with time, as massive stars that have segregated to the centre are progressively depleted. Such trend can be seen for the runs with 3 × 105 stars (in particular if one draws the complete merger tree past MVMS= 2000 M), but is not perceptible for a cluster of 3 × 106 stars whose central parts constitute a larger reservoir of massive stars. Indeed among 3 × 105 stars forming a W0= 3 cluster, there is only about 1600 M of stars more massive than 60 M in the core.

An essential result is that, in the vast majority of cases the K–H time is (much) longer that the average time between successive collisions. This means that the VMS is unlikely to relax to thermal equilibrium, that is, to the MS after it has merged with an impactor and before the next collision. In principle, then, it may be kept in a swollen state by repeated collisions and our assumptions of MS MR relation and continuous nuclear burning (to set ΔtMS) are probably incorrect during runaway growth. This suggests a very different picture concerning the evolution of the VMS in the runaway phase and the termination of the growth. The VMS may grow larger and larger until it is so diffuse that it lets the smaller stars fly through it without the ability to stop them, hence bringing back the ‘transparency saturation’ (Colgate 1967; Lightman & Shapiro 1978). Another important effect is that the central nuclear reactions are very likely switched off in a merger product which expands as a result of shock heating (e.g. Lombardi et al. 2002, 2003), so the stellar evolution of the VMS may be suspended. Of course these effects will persist only as long as Δtcoll < tKH. When the collisions eventually become rarer – from exhaustion of stars to collide with or because stellar evolution from normal massive stars drives core re-expansion at t≃ 3 Myr– the VMS returns to the MS and to nuclear burning. Given the complexity of interplay between repeated collisions and VMS structure (itself determining the probability and outcome of further collisions), it seems impossible to make any simple prediction as to what the final VMS mass would be when it eventually leaves the MS. In future works, this situation may be studied by combining a stellar dynamics code with on-the-fly SPH simulation to follow the VMS structure collision after collision (Baumgardt & Nakasato, in preparation). For the time being, these aspects have to be considered as one central source of uncertainty regarding the outcome of the runaway phase.

Putting aside the possible complications due to the non-MS structure of the VMS, we can look at the changes in the runaway process introduced by a more realistic treatment of collisions between MS stars, as allowed by the SPH results. Figs 18 and 19 show when and why the sticky sphere approximation may become questionable. These diagrams show the impact parameter and relative velocity for all collisions that happened in two cluster simulations, one with σ0≃ 110 km s−1 (Fig. 18), the other with σ0≃ 635 km s−1 (Fig. 19). For lower velocity dispersion, a large majority of encounters happened with parameters leading to merger according to equation (13) of Paper I. In particular, all collisions but one in which the runaway object took part were in this regime. They had relative velocities smaller than 0.1 V* and would have produced less than 1–3 per cent mass loss (Freitag & Benz 2005). Consequently, treating all collisions as perfect mergers is certainly a very good approximation for system with comparable or lower velocity dispersion, including any reasonable globular cluster. In the high-velocity situation, a significant fraction of the collisions (featuring the runaway object or not) should have left the stars unbound rather than produced a single merged object, suggesting that using the SPH prescription to set the impact parameter for merger may affect the results. On the other hand, most collisions occurred at less than V* and would have produced no more than a few per cent mass loss.

Figure 18

Collision parameters for a cluster with relatively low velocity dispersion (same simulation as in Figs 15a and 17a). For each collision (dot), we indicate the ‘impact parameter’, in units of the stellar half-mass radii, and the relative velocity at large separation, in units of V(h)*; see equation 13 of Paper I). Circled dots indicate mergers involving a star with a mass larger than 120 M (i.e. the runaway object). The lines indicate the condition for merger (equation 13 of Paper I) for [M1, M2]/M=[1, 1] (solid line), [1, 10], (dots), [0.1, 1], (short dashes), [0.1, 60], (long dashes) and [60, 60], (dot–dashed). Only below these lines should the collisions result in mergers.

Figure 18

Collision parameters for a cluster with relatively low velocity dispersion (same simulation as in Figs 15a and 17a). For each collision (dot), we indicate the ‘impact parameter’, in units of the stellar half-mass radii, and the relative velocity at large separation, in units of V(h)*; see equation 13 of Paper I). Circled dots indicate mergers involving a star with a mass larger than 120 M (i.e. the runaway object). The lines indicate the condition for merger (equation 13 of Paper I) for [M1, M2]/M=[1, 1] (solid line), [1, 10], (dots), [0.1, 1], (short dashes), [0.1, 60], (long dashes) and [60, 60], (dot–dashed). Only below these lines should the collisions result in mergers.

Figure 19

Same as Fig. 18 but for a cluster with high velocity dispersion: model K3-55. For this case, a large fraction of collisions would actually be fly-bys and the assumption of perfect mergers is questionable.

Figure 19

Same as Fig. 18 but for a cluster with high velocity dispersion: model K3-55. For this case, a large fraction of collisions would actually be fly-bys and the assumption of perfect mergers is questionable.

To assess the impact of mass loss on the runaway mechanism, we have performed a few runs for which we assumed that all collisions result in a merger with some constant fractional mass loss δM, varying the value of δM. We show runaway growth sequences for δM= 0, 0.05 and 0.10 in Fig. 20. One sees that δM≥ 0.10 is required to prevent the growth of a star more massive than a few 100 M. As typical impactors have a mass around Mimp≈ 100 M, it is obvious why a VMS of 1000 M can only be grown through a sequence of merger if δM < 0.1. The fact that we assume all mass lost to originate from the hydrogen envelopes also shortens the remaining MS lifetime of the collision product and prevents reaching MVMSMimpM.

Figure 20

Role of mass loss in the runaway process. For the three simulations presented here, we used the same initial model with W0= 3, N*= 3 × 105 and RNB= 0.3 pc (models K3-13, K3-14 and K3-15). All collisions are assumed to result in merger with a constant fractional mass loss δM≡δM/M, that is, Mmerger= (1 −δM)(M1+M2) where M1, M2 are the masses of the colliding stars. This prescription is artificial in that a collision with a light star will cause more damage than one with a more massive star. For case (a), with no mass loss, the simulation was interrupted before the runaway object evolved off the MS because mergers constantly rejuvenate the star. With some mass loss, stellar evolution can catch up with collisional growth and terminate it. For simplicity, in these cases, we did not assume an IMBH was formed at the end of the VMS MS, but a 7-M BH instead, which was not allowed to collide with other stars.

Figure 20

Role of mass loss in the runaway process. For the three simulations presented here, we used the same initial model with W0= 3, N*= 3 × 105 and RNB= 0.3 pc (models K3-13, K3-14 and K3-15). All collisions are assumed to result in merger with a constant fractional mass loss δM≡δM/M, that is, Mmerger= (1 −δM)(M1+M2) where M1, M2 are the masses of the colliding stars. This prescription is artificial in that a collision with a light star will cause more damage than one with a more massive star. For case (a), with no mass loss, the simulation was interrupted before the runaway object evolved off the MS because mergers constantly rejuvenate the star. With some mass loss, stellar evolution can catch up with collisional growth and terminate it. For simplicity, in these cases, we did not assume an IMBH was formed at the end of the VMS MS, but a 7-M BH instead, which was not allowed to collide with other stars.

SPH simulations show mass loss of at most a few percent for Vrel < V*, so reaching δM≈ 0.1 (on average) would require really extreme conditions. The cluster with the highest velocity dispersion for which we have performed a simulation using the SPH prescriptions is model K3-51, with σ0≃ 1270 km s−1. It has the same parameters as the sticky-sphere run K3-61 and we obtain the runaway VMS growth in a qualitatively similar fashion, although about four times more collisions are required to attain a given MVMS. The most noticeable difference is that the time required to enter the runaway phase tra is also about four times longer than in the sticky-sphere approximation (see left-most open triangle in Fig. 2). This is obviously a consequence of the collisions being about four times less efficient at driving collapse in this collision-dominated cluster. In cases with lower velocity dispersion, using the SPH prescription make very little change to the conditions required to achieve fast collapse and collisional runaway or to the runaway itself.

In Fig. 21, we plot the cumulative fractional mass loss for all collisions that took place in four simulations for clusters with initial central velocity dispersions ranging from σ0= 28 to 1420 km s−1. We see that all collisions yield δM < 10 per cent for σ0≲ 100 km s−1 with an average 〈δM〉 < 1 per cent. The cluster with the highest σ0 experiences a few completely disruptive collisions but only involving stars less massive than 2 M. Even in this extreme case, 〈δM〉 is only of the order of 5 per cent. For collisions involving a star more massive than 120 M, the average mass loss is as low as 0.3 per cent.

Figure 21

Cumulative distribution of fractional mass loss in collisions for four cluster simulations using the SPH prescription. In all four cases, runaway VMS formation happened. The average fractional mass losses for these runs are 0.2, 0.4, 0.5 and 4.9 per cent.

Figure 21

Cumulative distribution of fractional mass loss in collisions for four cluster simulations using the SPH prescription. In all four cases, runaway VMS formation happened. The average fractional mass losses for these runs are 0.2, 0.4, 0.5 and 4.9 per cent.

The only case for which we find a qualitatively different evolution when using SPH-generated collision recipes instead of pure mergers, is for W0= 3, N*= 108 and RNB= 0.2 pc. In K3-55, assuming pure mergers, we found runway, thanks to collisional rejuvenation that kept the VMS on the MS for more than 4 Myr. With collisional mass loss and a large fraction of collisions not resulting in mergers, the (orderly) growth of stars stopped at t≃ 3.5 Myr, before any runaway star detached from the mass spectrum (runs K3-65 and K3-56).

Because a realistic treatment of collisions has so little impact on most results pertaining to the cluster evolution and runaway VMS growth, it is no surprise that we find them essentially unchanged if, instead of our SPH-based prescription, we implement the simple fitting formulae published by Rauch (1999), which were derived from a smaller set of SPH simulations.

### 2.5 Non-standard simulations

#### 2.5.1 Other IMFs

GFR04 established that, for any realistic IMF and all cluster models we tried for which trc(0) > 0, the segregation-driven core collapse occurs at tcc|rlx≃ 0.15 trc(0). By ‘realistic’ IMF we mean one that is sufficiently broad, with μ≡M*,max/〈M*〉 > 50, and not too steep. This includes the ‘Kroupa’ IMF (Kroupa, Tout & Gilmore 1993), which contains significantly fewer high-mass stars than the Salpeter IMF (with dN*/dM*M−2.7* for M* > 1 M). Although this is actually an effective IMF for stars forming the Galactic field, with clusters probably having a Salpeter-like IMF beyond 1 M (Kroupa & Weidner 2003; Weidner & Kroupa 2004), we used it in order to investigate the possible effects of a smaller fraction of massive stars. We find that the runaway growth of a VMS happens essentially in the same way as with the Salpeter IMF. The only noticeable, and not unexpected, difference is that because stars around 100 M are so much rarer, the average impactor mass decreases much more quickly during the VMS growth, as can be seen in the merger tree shown in Fig. 22.

Figure 22

Merger tree for the runaway growth of a VMS in model K3-61, a cluster with steep Kroupa IMF extending from 0.01 to 120 M (see text).

Figure 22

Merger tree for the runaway growth of a VMS in model K3-61, a cluster with steep Kroupa IMF extending from 0.01 to 120 M (see text).

We have also simulated clusters with a reduced stellar mass range. The single-mass case lacks realism and was only considered in Paper I as a test-case to compare with Quinlan & Shapiro (1990).

If the Salpeter IMF only extends from 0.2 to 10 M, the average stellar mass is 0.58 M and μ= 10/0.58 ≃ 17. The results of GFR04 suggest tcc|rlx≃ 0.7 trc(0) for W0= 3. For such a cluster, we find tra≃ 1.2 trc(0) (run K3-20). The 10-M stars leave the MS after t*= 25 Myr; these two nearly equally longer time-scales combine to give conditions for runaway in the (N*, RNB) plane that are very close to those for our standard IMF. Hence, the details of the upper end of the IMF have little effect on the conditions for the onset of runaway collisions. Of course, if the maximum stellar mass in the IMF is lower, a larger number of mergers are required to grow a VMS but there is about 3 times more mass in stars in the range 5–10 M in a Salpeter IMF truncated at 10 M than in the range 60–120 M if the IMF extends to 120 M. In simulation K3-20, the cluster entered the runaway regime at tratcc|rlx= 20 Myr, allowing a VMS to reach 1500 M when we stopped the computation.

#### 2.5.2 Conditions for collisional runaway in MGG-9

Portegies Zwart et al. (2004) have performed N-body simulations of the clusters MGG-9 and MGG-11 in M 82 using ∼130 000 to ∼600 000 particles. Following McCrady, Gilbert & Graham (2003), they assume for MGG-9 an initial mass of 1.8 × 106 M and a half-mass radius of 2.6 pc. Unlike the case of MGG-11 for which there are indications of a lower cut-off of the IMF at ∼1 M (which yields 〈M*〉≃ 3 M), the observational data for MGG-9 is compatible with a normal cluster Kroupa IMF extending from 0.1 to 100 M (Kroupa 2001): dN*/dM*M−α* with α= 1.3 below 0.5 M and α= 2.3 for higher masses (〈M*〉≃ 0.64 M). To get the total cluster mass as estimated observationally from its size and velocity dispersion, some 2.7 × 106 stars are required. Because this number is too high for present-day direct N-body simulations, Portegies Zwart et al. (2004) adopted a 1 M cut-off. They considered King models with W0 ranging from 6 to 12; none of them led to collisional runaway, in contrast to their MGG-11 models.

In light of GFR04, this is a surprising result. From its position in Fig. 1 of Paper I, the MGG-9 model is expected to undergo quick core collapse and runaway collisions if its initial concentration corresponds to W0≳ 8.5. We have run three simulations of this cluster with Np=N*= 2.7 × 106 for W0= 8, 9 and 12 (MGG9-K8, MGG9-K9 and MGG9-K12a). The model with the lowest concentration does not experience runaway but the other two do. We did not try to simulate MGG-11. If the average stellar mass is as high as 3 M, the number of stars in that cluster is of the order of ∼105, too small for robust MC treatment with a broad IMF and W0≳ 3. We note, however that the condition for runaway suggested by GFR04, tcc|rlx < t*= 3 Myr, predicts that many more of the MGG-11 models considered by Portegies Zwart et al. should have experienced collisional runaway (see Fig. 23). This disagreement indicates again that our simple tcc|rlx criterion only applies to systems containing a sufficiently large number of stars. In particular, for relatively low N*, binaries can form through three-body processes and suspend core collapse before collisions occur. Furthermore, systems of very high concentration and (relatively) low number of stars, contain a very low number of high-mass stars in their core, initially. Hence, a clean mass-segregation-induced core collapse of the type found in our MC simulations may not be possible.

Figure 23

Conditions for collisional runaway for models of the clusters MGG-9 and MGG-11 in M 82 (McCrady et al. 2003). This diagram is adapted from Fig. 3 of Portegies Zwart et al. (2004). Their dynamical friction time tdf is approximately equal to 4.7 × 10−3TFP for models with 〈M*〉= 3 M and 1.3 × 10−3TFP for 〈M*〉= 0.64 M. c is the initial concentration parameter; corresponding W0 values are given on the right of the diagram. Round points are for simulations by Portegies Zwart et al. (〈M*〉= 3 M), squares for our MGG-9 MC runs (〈M*〉= 0.64 M). Solid symbols indicate that the model underwent collisional runaway. The models on the left side (low tdf values) are for MGG-11, the ones on the right for MGG-9. The solid and dash–dotted lines indicate where results of GFR04 predict mass-segregation-driven core collapse to happen at t=t*= 3 Myr for the two different IMFs.

Figure 23

Conditions for collisional runaway for models of the clusters MGG-9 and MGG-11 in M 82 (McCrady et al. 2003). This diagram is adapted from Fig. 3 of Portegies Zwart et al. (2004). Their dynamical friction time tdf is approximately equal to 4.7 × 10−3TFP for models with 〈M*〉= 3 M and 1.3 × 10−3TFP for 〈M*〉= 0.64 M. c is the initial concentration parameter; corresponding W0 values are given on the right of the diagram. Round points are for simulations by Portegies Zwart et al. (〈M*〉= 3 M), squares for our MGG-9 MC runs (〈M*〉= 0.64 M). Solid symbols indicate that the model underwent collisional runaway. The models on the left side (low tdf values) are for MGG-11, the ones on the right for MGG-9. The solid and dash–dotted lines indicate where results of GFR04 predict mass-segregation-driven core collapse to happen at t=t*= 3 Myr for the two different IMFs.

#### 2.5.3 Loss-cone effects

The MC treatment of collisions is based on the assumption that the two neighbouring particles picked up at a given step as potential collision partners are members of the same (local) distribution. In this approach, one does not actually test whether the orbits of the two particles would lead them to collide with each other. Instead, it is based on the average collision probability between such stars in an imaginary volume over which the properties of the stellar population (density, mass spectrum, velocity distributions…) should be homogeneous. For this method to yield accurate collision rates, the cluster properties must be relatively constant over the range of radii used to estimate the local stellar density.3 A more subtle condition, in cases where one particle has unique properties and cannot be thought of as a typical member of some continuous distribution, is that its orbit should allow it to explore a significant fraction of this imaginary volume (for which the stellar density and collisions times are computed) so that its collision probability is also relatively homogeneous.

Unfortunately, these conditions may cease to be fulfilled when the runaway star has grown to such a mass that it basically stays at the centre of the cluster with little radial motion in comparison with the size of the orbits of nearby stars. The situation then becomes similar to that of a massive (or intermediate-mass) BH at a the centre of a galactic nucleus (or globular cluster). While the VMS grows by merging with stars coming into physical contact with it, the (I) MBH destroys stars venturing into the Roche zone and accretes (some of) their material. In both cases, stars must come very close to the cluster centre, that is, have very small pericentre distance Rperi, so the collision or disruption rate is dominated by stars on high-eccentricity orbits. How these orbits are repopulated by two-body relaxation is known as the ‘loss-cone replenishment problem’ (Frank & Rees 1976; Lightman & Shapiro 1977). The basic difficulty is that relaxation may cause large relative changes in Rperi in only a small fraction of the relaxation time, that is, the time over which the energy of the orbit would be strongly affected. To give an illustration, for a Keplerian orbit of eccentricity e≈ 1, The ‘pericentre relaxation time’ is of the order of tr,peri≃ (1 −e)trlx where trlx is the usual (energy) relaxation time. Let Porb be the orbital period of the star (from one pericentre passage to the next). For regions of phase space where Porbtr,peri, the loss cone is empty, while in the opposite regime it is full (since relaxation replenishes loss-cone orbits as quickly as stars are removed through interactions with the central object). According to standard loss-cone theory, in most situations, the disruption rate is dominated by stars at the boundary between the two regimes. Hence, to include loss-cone physics correctly into stellar dynamical models, one should in principle resolve the effects of relaxation on time-scales much shorter than trlx, actually as short as Porb for stars in the critical regime. The MC scheme used here does not allow this directly because the time-step δt of a given particle can only be a function of its present radial position R, not of its eccentricity. Furthermore, δt(R) must be an increasing function of R so that setting δt=Porb at any given R would impose time-steps at least as short as this for all particles interior to R. To handle the loss-cone problem in a satisfying way despite of this, an approximate treatment of the effect of relaxation on time-scales shorter than δt has been devised; it is described in detail in Freitag & Benz (2002). But this approach, as all standard loss-cone theories, is based on the assumption that the central object is rigorously fixed at R= 0; in other terms, no account is made of its ‘wandering’ through Brownian motion.

In summary, our present code allows us to follow the growth of a VMS using two opposite simplifying assumptions. (i) The ‘standard collisional treatment’ which computes collision probabilities with the ‘nScollVrel’ formula assuming a small central homogeneous box, ignoring small-number effects and the small amplitude of the VMS motion. (ii) The ‘loss-cone approach’ which neglects this motion altogether but looks at the possibility of collision between the VMS and a given star with (approximate) account of the specific orbit of the star and its evolution (due to relaxation) on appropriately small time-scales. A limitation of the second approach, as implemented in the MC code is that it only works for interactions with the central object leading to disappearance of the star in one pericentre passage. In consequence, here, we can only allow mergers but not fly-by collisions.

In Fig. 24, a comparison is made of the results of both treatments for a cluster with W0= 8, N*=Np= 1.45 × 106 and RNB= 0.96 pc. Obviously, after the first few collisions, the growth of the VMS becomes much slower with the loss-cone approach, probably as a result of loss-cone orbits being rapidly depleted in the innermost region. In this situation, stellar evolution has a chance to catch up with mergers and to terminate the growth of the VMS as the interval of time between successive collisions decreases. For this particular cluster model, when stellar evolution is allowed, the runaway star still grows to more than 3000 M before it leaves the MS (models K8-16 and K8-17). This happens around t= 3 Myr but, thanks to collisional rejuvenation, the effective age of the VMS (its MS lifetime) is only 1.4 Myr. If the collision rate is so reduced after the onset of runaway growth, clusters which experience core collapse within just a little less than 3 Myr can probably not form a VMS more massive than a few hundred M before growth is interrupted by the evolution of the VMS or normal high-mass stars off the MS.

Figure 24

Growth of the runaway star. We compare the results of simulations using standard collision treatment (with sticky sphere approximation) with runs for which the runaway star is considered as a fixed central object and the loss-cone effects are taken into account following the treatment of Freitag & Benz (2002). All runs are for a cluster with W0= 8, N*=Np= 1.45 × 106 and RNB= 0.96 pc. The curves in solid line are for the normal collision treatment (K8-14 and K8-15). The dashed curves are cases for which a central fixed star of 120 M is present from the beginning (K8-11, K8-12, K8-17). To obtain the two dash–dotted curves, we introduced a central stellar object in a cluster which had evolved to deep core collapse using normal collision prescription (and no central object). The simulations with normal collision treatment were interrupted during the runaway phase (models RA095 and RA096). Stellar evolution was not taken into account except in one run with loss-cone treatment (K8-17), in which case the VMS left the MS at t= 2.8 Myr when its mass had reached 3700 M.

Figure 24

Growth of the runaway star. We compare the results of simulations using standard collision treatment (with sticky sphere approximation) with runs for which the runaway star is considered as a fixed central object and the loss-cone effects are taken into account following the treatment of Freitag & Benz (2002). All runs are for a cluster with W0= 8, N*=Np= 1.45 × 106 and RNB= 0.96 pc. The curves in solid line are for the normal collision treatment (K8-14 and K8-15). The dashed curves are cases for which a central fixed star of 120 M is present from the beginning (K8-11, K8-12, K8-17). To obtain the two dash–dotted curves, we introduced a central stellar object in a cluster which had evolved to deep core collapse using normal collision prescription (and no central object). The simulations with normal collision treatment were interrupted during the runaway phase (models RA095 and RA096). Stellar evolution was not taken into account except in one run with loss-cone treatment (K8-17), in which case the VMS left the MS at t= 2.8 Myr when its mass had reached 3700 M.

The correct VMS growth rate, assuming the same MR relation and collision physics, is very likely between what is found with our two different treatments of the central object. Which approach is the most accurate is difficult to tell. The effects of the wandering of the central object on rates of interaction with stars remain to be combined with the loss-cone theory (see Young 1977; Magorrian & Tremaine 1999; Sigurdsson 2003, for short discussions). Considering length scales pertaining to the interactions of the VMS with the cluster indicates that they take place in a complex regime. We take as an example, simulation K8-29 for which we have the best data about the central conditions in the runaway phase. The MC simulation indicates that the VMS has a wandering radius Rwand of a few 10−4 (in N-body units or pc) by the time it has grown heavier than 1000 M. This value is in agreement with a simple analysis based on the assumption of energy equipartition with the surrounding stars, MVMSVwand≈〈M*〉σv, with 〈M*〉≈ 50 M, and equal amount of kinetic and potential energy in the cuspy density background of these stars, ρ(R) ∝R−1.5. However, even though the motion of a massive object in a background of lighter particles interacting with it gravitationally obeys energy equipartition approximately (Chatterjee, Hernquist & Loeb 2002; Dorband, Hemsendorf & Merritt 2003; Laun & Merritt 2004), the energy dissipation introduced by collisions can lead to non-equipartition, so the agreement may be fortuitous. In any case, the wandering radius of the VMS is still much larger than its radius (∼50 R≃ 10−6 pc) so it cannot, in principle, be considered strictly fixed at the centre. On the other hand, the central number density (estimated inside R= 1.2 × 10−3 for 18 successive snapshots) is n0≃ 1.5 × 108, corresponding to an average distance between stars of 〈d〉=n−1/30≃ 2 × 10−3 and the influence radius of the VMS is formally Rh=GMVMS/(3σ2v) ≈ 10−3 (Frank & Rees 1976). These scales are significantly larger than Rwand, pointing to the breakdown of our standard collision treatment already when MVMS≃ 1000 M. The two approaches we have tried make strong simplifications of a situation that can probably only be treated correctly by direct orbital integration using, for instance, some hybrid numerical code in which the central region of the cluster is integrated with a direct N-body scheme, while the outer region is modelled in the FP approximation.

## 3 SUMMARY AND DISCUSSION

### 3.1 Summary of results

This work is a continuation of our study of the dynamical evolution of young, dense stellar clusters. Our goal is to determine under which conditions an IMBH can form at the centre of such systems. We are exploring the collisional runaway route. This scenario applies to clusters in which the accumulation of massive MS stars (M*≃ 10–120 M) at the centre, through the relaxation-driven process of mass segregation, leads to the Spitzer instability, thus triggering core collapse, before these stars evolve to become compact objects, that is, within ∼3 Myr. As the central density of massive stars increases and their velocity dispersion decreases, direct collisions between stars must occur, leading to the possibility of successive mergers and the runaway formation of one VMS (M*≳ 400 M). Such a VMS is a possible IMBH progenitor.

In GFR04, we have established the cluster conditions required for the core collapse to be sufficiently rapid. In the present work, we have investigated the next stage by implementing stellar collisions. We have followed the core collapse and collisional runaway using a MC simulation code specially devised to account for stellar collisions in a realistic way [me(ssy)**2]. In GFR04, two-body relaxation was the only physical process considered in most simulations so that, as expected, up to statistical fluctuations, the results did not depend on the physical size of the cluster or on the number of stars (provided there are enough for the core to be resolved). The significant free parameters were the initial structure of the cluster (profiles of density, velocity and, in some cases, mass segregation) and the IMF. Here, two more ingredients, collisions and stellar evolution, are introduced, which break this degeneracy. For instance, the ratio of the physical central velocity dispersion to the escape velocity from a star (hence N*R*/Rh) determines gravitational focusing and the outcome of collisions and, more importantly, the ratio tcoll/trlx.

We have computed more than 100 models, varying the cluster size (Rh in the range 0.03–5 pc), mass (N*= 105− 108) and concentration. For the latter parameter, we have considered King models with W0= 3 and W0= 8. For most simulations we have used a Salpeter IMF extending from 0.2 to 120 M. We treated collisions either in the sticky-sphere approximation or by making use of prescriptions for merger conditions and mass loss derived from a large set of SPH simulations. We used up to 9 × 106 particles (one simulation) but, because of the large number of simulations and the relatively large amount of CPU time required to track the runaway growth (requiring very short time-steps), we had to perform most high-N* runs using a number of particles lower than N*. Taking advantage of the statistical nature of the MC algorithm, our code scales physical processes so that each particle may represent a given (fixed) number of stars, not necessarily one. This feature may legitimately be questioned when a runaway happens as one expects only one star to experience it. However, by varying the number of particles (but keeping N* fixed) for a few cluster simulations, we have checked that using Np < N* does not appear to produce spurious results (such as, e.g. preventing the runaway).

We found that runway VMS formation occurs in all cases for which the core-collapse time is shorter than the MS lifetime of massive stars, as predicted in GFR04. Furthermore, for very massive clusters (≳107 M), the core collapse is facilitated (and sometimes driven) by collisions themselves, an effect which extends the runaway domain to clusters of larger size than predicted by an analysis based solely on relaxation. Such systems have velocity dispersions in excess of 300 km s−1 but we find that the collisional mass loss and the reduced merger cross section cannot prevent runaway, even at ∼1000 km s−1.

We never observe the formation of more than one VMS. Only in clusters that are initially collisional do collisions produce a high-mass tail in the stellar mass spectrum. Most of the mass eventually incorporated into the VMS originates in stars near the top of the IMF. Typically, for clusters that are not initially collisional, ≳90 per cent of the VMS mass is contributed by stars in the 60–120 M range. This is due to the larger cross section of massive stars and, more importantly, to the fact that collisions take place in the collapsed core, which is dominated by these stars. In such clusters, most stars that experience a collision do so only once, namely when they merge with the growing VMS. In very dense clusters where collisions are happening from the beginning of the evolution, the contribution of lighter stars is slightly higher and a significant fraction of stars that merge with the VMS have experienced at least one previous collision.

Changing the parameters of the IMF modifies the condition for rapid core collapse but has little qualitative effect on the runaway process itself. In particular, clusters with a field Kroupa IMF (Kroupa et al. 1993), which has an exponent α= 2.7 for massive stars, do experience runaway in essentially the same way as those with a Salpeter IMF (α= 2.35), despite a much smaller total number of massive stars. If the upper cut-off of the IMF is significantly lower than the assumed Mmax= 120 M, core collapse will take longer but for Mmax≲ 30 M, this is over-compensated by a longer MS lifetime. From one simulation we performed for Mmax= 10 M, we predict that the region of parameter space leading to runaway is similar to what we found for Mmax= 120 M.

We found that in most circumstances, once the runaway has started, the average time between successive mergers becomes much shorter than the MS thermal relaxation time-scale for the VMS, tKH. Hence, its structure and response to further collisions may be greatly affected by its prior collision history. Modelling the complex interplay between stellar dynamics, collisions, stellar structure and stellar evolution which should eventually determine the final mass of the VMS is presently not possible. Among the effects with potentially strong bearing on the runaway process is the depletion of loss-cone orbits that may come into play before the VMS has reached ∼1000 M as, being more massive than surrounding stars, it becomes confined to a very small central region of the cluster. Assuming a strictly fixed central object and using the approximate treatment of loss-cone physics built into our MC code (for the study of MBHs in galactic nuclei), we have found a significant reduction in the VMS growth rate (after the first 2–4 mergers) but, in the few cases considered, the VMS still grew to a few thousand M before evolving off the MS and the average time between collisions, although 3–10 times longer, remained below tKH most of the time.

Our simulations allow us to establish the conditions for the onset of runaway collisions and to study the early VMS growth. If globular clusters are formed with a concentration equal or higher than a W0= 8 King model (and if primordial binaries cannot prevent the runaway process), a significant fraction ≳20 per cent; see Fig. 1 of Paper I) may experience this phase in the first few million years of their dynamical evolution. Typically, a cluster containing ∼106 stars at this concentration needs to have a half-mass radius smaller than Rh≃ 2 pc to experience runaway, but a proto-galactic nucleus with 108 stars has to be more compact than Rh≃ 0.8 pc.

### 3.2 Discussion

Many aspects of the runaway collision scenario and how we have approached it in our numerical simulations have been discussed already. In particular, the choice of physical ingredients considered and the uncertainties and simplifications involved in their treatment have been presented in Paper I. In the results section of the present paper, we have discussed other important uncertain aspects of the process, namely the impact of collisions on the structure and growth of the VMS and the possible effects of loss-cone depletion. In this last subsection, we consider the role of interstellar gas.

#### 3.2.1 Role of residual gas in young clusters

Because the runaway scenario presented here has to operate within the first few Myr of the cluster life, it has, in fact, to be considered in the framework of cluster formation, which takes place on a similar time-scale. No only does it mean, as mentioned above, that stars below 2–3 M may still be on the pre-MS, with much larger radii than adopted here, but, more importantly, that a significant amount of residual gas may be present. Observations show that when a cluster forms, not more than 30 per cent of the gas is eventually turned into stars (Lada 1999) but that clusters like R136 or NGC 3603, as young as 1–2 Myr, are already devoid of gas (e.g. Stolte et al. 2004). In such clusters, the expulsion of gas (originally the dominant mass component) is driven by the ionizing radiation and winds of OB stars and occurs impulsively. A large fraction of the stars are then lost from the cluster and the remaining bound cluster swells quite dramatically (Kroupa, Aarseth & Hurley 2001; Boily & Kroupa 2003a,b), an event likely to terminate the collisional phase after ∼1 Myr rather than ∼3 Myr. In clusters with an escape velocity much higher than the typical expansion velocity of the ionized gas, that is, ∼10 km s−1, complete expulsion of the gas probably only occurs when the first supernova explodes.

The faster evolution due to the much smaller size of the cluster in the embedded phase counterbalances the possible reduction in the time available for collisional runaway. Assuming that, when gas is expelled, the mass of the cluster Mcl (including gas) is reduced by a factor of ∼4 and its size Rh expands ∼10 times, the true initial relaxation time (trlx∝σ3vn−1M3/2clR3/2h) should be of the order of four times shorter than what is inferred from observations of the gas-free cluster. Furthermore, the collision time (tcoll∝σvn−1M−1/2clR5/2h) may have been some 1000 times shorter, suggesting that a large fraction of clusters may actually be collisional at birth.

The standard theory of star formation, which assumes spherically symmetric accretion of gas on to the protostar, fails for stars more massive than 10 M (Yorke 2004; Bally & Zinnecker 2005). This problem led to the proposal that the truly primordial IMF may not extend past 10 M and that mergers may be – at least partially – responsible for the formation of more massive stars (Bonnell, Bate & Zinnecker 1998; Bonnell & Bate 2002; Bally & Zinnecker 2005). This would also naturally explain why these massive stars are preferentially found in the central regions of clusters that are possibly too young to have experienced relaxational mass segregation (Bonnell & Davies 1998; de Grijs et al. 2002). To occur on a sufficiently short time-scale, collisions require central stellar densities orders of magnitude higher than what is observed in those young clusters but similar to what it may have been in the embedded phase.

Realistic simulations of embedded clusters should include the role of gas, since its mass is likely to be more than twice that in stars. Especially important are the increased velocity dispersion and accretion from the ambient gas. The pioneering SPH/N-body simulations of very early cluster evolution could only handle up to a few hundred stars (Bonnell, Vine & Bate 2004; Bate & Bonnell 2005). Our MC results (run K3-20 presented in Section 2.5.1) suggest that, in much larger systems, collisions would lead to the formation of just one (very) massive star rather than a continuous spectrum spreading from 10 to ∼100 M. However, this result applies to clusters driven to core collapse by mass segregation. When the cluster is initially collisional, more stars should take part in collisions and a high-mass spectrum may develop in a more orderly way (see Fig. 13). Also, the collisional phase may be terminated before a runaway collision sequence starts through gas expulsion by massive stars.

#### 3.2.2 Gas retention in galactic nuclei

Our predictable result that evolution of the massive stars off the MS terminates the evolution towards core collapse stems from our assumption that all gas released by the stars is completely (and instantaneously) lost from the cluster. If the gas is not totally expelled, there may be no core re-expansion. The core may instead remain very dense. This opens the possibility for a longer-term collisional evolution but now there will be plenty of stellar BHs around that may destroy MS stars.

Over a long time-scale, it is possible that star formation will convert the gas into new stars (Sanders 1970; Quinlan & Shapiro 1990). All this becomes obviously quite intricate but would have to be included to treat realistically proto-galactic nuclei. Nowadays, multi-phase dynamical simulations, including interstellar gas (possibly in two phases), stars and dark matter, and the matter and energy exchanges between these components, are commonly performed to investigate the formation and evolution of galaxies (e.g. Mihos & Hernquist 1996; Abadi et al. 2003a,b; Springel & Hernquist 2003). These works are based on the use of collisionless N-body codes coupled with SPH. Including detailed gas dynamics in models of relaxing clusters is extremely challenging because the gas should possess a complex, non spherically symmetric structure, evolving on time-scales much shorter than trlx (Coker, Melia & Falcke 1999; Williams, Baker & Perry 1999; Rockefeller et al. 2004; Cuadra et al. 2005).

Note that, in clusters with very high velocity dispersions, where collisions may significantly contribute to gas production, some of it will probably be retained. Indeed this can be seen in the extreme situation of a completely disruptive collision (a rare occurrence; see Freitag & Benz 2005). In this case, the total energy of the colliding stars is their kinetic energy at (relatively) large separation minus the sum of their individual binding energy (a positive quantity, by definition). None of the stars can have a velocity higher than the cluster's escape velocity Vesc and because the binding energy is positive, the energy available per unit mass of stellar gas has to be smaller than 1/2 V2esc so at least some gas must stay in the cluster to conserve energy.

We thank Holger Baumgardt, Melvyn Davies, John Fregeau, Pavel Kroupa, Jamie Lombardi, Simon Portegies Zwart, Rainer Spurzem and Hans Zinnecker for useful discussions. MF is indebted to Pau Amaro Seoane for his extraordinary support during his stay in Heidelberg. The work of MF was funded in part by the Sonderforschungsbereich (SFB) 439 ‘Galaxies in the Young Universe’ (subproject A5) of the German Science Foundation (DFG) at the University of Heidelberg. This work was also supported by NASA ATP Grants NAG5-13236 and NNG04G176G, and NSF Grant AST-0206276 to Northwestern University.

## REFERENCES

Abadi
M. G.
Navarro
J. F.
Steinmetz
M.
Eke
V. R.
,
2003
,
ApJ
,
591
,
499
Abadi
M. G.
Navarro
J. F.
Steinmetz
M.
Eke
V. R.
,
2003
,
ApJ
,
597
,
21
Amaro-Seoane
P.
Freitag
M.
Spurzem
R.
,
2004
,
MNRAS
,
352
,
655
Bahcall
J. N.
Wolf
R. A.
,
1976
,
ApJ
,
209
,
214
Bally
J.
Zinnecker
H.
,
2005
,
AJ
,
129
,
2281
Baraffe
I.
Heger
A.
Woosley
S. E.
,
2001
,
ApJ
,
550
,
890
Bate
M. R.
Bonnell
I. A.
,
2005
,
MNRAS
,
356
,
1201
Baumgardt
H.
Heggie
D. C.
Hut
P.
Makino
J.
,
2003
,
MNRAS
,
341
,
247
Baumgardt
H.
Makino
J.
Ebisuzaki
T.
,
2004
,
ApJ
,
613
,
1133
Baumgardt
H.
Makino
J.
Ebisuzaki
T.
,
2004
,
ApJ
,
613
,
1143
Begelman
M. C.
Rees
M. J.
,
1978
,
MNRAS
,
185
,
847
Binney
J.
Tremaine
S.
,
1987
,
Galactic Dynamics
.
Princeton Univ. Press
,
Princeton, NJ
Boily
C. M.
Kroupa
P.
,
2003
,
MNRAS
,
338
,
665
Boily
C. M.
Kroupa
P.
,
2003
,
MNRAS
,
338
,
673
Bond
J. R.
Arnett
W. D.
Carr
B. J.
,
1984
,
ApJ
,
280
,
825
Bonnell
I. A.
Bate
M. R.
,
2002
,
MNRAS
,
336
,
659
Bonnell
I. A.
Bate
M. R.
Zinnecker
H.
,
1998
,
MNRAS
,
298
,
93
Bonnell
I. A.
Vine
S. G.
Bate
M. R.
,
2004
,
MNRAS
,
349
,
735
Bonnell
I. A.
Davies
M. B.
,
1998
,
MNRAS
,
295
,
691
Chatterjee
P.
Hernquist
L.
Loeb
A.
,
2002
,
ApJ
,
572
,
371
Coker
R.
Melia
F.
Falcke
H.
,
1999
,
ApJ
,
523
,
642
Colgate
S. A.
,
1967
,
ApJ
,
150
,
163
Cuadra
J.
Nayakshin
S.
Springel
V.
Di Matteo
T.
,
2005
,
MNRAS
,
360
,
L55
De Grijs
R.
Gilmore
G. F.
Mackey
A. D.
Wilkinson
M. I.
Beaulieu
S. F.
Johnson
R. A.
Santiago
B. X.
,
2002
,
MNRAS
,
337
,
597
Dorband
E. N.
Hemsendorf
M.
Merritt
D.
,
2003
,
J. Comput. Phys.
,
185
,
484
Duncan
M. J.
Shapiro
S. L.
,
1982
,
ApJ
,
253
,
921
Ebisuzaki
T.
et al
,
2001
,
ApJ
,
562
,
L19
Frank
J.
Rees
M. J.
,
1976
,
MNRAS
,
176
,
633
Freitag
M.
Benz
W.
,
2001
,
A&A
,
375
,
711
Freitag
M.
Benz
W.
,
2002
,
A&A
,
394
,
345
Freitag
M.
Benz
W.
,
2005
,
MNRAS
,
358
,
1133
Freitag
M.
Rasio
F. A.
Baumgardt
H.
,
2006
,
MNRAS
, this issue (doi: DOI: 10.1111/j.1365-2966.2006.10095.x) (Paper I)
Gürkan
M. A.
Freitag
M.
Rasio
F. A.
,
2004
,
ApJ
,
604
,
653
(GFR04)
Hénon
M.
,
1971
,
Ap&SS
,
13
,
284
Hénon
M.
,
1973
, in
Martinet
L.
Mayor
M.
, eds,
Dynamical structure and evolution of stellar systems, Lectures of the 3rd Advanced Course of the Swiss Society for Astronomy and Astrophysics
.
SSAA
,
Sauverny
, p.
183
Kroupa
P.
,
2001
,
MNRAS
,
322
,
231
Kroupa
P.
Weidner
C.
,
2003
,
ApJ
,
598
,
1076
Kroupa
P.
Tout
C. A.
Gilmore
G.
,
1993
,
MNRAS
,
262
,
545
Kroupa
P.
Aarseth
S.
Hurley
J.
,
2001
,
MNRAS
,
321
,
699
Lada
E. A.
,
1999
, in
Lada
C. J.
Kylafis
N. D.
, eds, NATO ASIC Proc. 540,
The Origin of Stars and Planetary Systems
.
Kluwer Academic Publishers
,
Dordrecht
, p.
441
Laun
F.
Merritt
D.
,
2004
,
Brownian Motion of Black Holes in Dense Nuclei
. preprint (0408029)
Lee
H. M.
,
1987
,
ApJ
,
319
,
801
Lee
M. H.
,
1993
,
ApJ
,
418
,
147
Lee
M. H.
,
2000
,
Icarus
,
143
,
74
Lightman
A. P.
Shapiro
S. L.
,
1977
,
ApJ
,
211
,
244
Lightman
A. P.
Shapiro
S. L.
,
1978
,
Rev. Mod. Phys.
,
50
,
437
Lombardi
J. C.
Warren
J. S.
Rasio
F. A.
Sills
A.
Warren
A. R.
,
2002
,
ApJ
,
568
,
939
Lombardi
J. C.
Thrall
A. P.
Deneva
J. S.
Fleming
S. W.
Grabowski
P. E.
,
2003
,
MNRAS
,
345
,
762
Magorrian
J.
Tremaine
S.
,
1999
,
MNRAS
,
309
,
447
Malyshkin
L.
Goodman
J.
,
2001
,
Icarus
,
150
,
314
McCrady
N.
Gilbert
A. M.
Graham
J. R.
,
2003
,
ApJ
,
596
,
240
McMillan
S.
Portegies Zwart
S.
,
2004
,
St-Louis
N.
Moffat
A.
, eds,
Stellar Collisions and Black Hole Formation in Dense Star Clusters. Massive Stars in Interacting Binaries
. Preprint (0412622)
Mihos
J. C.
Hernquist
L.
,
1996
,
ApJ
,
464
,
641
O'Leary
R. M.
Rasio
F. A.
Fregeau
J. M.
Ivanova
N.
O'Shaughnessy
R.
,
2006
,
ApJ
,
637
,
937
Portegies Zwart
S. F.
McMillan
S. L. W.
,
2002
,
ApJ
,
576
,
899
Portegies Zwart
S. F.
Makino
J.
McMillan
S. L. W.
Hut
P.
,
1999
,
A&A
,
348
,
117
Portegies Zwart
S. F.
Baumgardt
H.
Hut
P.
Makino
J.
McMillan
S. L. W.
,
2004
,
Nat
,
428
,
724
Quinlan
G. D.
Shapiro
S. L.
,
1990
,
ApJ
,
356
,
483
Rauch
K. P.
,
1999
,
ApJ
,
514
,
725
Rockefeller
G.
Fryer
C. L.
Melia
F.
Warren
M. S.
,
2004
,
ApJ
,
604
,
662
Sanders
R. H.
,
1970
,
ApJ
,
162
,
791
Schaerer
D.
,
2002
,
A&A
,
382
,
28
Shapiro
S. L.
,
1977
,
ApJ
,
217
,
281
Sigurdsson
S.
,
2003
,
Classical and Quantum Gravity
,
20
,
45
Spitzer
L.
,
1987
,
Dynamical Evolution of Globular Clusters. Princeton Univ. Press, Princeton

Spitzer
L. J.
Saslaw
W. C.
,
1966
,
ApJ
,
143
,
400
Springel
V.
Hernquist
L.
,
2003
,
MNRAS
,
339
,
289
Stolte
A.
Brandner
W.
Brandl
B.
Zinnecker
H.
Grebel
E. K.
,
2004
,
AJ
,
128
,
765
Vishniac
E. T.
,
1978
,
ApJ
,
223
,
986
Weidner
C.
Kroupa
P.
,
2004
,
MNRAS
,
348
,
187
Williams
R. J. R.
Baker
A. C.
Perry
J. J.
,
1999
,
MNRAS
,
310
,
913
Yorke
H. W.
,
2004
, in
Burton
M.
Jayawardhana
R.
Bourke
T.
, eds, IAU Symp. 221,
Star Formation at High Angular Resolution
.
Astron. Soc. Pac.
,
San Francisco
, p.
141
Young
P. J.
,
1977
,
ApJ
,
217
,
287
1
This had to be expected because the ratio of relaxation time to collision time is trlx/tcoll∼ (ln Λ)−1 (1 +θ)θ−2 with θ≃ (V*v)2.
2
See Freitag & Benz (2001) for comparisons of the tcc|rlx values found with different codes, in the case of the collapse of a single-mass model. Models with a broad mass spectrum are more affected by particle noise than single-mass clusters (at a given Np) and display more intrinsic and code-to-code variation.
3
In the innermost parts of the cluster, this density is computed from the volume of the spherical shell containing 17 nearby particles (Hénon 1973).

## Author notes

*
Present address: Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA.