The Origin of OB Runaway Stars
About 20% of all massive stars in the Milky Way have unusually high velocities, the origin of which has puzzled astronomers for half a century. We argue that these velocities originate from strong gravitational interactions between single stars and binaries in the centers of star clusters. The ejecting binary forms naturally during the collapse of a young ( Myr) star cluster. This model replicates the key characteristics of OB runaways in our galaxy and it explains the runaway stars around young star clusters, e.g. R136 and Westerlund 2. The high proportion and the distributions in mass and velocity of runaways in the Milky Way is reproduced if the majority of massive stars are born in dense and relatively low-mass ( ) clusters.
Most stars in our galaxy have a relatively low velocity. However, there is a population of fast moving stars, called OB runaways (?), which have considerably higher space motions of km/s (?). The origin of such velocities can be attained in two very distinct ways: A runaway can be launched when its binary companion explodes in a supernova (?), or by the ejection via a dynamical slingshot (?). The relative importance of both mechanisms has remained elusive, mainly because both are associated with young stellar populations and the high speed of a star is generally observed long after it has moved away from its birth place.
A massive star can be accelerated effectively by a three-body dynamical interaction (?). Every hard interaction eventually results in a collision between two or all three participating stars, or in the escape of one star and one binary (?). The velocity acquired by the ejected star easily exceeds the escape speed of a star cluster: For a binary with total mass and semi-major axis , the typical velocity with which the single star is ejected is .
It is typically the least massive star that is ejected (?), and both components of the ejecting (bully) binary are therefore likely to be more massive than the escaping star. For higher-mass binaries, the maximum cross section for ejecting an intruder with high speed shifts to a larger orbital period, and to a higher binding energy (Supplement § A). The binary that most effectively produces massive runaways is relatively wide, composed of the most massive stars in the cluster and located in the core. The gravothermal collapse of a cluster core naturally produces such a binary (?), which subsequently hardens by ejecting stars (?), until it experiences a collision or hardens to kT (?), upon which an encounter causes its ejection from the cluster (?)
A binary with a binding energy of 10 to 10,000 kT is therefore a natural consequence of the core collapse of a star cluster (?), and leads to the frequent ejection of other cluster members, until it eventually ejects itself or engages in a phase of Roche-lobe overflow. During the hardening process, each encounter typically drains % (?) of the binary’s binding energy. In this process it ejects on the order of stars from the cluster before the binary ejects itself or becomes semi-detached. The collapsed cluster core can support only one bully binary at a time and as a consequence the number of runaways produced does not depend on cluster mass.
We quantized these results by performing a series of -body simulations for clusters of 6,300 to (Tab. S1). For each simulation, we recorded the moment of core collapse and the time of ejection, mass and velocity of the escaping stars. Each cluster produced a single binary that ejected stars with km/s of which stars have a mass , irrespective of the mass of the cluster (Supplement § B); this number is consistent with the 6 runaways with ejected from R136 (?, ?). With a mass of and a core radius of pc (?, ?) R136 represents the population of young ( Myr), dense and massive star clusters in the local group of galaxies (Fig. S2). This star cluster is ideally suited for comparing with our simulations because of the excellent quality of the copious observations.
In our simulations, the large mass and relatively wide orbits of the binaries produced during the gravothermal collapse of a cluster core completely dominate the cross section for stellar encounters (Supplement § A). The encounter probability then should depend only weakly on the mass and size of the single stars in the cluster core. As a result, the mass distribution of the ejected stars roughly follow the mass function in the cluster core which should be relatively flat compared to the global cluster mass function (?), because equipartition of kinetic energy has resulted in the massive stars to sink to the cluster center (Supplement § B).
In Fig. 1 we compare the distribution of the masses of the runaway stars with the mass function in the cluster core (see also Supplement § B). The mass distribution of the runaway stars is statistically indistinguishable from the core distribution. For stars of the mass distribution of the runaways are flatter than the Salpeter slope, whereas it is considerably steeper for . The steep high-mass end is consistent with the recently observed distribution of massive runaway stars and field stars in the small Magellanic cloud (?).
The binaries from the -body simulations populate the area around the dashed curve in Fig. 2. These binaries came from a longer orbital period and lower mass to reach their minimum orbital period mediated by encounters and collisions (Supplement § A). Binaries to the right of this curve are in the process of dynamical hardening by ejecting stars from the cluster; they most efficiently produce runaways. When they move away from the dashed blue curve their encounter frequency drops and stellar mass loss in the winds of the stars start to dominate their evolution towards lower masses and longer orbital periods (Fig. 2).
The short day period binaries observed inside young star clusters like R136, must be primordial and probably never experienced a dynamical encounter with another cluster member (Fig. 2). From a dynamical perspective these binaries can be considered as single stars (?), because their interaction cross section is very small (Fig. S1). A strong interaction with such a tight binary tends to cause a collision between two or all three stars (?, ?); % of the collision products remain in the cluster, the others escape (Supplement § A).
The relatively wide ( year) binary HD93129AB (?, ?) (green star in Fig. 2) probably formed during the early collapse of Trumpler 14 (?, ?). With the observed mass and density of Trumpler 14 (?) this binary would be about 2,700 kT, consistent with being formed dynamically during core collapse and currently in the process of hardening, but insufficiently hard to eject itself from the cluster. We expect that HD93129AB has ejected stars, a few of which with .
The outskirts of the star cluster R136 also exhibits evidence of its violent dynamical history. The rapidly rotating 150 star VFTS 682 (?) and the single 90 star 30Dor 016 (?), ware recently demonstrated to once have belonged to the cluster. The proper motion and projected distance from R136 indicates that VFTS 682 and 30Dor 016 were ejected 29[pc]/30[km/s] Myr and 120[pc]/85[km/s] Myr ago, respectively, shortly after the birth of the Myr old cluster. By this time R136 was much too young to have experienced a supernova, challenging the favorite explanation for the origin of OB runaway stars by the supernova explosion of a stellar companion (?).
The binary R145 is sufficiently hard to be responsible for having ejected VFTS 682 and 30Dor 016 from R136, but so massive ( (?, ?)) that conservation of linear momentum would not impart a high enough velocity kick to its center of mass to eject itself from R136.
Diagonally opposite R145 with respect to the center of R136 is the x-ray source 16, which is an ACIS (?) source with a visual as well as a 2MASS counterpart (?). Conservation of linear momentum may have caused both objects to be ejected in opposite directions with respect to R136. Based on the respective projected distances of these two objects from the star cluster we predict that R145 has ejected itself and the x-ray source 16 from R136. The latter then exhibits a single star or a tight binary.
R145 and the other massive binary in Fig. 2, WR25 (HD93162) in the Carina region (?), are representative of the population formed dynamically during the core collapse of the nearby parent cluster, and flung out by a dynamical interaction. Neither are currently in the central portion of the nearby clusters, R136 and Trumpler 16, respectively. Both clusters must therefore be in a later stage of core collapse where the binary burning has stopped at the cost of the ejection of a rich population of runaways, whereas the core of Trumpler 14 collapsed only recently.
Because each cluster generates one runaway-producing binary during core collapse the relative fraction of runaways is inversely proportional to the mass of the cluster (Fig. 3). A more complete census of OB runaways can be obtained in the disk of the Milky Way in which more than 634 runaway stars have been identified (?). The mass function of these runaways is consistent with that of our simulations of young core-collapsed star clusters of 6300 (Fig. 3).
The majority of the galactic OB runaways seem to originate from star clusters that experienced core collapse within the first 1 Myr of their existence. The relative fraction of field runaways (Fig. 3), is somewhat higher than in R136. For stars with a mass this fraction is consistent with the stellar ejecta produced in relatively low mass star clusters that experience core collapse before they turn 1 Myr old, whereas the R136 results matches the simulations with a mass of . Such star clusters can experience core collapse within about a million years if they are born highly concentrated.
Supporting Online Material
Figs. S1 to S4
References (39 — 57)
In the main paper we discuss the ejection of stars through a gravitational interaction from a binary system. Here we present the more technical part of our discussion. In particular the discussion on three-body interactions in the cluster core and the macroscopic dynamics of the star cluster. This combination of microscopic gravitational dynamics together with the global evolution of the star cluster allows us to quantify the processes that lead to the ejection of massive runaway stars and at the same time qualify these processes by means of global structure evolution calculations of the stellar clusters of interest. Our focus is generally on nearby young and massive clusters. From an observational point of view these clusters are well studied, and form a pivotal role in the understanding of the formation of star clusters. We aim our analysis towards populations of young and massive clusters from Tab. 2 of (?).
The cluster R136 is of particular interest to us, because it is well studied and a rich population of runaways was recently identified to roam the neighborhood. With a mass of , a core radius of pc (?, ?) and a virial radius pc (?) the density profile of R136 fits a King (?) model with a dimension-less central potential . As a consequence the density in the cluster core /pc, which is quite typical for a core collapsed star cluster (?) (see also Fig. 5). Inside this cluster close encounters between stars occur quite frequently. Most common are the encounters between a single star and a binary. In § A we discuss the results of a series of numerical scattering experiments between a binary and a single star. We attempted to match the initial conditions for the scattering experiments to mimic the environment in the core of R136 and in particular the runaway 30Dor 016. In § B we report on N-body simulations of star clusters which bracket the range of initial masses of observed star clusters, and the cluster size is chosen such that they experience core collapse within about one million years.
Appendix A Numerical three-body scattering experiments
Here we specifically address the relation between the interaction cross section and the binary parameters in order to understand the orbital characteristics of the binaries that are likely to eject massive stars from a star cluster. The conclusions drawn from the results presented in this supplement support the claim in the main paper that relatively wide ( days) and massive binaries are most efficient in ejecting stars from the cluster, whereas tight ( day) binaries are unlikely to interact and therefore likely to be primordial. If binaries with orbital parameters similar to the observed population of short period binaries in young and massive star cluster interact dynamically they either cause a collision between stars or eject themselves from the cluster.
We determine the interaction cross section by performing a series of three-body scattering simulations using the sigma3 package within the starlab software environment (?). In sigma3 cross sections are determined by executing a number of three-body simulations in which a single star is injected into a specific binary. The encounter is initialized using the masses and radii of all the three stars, and the relative velocity between the binary center of mass and the intruder. The direction from which the intruding star approaches the binary and the binary eccentricity are selected randomly. The latter from the thermal distribution ranging from a circular orbit to a maximum eccentricity such that the stars do not touch at pericenter. The radii of the stars are calculated adopting zero-age main-sequence stars with solar composition using MESA (?) called from the AMUSE framework (?). The orbital phase and the direction from which the intruder approaches the binary are properly randomized (?). The three-body cross section is subsequently determined, as described in (?). All calculations are run with a relative energy error .
In Fig. 4 we present the results of a series of scattering experiments in which we launch a main-sequence star into a binary consisting of two stars of the same mass as the incoming star. The left panel in Fig. 4 is calculated using 16 stars with a radius of 4.9 whereas the right panel is calculated using 90 stars with a radius of 12.8 . We deliberately selected the single stars and the binary components to have the same mass to make the interpretation easier, and to prevent selection effects due to the sampling of the other binary parameters. We performed additional simulations using stellar masses of 150 200 and 400 The blue dashed curve in Fig.1 of the main paper is a least squares fit to the stellar mass and the maximum cross section for ejecting a star with a velocity of km/s, which we calculated by means of scattering experiments. This curve is consistent with a constant cross section through mass and orbital period.
The bottom x-axis and left-hand y-axis in Fig. 4 give initial orbital period of the target binary (in days) and the cross section (in units of ) resulting from the scattering experiments. Each bullet point in Fig. 4 results from a encounter density of 3000 per surface area of the binary, totaling to about scattering experiments. The various curves are linear interpolations between two bracketing points. The top (blue) curve in Fig. 4 gives the total cross section for exchanges and collisions (between two or all three stars). The green curve presents the cross section for collisions only. After the encounter is resolved, we record the velocity of the escaping star. The red bullets give the cross section for an encounter which leads to an escaper with an ejection velocity of at least km/s, which is sufficient to escape from a cluster like R136 and to be considered a runaway star. The blue (“+” signs) and magenta (“” signs) curves present the cross section for a merger product and a binary with an ejection velocity km/s, respectively.
Along the top x-axis of Fig. 4 we present the hardness of the binary (in terms of kT, where kT is the thermodynamic unit of kinetic energy in the cluster with being the total binding energy for a cluster of stars) with an orbital separation consistent with the bottom x-axis. The binding energy of the binary is calculated adopting cluster parameters consistent for R136. Along the right-hand y-axis we express the cross section (left hand y-axis) in terms of the encounter rate between the binary and any other star (in units of Myr), which we calculated by adopting a central velocity dispersion of km/s.
A binary of two 16 stars (the minimum for a main-sequence spectral type O star) with a year orbital period can already eject an incoming 16 star with a space velocity of km/s. In the core of a typical young and dense star cluster a binary with these parameters exhibits a binding energy of several 10 kT. Although such a binary is wide, the cross section for producing a runaway in an interaction is relatively small, but rapidly increases with decreasing orbital period, until a maximum is reached at a period of a few years (see Fig. 4). The maximum cross section for this binary for producing a runaway with km/s peaks at an orbital period of 800 days ( kT) with an encounter rate of Myr. For higher mass stars (90 ) the peak shifts to 2100 days for an encounter rate of Myr. These rates are consistent with our global estimate in the main paper for the orbital separation and mass of a binary that most efficiently ejects an incoming star from the cluster.
For larger orbital period (softer binaries) the cross section for producing a runaway drops. For short orbital period ( day) the cross section is dominated by collisions, which is consistent with the result presented by (?), although they adopted with stellar masses comparable to those adopted in the bottom panel of Fig. 4.
The cross section for a collision exceeds the cross section for producing a runaway over the entire range of orbital separations, and this difference becomes more prominent for wider binaries. The majority of these collision products, remain in the cluster, but a small fraction of about 1/10 escape. The fraction of merger products among the escaping stars hardly depends on the orbital separation, at least up to day (compare the red bulleted curve with the lower magenta curve in Fig. 4). The fraction of collision products among the runaways is then expected to be rather constant, in particular because this trend appears to be rather insensitive to the adopted initial conditions. For each five escaping normal stars we then expect about double that number of mergers to remain in the cluster and one merger product to escape as a runaway.
The cross section for escaping wide binaries is comparable to that for escaping merger products. This cross section does not include day primordial binaries which dynamically behave similar fashion as single stars (?).
Appendix B Simulations of massive young star clusters
To quantify the arguments concerning the orbital characteristics of the binaries that form during core collapse, their efficiency of ejecting stars with a high velocity and the distribution of stellar masses ejected, we performed a series of direct N-body simulations of star clusters over a range of cluster mass and size.
b.1 Method and Initial Conditions
In our -body simulations we varied the total cluster mass, its size and the stellar mass function in order to study the characteristics of the population of runaway stars. In Tab. 1 we present an overview of the simulations performed.
For all simulations we adopted the Salpeter initial mass function (?) between a lower limit , which we varied between 0.2 and 1 and a fixed upper limit of 100 . The upper limit and the slope of the selected power-law mass function resulted in most of the low-mass ( ) clusters to have a maximum mass that rarely exceeded 30 .
In our simulated star clusters we varied the mass, the virial radius and the concentration of the density profile, in such a way that the central density remained roughly the same. Each cluster was initialized using a King (?) model density distribution with a dimension-less central potential, , of 2 for the low-mass clusters, 6 for the intermediate mass clusters and 8 for the most massive clusters. As a consequence, the central density for all clusters pc.
In the table (Tab. 1) we present the initial conditions and main results of our model simulations.
In Fig. 5 we present the evolution of the density within the half-mass radius for one of our simulations (C, see Tab. 1) and for comparison present the initial densities for three globular clusters, and the observed densities of several nearby young and massive clusters. The half-mass densities for the three globular clusters have been derived by iteratively performed -body simulations with primordial binaries and stellar evolution until the final simulation results matches the observed star clusters (?, ?, ?). The initial density used in our model is bracketed by those derived for the three globular clusters, and is quite consistent with the observed values.
We run the simulations using a sixth-order Hermite predictor-corrector integrator (?), with a dimension-less accuracy parameter of –0.3 (?). We incorporated stellar collisions using the sticky-sphere approach for which we adopted the stars to have a radius from solar composition zero-age main-sequence stars (?). Initially our simulations have no stars with a mass , but due to collisions more massive stars started to appear shortly after core collapse. For stars of , we adopted a mass loss by stellar wind of [yr] (?), for lower mass stars we ignored the mass loss.
The simulations were curried out using Cray XT4 in National Astronomical Observatory of Japan with 64–512 cores for 16k and 64k clusters and Little Green Machine in Leiden Observatory with 8 cores for 2k clusters.
b.2 Results from the -body simulations
We ran each simulation up to an age of 3 Myr which is, for all the initial conditions presented in Tab. 1, well beyond the moment of core collapse. Core collapse in our models is a necessity to develop a mass-segregated mass function and the dynamically active bully binary in the core. For each simulation we record the moment of core collapse and the time of ejection, mass, and velocity of the escaping stars.
In Fig. 1 (of the main paper) we present the mass function of the stars in the cluster core at an age of 3 Myr. Except for the natural run-to-run variations, the range of initial conditions did not show a significant difference in the runaway or the core mass function. We demonstrate this by presenting the core mass function for all runs combined, but present the separate runaway mass function for the three models A, C and D (see Tab. 1). We compare this mass function with the distribution of stellar masses that escaped the cluster, with a Selpeter mass function and with the recently observed distribution of massive runaway and field stars in the small Magellanic cloud (?).
The distribution of relatively low mass runaway stars is consistent with the Salpeter slope. More massive stars ( ) are significantly over represented among the runaways and among the population in the cluster core. The mass function of the runaway stars with a mass is indistinguishable from the distribution of masses in the cluster center. Although the massive stars are over-represented among the runaways and field stars (see Fig. 1 of the main paper), the high-mass end ( ) distribution is considerably steeper, as was also observed in the small Magellanic cloud (?). The consistency between the observed and simulation results further support our finding that the majority of massive stars are formed in relatively small (in size as well as in mass) star clusters and that the ejection mechanism is dynamic.
As we discussed in the main paper, both distributions are statistical indistinguishable because the cross-section for gravitational interactions is dominated by the relatively wide binary formed dynamically during the collapse of the cluster core; the distribution of ejected stars depends only weakly on the mass of the single stars. The mass function of the runaway stars depends slightly on the mass of the cluster in the sense that the least massive clusters tend to have a slightly steeper mass function. This is caused by the lower probability of each individual low-mass cluster to contains a massive star.
b.3 Velocity distribution
In figure 6 we present the velocity distribution of runaway stars of our simulations with a minimum mass of 1 using 6300 and . Opting for 0.2 as a minimum mass to the initial mass function does not affect the results significantly.
The velocity distribution of ejected stars is skewed towards the cluster escape speed, but when the binary becomes harder and the ejected star less massive, higher escape velocities become more common. The number of runaway star depends on the number of binaries formed during core collapse, and not on the mass of the cluster. The distribution of stars ejected reflects the hardening process of a single tailor-made binary that formed during core collapse. Because it is a single binary that produces all runaways, the scattering process in the cluster core drives the energy and rate of escapers. The mass distribution of the escapers is then consistent with the stellar population in the cluster core. In a single cluster, only one hard binary forms, and the number of ejected stars reflects its hardening process. The number of runaways produced, their distribution in mass and in velocity then does not depend on the mass of the cluster.
The characteristic runaway velocity of massive stars, however, does depend on cluster mass, even though it is still the single binary process that drives the runaways. The more massive stars tend to have a higher velocity in massive clusters than in the lower mass clusters. This trend continues also for a more massive selection of stars, of . The lower velocities, on average, of the more massive stars in the lower-mass clusters is caused by the smaller proportion of massive stars. Ejecting a massive star from a cluster is achieved most easily by a binary consisting of massive stars. However, the lower mass clusters have a smaller probability of containing a sufficiently large population of massive stars to achieve this. The more massive cluster, therefore, are able to eject massive stars with higher speed. A similar effect was noted in relation to stellar collisions (?).
In Fig. 7 we present the space velocity as a function of stellar mass. The dark blue solid curve gives the mean runaway velocity for the simulations C and D as a function of mass in a co-moving bin of 30 stars. The horizontal bars (red) give the observed mean velocity of stars in the small Magellanic cloud (?). For spectral type O stars ( ) our simulations are quite consistent with the observed velocities, although for later type spectral types our simulations tends to somewhat under produce the velocities. Although using the Hipparcos catalogue (?) have derived a mean velocity of km/s for runaways with a mean mass of . The average runaway velocity for stars with seem to be systematically higher in the SMC than those observed in the Galaxy and our simulations.
References and Notes
- 1. Blaauw, A. & Morgan, W. W. The Space Motions of AE Aurigae and μ Columbae with Respect to the Orion Nebula. ApJ 119, 625 (1954).
- 2. Gies, D. R. & Bolton, C. T. The binary frequency and origin of the OB runaway stars. ApJS 61, 419–454 (1986).
- 3. Blaauw, A. On the origin of the O- and B-type stars with high velocities (the ”run-away” stars), and some related problems. Bul. Astron. Inst. Neth. 15, 265–+ (1961).
- 4. Leonard, P. J. T. & Duncan, M. J. Runaway stars from young star clusters containing initial binaries. I - Equal-mass, equal-energy binaries. AJ 96, 222–232 (1988).
- 5. Heggie, D. C. Binary evolution in stellar dynamics. MNRAS 173, 729–787 (1975).
- 6. Tanikawa, A., Hut, P. & Makino, J. Unexpected Formation Modes of the First Hard Binary in Core Collapse. ArXiv e-prints 1107.3866, (2011).
- 7. We adopted kT as the thermodynamic unit of kinetic energy in the cluster; the binding energy for a cluster of stars is kT.
- 8. Makino, J. & Hut, P. On core collapse. ApJ 383, 181–191 (1991).
- 9. Portegies Zwart, S. F. & McMillan, S. L. W. Black Hole Mergers in the Universe. ApJL 528, L17–L20 (2000).
- 10. Tanikawa, A. & Fukushige, T. Effects of Hardness of Primordial Binaries on the Evolution of Star Clusters. Publ. Astr. Soc. Japan 61, 721– (2009).
- 11. Gvaramadze, V. V., Kroupa, P. & Pflamm-Altenburg, J. Massive runaway stars in the Large Magellanic Cloud. A&A 519, A33+ (2010).
- 12. Bestenlehner, J. M. et al. The VLT-FLAMES Tarantula Survey. III. A very massive star in apparent isolation from the massive cluster R136. A&A 530, L14+ (2011).
- 13. Andersen, M. et al. The Low-Mass Initial Mass Function in the 30 Doradus Starburst Cluster. ApJ 707, 1347–1360 (2009).
- 14. Portegies Zwart, S. F., McMillan, S. L. W. & Gieles, M. Young Massive Star Clusters. ARA&A 48, 431–493 (2010).
- 15. Fregeau, J. M., Joshi, K. J., Portegies Zwart, S. F. & Rasio, F. A. Mass Segregation in Globular Clusters. ApJ 570, 171–183 (2002).
- 16. Lamb, J. B., Oey, M. S., Graus, A. S. & Segura-Cox, D. M. The IMF of Field OB Stars in the Small Magellanic Cloud. ArXiv e-prints, 1109.6655 (2011).
- 17. Portegies Zwart, S., Gaburov, E., Chen, H.-C. & Gürkan, M. A. The present day mass function in the central region of the Arches cluster. MNRAS 378, L29–L33 (2007).
- 18. Fujii, M., Iwasawa, M., Funato, Y. & Makino, J. Evolution of Star Clusters near the Galactic Center: Fully Self-Consistent N-Body Simulations. ApJ 686, 1082–1093 (2008).
- 19. Massey, P., Penny, L. R. & Vukovich, J. Orbits of Four Very Massive Binaries in the R136 Cluster. ApJ 565, 982–993 (2002).
- 20. Schnurr, O., Casoli, J., Chené, A.-N., Moffat, A. F. J. & St-Louis, N. The very massive binary NGC 3603-A1. MNRAS 389, L38–L42 (2008).
- 21. Schnurr, O., Moffat, A. F. J., Villar-Sbaffi, A., St-Louis, N. & Morrell, N. I. A first orbital solution for the very massive 30 Dor main-sequence WN6h+O binary R145. MNRAS 395, 823–836 (2009).
- 22. Taylor, W. D. et al. The VLT-FLAMES Tarantula Survey. II. R139 revealed as a massive binary system. A&A 530, L10+ (2011).
- 23. Moffat, A. F. J. Massive Binaries. In IAU Symposium (ed. F. Bresolin, P. A. Crowther, & J. Puls), vol. 250 of IAU Symposium, 119–132 (2008).
- 24. Koumpia, E. & Bonanos, A. Z. Fundamental Parameters of 4 Massive Eclipsing Binaries in Westerlund 1. ArXiv e-prints, 1009.4709 (2010).
- 25. Piatti, A. E., Bica, E. & Claria, J. J. Fundamental parameters of the highly reddened young open clusters Westerlund 1 and 2. A&AS 127, 423–432 (1998).
- 26. Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J. & Helsel, J. W. The High Angular Resolution Multiplicity of Massive Stars. AJ 137, 3358–3377 (2009).
- 27. Chené, A.-N., Schnurr, O., Crowther, P. A., Fernández-Lajús, E. & Moffat, A. F. J. Very massive binaries in R 136. In IAU Symposium (ed. C. Neiner, G. Wade, G. Meynet, & G. Peters), vol. 272 of IAU Symposium, 497–498 (2011).
- 28. Fregeau, J. M., Cheung, P., Portegies Zwart, S. F. & Rasio, F. A. Stellar collisions during binary-binary and binary-single star interactions. MNRAS 352, 1–19 (2004).
- 29. Gvaramadze, V. V. & Gualandris, A. Very massive runaway stars from three-body encounters. MNRAS 410, 304–312 (2011).
- 30. Nelan, E. P. et al. Resolving OB Systems in the Carina Nebula with the Hubble Space Telescope Fine Guidance Sensor. AJ 128, 323–329 (2004).
- 31. Sana, H. et al. A MAD view of Trumpler 14. A&A 515, A26+ (2010).
- 32. Maiz-Apellaniz, J. A Complete Multiplicity Survey of Galactic O2/O3/O3.5 Stars with ACS. In HST Proposal, 10602–+ (2005).
- 33. Ortolani, S., Bonatto, C., Bica, E., Momany, Y. & Barbuy, B. The embedded cluster DBSB 48 in the nebula Hoffleit 18: Comparison with Trumpler 14. New Astron. 13, 508–518 (2008).
- 34. Evans, C. J. et al. A Massive Runaway Star from 30 Doradus. ApJL 715, L74–L79 (2010).
- 35. But see (?) according two whom .
- 36. ACIS stands for AXAF CCD Imaging Spectrometer, see http://acis.mit.edu/acis.
- 37. Townsley, L. K., Broos, P. S., Feigelson, E. D., Garmire, G. P. & Getman, K. V. A Chandra ACIS Study of 30 Doradus. II. X-Ray Point Sources in the Massive Star Cluster R136 and Beyond. AJ 131, 2164–2184 (2006).
- 38. van der Hucht, K. A. et al. XMM-Newton Studies of the Wolf-Rayet Colliding-Wind Binaries WR 25 (WN6h+O4f) and WR 11 (WC8+O7.5III). In Massive Stars in Interactive Binaries (ed. N. St.-Louis & A. F. J. Moffat), vol. 367 of Astronomical Society of the Pacific Conference Series, 159–+ (2007).
- 39. Mackey, A. D. & Gilmore, G. F. Surface brightness profiles and structural parameters for 53 rich stellar clusters in the Large Magellanic Cloud. MNRAS 338, 85–119 (2003).
- 40. King, I. R. The structure of star clusters. III. Some simple dvriamical models. AJ 71, 64–75 (1966).
- 41. Harris, W. E. Globular Clusters in the Milky Way (Harris, 1996). VizieR Online Data Catalog 7195, 0–+ (1996).
- 42. Portegies Zwart, S. F., McMillan, S. L. W., Hut, P. & Makino, J. Star cluster ecology - IV. Dissection of an open star cluster: photometry. MNRAS 321, 199–226 (2001).
- 43. Paxton, B. EZ to Evolve ZAMS Stars: A Program Derived from Eggleton’s Stellar Evolution Code. PASP 116, 699–701 (2004).
- 44. Portegies Zwart, S. et al. A multiphysics and multiscale software environment for modeling astrophysical systems. New Astronomy 14, 369–378 (2009).
- 45. Hut, P. & Bahcall, J. N. Binary-single star scattering. I - Numerical experiments for equal masses. ApJ 268, 319–341 (1983).
- 46. McMillan, S. L. W. & Hut, P. Binary–Single-Star Scattering. VI. Automatic Determination of Interaction Cross Sections. ApJ 467, 348 (1996).
- 47. Salpeter, E. E. The Luminosity Function and Stellar Evolution. ApJ 121, 161 (1955).
- 48. Heggie, D. C. & Giersz, M. Monte Carlo simulations of star clusters - V. The globular cluster M4. MNRAS 389, 1858–1870 (2008).
- 49. Giersz, M. & Heggie, D. C. Monte Carlo simulations of star clusters - VI. The globular cluster NGC 6397. MNRAS 395, 1173–1183 (2009).
- 50. Giersz, M. & Heggie, D. C. Monte Carlo simulations of star clusters - VII. The globular cluster 47 Tuc. MNRAS 410, 2698–2713 (2011).
- 51. Nitadori, K. & Makino, J. Sixth- and eighth-order Hermite integrator for N-body simulations. New Astronomy 13, 498–507 (2008).
- 52. Hurley, J. R., Pols, O. R. & Tout, C. A. Comprehensive analytic formulae for stellar evolution as a function of mass and metallicity. MNRAS 315, 543–569 (2000).
- 53. Fujii, M., Iwasawa, M., Funato, Y. & Makino, J. Trojan Stars in the Galactic Center. ApJ 695, 1421–1429 (2009).
- 54. Gaburov, E., Gualandris, A. & Portegies Zwart, S. On the onset of runaway stellar collisions in dense star clusters - I. Dynamics of the first collision. MNRAS 384, 376–385 (2008).
- 55. Evans, C. J. & Howarth, I. D. Kinematics of massive stars in the Small Magellanic Cloud. MNRAS 386, 826–834 (2008).
- 56. Tetzlaff, N., Neuhäuser, R. & Hohle, M. M. A catalogue of young runaway Hipparcos stars within 3 kpc from the Sun. MNRAS 410, 190–200 (2011).
- 57. Pfalzner, S. Universality of young cluster sequences. A&A 498, L37–L40 (2009).
We thank B. Brandl, A. Brown, A. Gürkan and H. Sana for discussions, and the anonymous referees for their constructive criticism. This work was supported by Research Fellowships of the Japan Society for the Promotion of Science (JSPS) for Young Scientists and the Netherlands Research Council NWO (grants VICI [#639.073.803], AMUSE [#614.061.608] and Little Green Machine) and by the Netherlands Research School for Astronomy (NOVA). The simulation software is available at http://amusecode.org and the data can be found at http://initialconditions.org.