RUNAWAY MASSIVE STARS FROM R136: VFTS 682 IS VERY LIKELY A "SLOW RUNAWAY"

Fast-moving massive field stars in our Galaxy and in external galaxies, which apparently are not associated with any nearby stellar assemblies, have been a focal topic for the past two decades. Such runaway OB stars are widely thought to be members of star clusters that are ejected following dynamical encounters in their parent clusters. If the proper motions of these massive stars can be reliably estimated, the vast majority of them can then be traced back to their parent clusters. Using this kinematic method, Schilbach & Röser (2008) have demonstrated that most of the Galactic OB runaway stars indeed originated from star clusters. Another way to determine the parents of runaway massive stars is to image their bow shocks, if present, which has been possible with modern infrared telescopes (Gvaramadze & Bomans 2008; Gvaramadze et al. 2011a). A detectable bow shock is generated if a star moves supersonically (typically ≳ 10 km s⁻¹) through a dense enough medium. The geometry of a bow shock allows one to determine the direction of motion of the shock-generating star and hence its possible parent cluster (Gvaramadze et al. 2010, 2011b). Bow shocks are particularly useful when the proper motions of the runaway OB stars are unavailable or are poorly known, which is often the case for distant stars. However, even if every OB field star has been ejected from a young cluster, 1%–4% of the OB stars cannot be traced back to their parent clusters because of the two-step ejection process (Pflamm-Altenburg & Kroupa 2010): a massive binary is first ejected from its cluster and when the primary explodes as a supernova, the secondary OB star continues on a diverted trajectory.

High-velocity single and binary stars are launched from a star cluster due to super-elastic single-star–binary and binary–binary dynamical encounters (Heggie 1975) that occur predominantly in the cluster's dense core. Therefore, the spectrum of the ejected bodies is governed by the properties of the primordial binaries in the cluster (Leonard & Duncan 1988; Clarke & Pringle 1992) and its stellar initial mass function (IMF), and hence serves as a useful tracer of these fundamentally important cluster properties. As the stars and the binaries are mostly ejected from the cluster's core, where the most massive members remain segregated (mass stratification; see Spitzer 1987), the runaway population contains a large fraction of OB single stars/binaries. The ejection of the massive members also alters the stellar mass function (MF) of the bound cluster as it evolves.

A particular runaway very massive star (VMS) that has recently drawn significant attention is VLT-FLAMES Tarantula Survey (VFTS) 682. This VMS, visible in the Tarantula Nebula (30 Doradus) of the Large Magellanic Cloud (LMC), has an inferred present-day mass of ≈150 M_☉ and is located at a projected distance of ≈30 pc in the northeast from the central massive star cluster R136 of the 30 Doradus complex (Bestenlehner et al. 2011). From radial velocity measurements and its projected separation from R136, Bestenlehner et al. (2011) inferred that the true velocity of VFTS 682 should be ≈40 km s⁻¹ if it is indeed a runaway from R136. On the other hand, the apparent isolation of VFTS 682 from any star cluster, its moderate peculiar motion, and the non-detection of a bow shock (but see Gvaramadze et al. 2010; Gvaramadze & Gualandris 2011) might also indicate that this VMS has formed isolated and is unrelated to R136, as Bestenlehner et al. (2011) argue. The origin of VFTS 682 is therefore currently unclear, i.e., whether it is a "slow runaway" from R136 or has formed alone, outside any clusters. This, in turn, freshens up the open question whether massive stars form only in dense environments (Bonnell et al. 2004) or whether isolated massive star formation is possible through pure circumstellar disk accretion (Kuiper et al. 2010). Notably, the stellar content of R136 itself presents a challenge to the current understanding of star formation phenomenology as it contains several single-star members with inferred initial masses up to ≈320 M_☉, i.e., well above the widely accepted 150 M_☉ canonical upper limit of the stellar IMF (Crowther et al. 2010). As coined by Kroupa et al. (2011), the stellar population of R136 is therefore "super-saturated" by containing super-canonical stars.

Because of the important implications in massive star formation theory, it is essential to determine the likeliness of a VFTS 682-like runaway from R136. To that end, we compute the evolution of model star clusters with properties resembling R136 using the direct N-body integration method and focus on the kinematics of the massive stars ejected from them. We find that the ejection of VMSs, within a few Myr, with kinematic properties similar to VFTS 682 is quite plausible: the latter star is then quite possibly a runaway from R136 and the substantially controversial assumption of isolated massive star formation is unnecessary to explain its existence. These computations also demonstrate that super-canonical stars appear naturally from merged binaries.

In Section 2, we describe our model clusters (Section 2.1) and the method of our computations (Section 2.2). In Section 3, we present our results focusing on VMS runaway members that resemble VFTS 682 (Section 3.1) and also on 30 Dor 016-like runaways (Section 3.2). We also discuss how the ejected fraction of stars depends on the stellar mass (Section 3.3). We conclude this paper in Section 4 with a discussion and by highlighting possible future improvements.

2.1. Initial Conditions

As our primary objective is to study the runaway massive stars from an R136-like cluster, we begin our N-body computations with star clusters modeled as Plummer spheres (Kroupa 2008) having properties conforming with the observed parameters of R136. The present-day parameters of R136 still remain ambiguous. We take the total initial mass of each of our model Plummer spheres to be M_cl(0) ≈ 10⁵ M_☉, which is an upper limit of the mass of R136 (Crowther et al. 2010). As for its size, a variety of half-mass radii/core radii has so far been assigned to this cluster by different authors (Hunter et al. 1995; Portegies Zwart et al. 2002; Mackey & Gilmore 2003; Andersen et al. 2009; Portegies Zwart et al. 2010). We take the initial half-mass radii of our models to be r_h(0) ≈ 0.8 pc. The core radii of our models then turn out to be r_c ≲ 0.3 pc throughout their evolutions, 0.3 pc being an observed upper limit (Mackey & Gilmore 2003). The core density is ≳ 1.3 × 10⁴ stars pc⁻³.

The initial masses m_s of the clusters' member stars are drawn from the canonical IMF (Kroupa 2001) in the range 0.08 M_☉ < m_s < 150 M_☉. This IMF is a two-part power law with indices α₁ = 1.3 for 0.08 M_☉ < m_s < 0.5 M_☉ and α₂ = 2.3 for more massive stars. Throughout this paper, we denote the (instantaneous) mass of a stellar member by m_s in general and M is exclusively reserved to denote the mass of a runaway (see Section 3) member. The metallicity of the stars are taken to be low: Z = 0.5 Z_☉, which is appropriate for R136.

Most stars form in binaries or higher-order multiple systems. Furthermore, close binary–binary encounters play an important role in the dynamical ejection of stars from a cluster (Leonard & Duncan 1988, 1990). However, at the current state of the art, binaries are computationally the most time-consuming members in a direct N-body calculation of a star cluster because of the particular numerical algorithms involved (see Section 2.2). Thus, in order to achieve a reasonable calculation speed, we initially arrange all the stars more massive than 5 M_☉ in binaries and all the less massive members are kept single. As our primary aim is to probe the massive runaway stars, which can be efficiently ejected only via close encounters with massive, hard binaries (see Section 3.1), only massive binaries are important in this study.

To create the initial binary population, we first sort the m_s > 5 M_☉ stars with decreasing mass and then pair them in order, so that the resulting binaries have their mass ratios close to unity. It is already understood that a mass-ratio distribution biased toward unity is required in order to get a significant ejection rate of massive stars (Clarke & Pringle 1992). Furthermore, observations indicate that O-stars favor OB stars as their companions (Sana & Evans 2011). The orbital periods, P, of the binaries with primary masses >20 M_☉ are chosen from a uniform distribution between 0.5 < log₁₀P < 4 (Öpik's law; P in days) since the O-type stars are preferentially found in short-period systems. This range of P is similar to that obtained by Sana & Evans (2011).

For 5 M_☉ < m_s < 20 M_☉, Equation (8) of Kroupa (1995b) is chosen as the period distribution, which is given by

$$\begin{eqnarray*} &&f_{P,\rm birth} = \eta \frac{\left(\log _{10}P - \log _{10}P_{\rm min}\right)}{\delta + \left(\log _{10}P - \log _{10}P_{\rm min}\right)^2}, \nonumber \end{eqnarray*}$$

with η = 2.5, δ = 45, and which covers a much wider range of the period between 1.0 < log₁₀P < 8.43 (Kroupa 1995a, 1995b; Marks et al. 2011). Here, f_{P, birth} is the "birth period distribution" of binaries as obtained by Kroupa (1995a) through "inverse dynamical population synthesis." It is the primordial binary period distribution unmodified by the dynamical destructions of these binaries (i.e., without a depletion function; Kroupa 1995a) or by the mutual interactions between the binary members in their pre-main-sequence stage (i.e., without an "eigenevolution"; Kroupa 1995b). The modification of a primordial binary's orbital period, mass-ratio, and eccentricity following its eigenevolution is known only for low-mass stars (Kroupa & Petr-Gotzens 2011; Marks et al. 2011). We choose the orbital eccentricities e following a thermal distribution, i.e., f(e) = 2e (Spitzer 1987; Kroupa 2008). Figure 1 shows the resulting distribution of binary binding energies as a function of the primary mass. The two distinct energy distributions across m_s = 20 M_☉ are clearly visible; the binaries with m_s > 20 M_☉ primaries are harder than those with m_s < 20 M_☉ primaries.

Zoom In Zoom Out Reset image size

The initial systems are completely mass-segregated (Spitzer 1987) using a method by Baumgardt et al. (2008) so that the heaviest stars are initially concentrated within the clusters' cores. Such an initial condition mimics primordial mass segregation that is inferred to be true for several Galactic globular clusters (Baumgardt et al. 2008; Marks & Kroupa 2010; Zonoozi et al. 2011; Hasan & Hasan 2011) and young star clusters (Littlefair et al. 2003; Chen et al. 2007). Notably, the initial setup employed here (all massive stars in binaries, fully mass-segregated cluster) are the very first ones ever computed under such extreme conditions.

2.2. N-body Integration, Stellar Evolution, and Stellar Hydrodynamics

From the above computations, we trace the bodies that are ejected from the clusters during their evolution (within ≈3 Myr). We consider a single-star/binary/multiplet to be a runaway member from its host cluster if it is found moving away from the cluster beyond a R > 10 pc distance from the cluster's center of density. Although our primary focus is on the runaway VMSs, we consider the whole mass spectrum of ejected stars as well.

Figure 2 shows the projected snapshots of the runaways with masses M > 3 M_☉, combined from the four computations, at t = 1 Myr and 3 Myr evolutionary times. It can be seen that by t = 3 Myr there are a significant number of fast runaway VMSs with total (three-dimensional) velocities up to ≈300 km s⁻¹ and a few fast VMSs are already present at t = 1 Myr. All of these runaways are on their MSs and hence are OB stars. We note that the vast majority of the massive ejected members are single stars—only two of the massive ejecta from our computations are found in hard binaries. We shall discuss the multiplicity properties of the ejected stellar population in detail in a future paper.

Zoom In Zoom Out Reset image size

It is worthwhile to note that although our adapted canonical IMF has a 150 M_☉ upper limit, our models yield single-star runaways with masses up to ≈250 M_☉ within t < 3 Myr, considerably exceeding the widely accepted 150 M_☉ limit. A few of such members are also found bound to each of our model clusters within the same evolutionary period. These VMSs form when massive binary components merge as a result of encounter hardening (Heggie 1975; Banerjee & Ghosh 2006) of these binaries and/or eccentricities induced to them by close encounters. It is the adapted relatively small and narrow orbital period distribution of the massive binaries (m_s > 20 M_☉), conforming with the observed properties of O-type stellar binaries (see Section 2.1), that makes such mergers probable, leading to the formation of single stars with m_s > 150 M_☉. On the other hand, it is these hard, massive binaries that can efficiently drive the VMS runaways (see Section 3.1). In other words, the possibility of the presence of runaway VMSs naturally leads to the formation of VMSs (bound or ejected), in the course of the evolution of the cluster, with masses significantly exceeding the canonical upper limit, even though the cluster's IMF is intrinsically canonical. Therefore, the observed super-canonical VMSs in R136, with inferred individual initial masses up to ≈320 M_☉, do not necessarily imply an IMF with an upper-limit greatly exceeding the canonical limit as argued by Crowther et al. (2010): the m_s > 150 M_☉ single-star members can as well be accounted for as recent massive binary merger products.

Each of our computed models yields a few of such super-canonical single-star members. Because of their substantial gravitational focusing (Spitzer 1987), they are involved in very frequent close encounters with other massive binaries, thereby being vulnerable to being ejected from the cluster. Three of our four computed clusters have produced super-canonical single-star runaways with M > 150 M_☉ and the most massive merger product has been ejected from two of the models.

3.1. VFTS 682: A Very Likely Runaway from R136

To have an understanding of the spectrum of velocities with which the stars run away, we plot their three-dimensional velocities V as a function of their instantaneous masses M at t = 1 Myr and 3 Myr as shown in Figure 3. The data points from all the computations have been compiled in this figure, which are distinguished by the different symbols. We note two trends in this plot, viz., (1) there is a lower boundary of the scatter in V that increases moderately with M; and (2) the upper boundary of this scatter is independent of M. These trends are vividly apparent from the plot at t = 3 Myr (Figure 3, lower panel) where there is a significantly larger number of runaways. No well-defined boundaries can, of course, be drawn in this plot and the above-mentioned boundaries are meant to indicate net trends.

Zoom In Zoom Out Reset image size

To understand these trends, we first note that an ejected single star of mass M and velocity V can appear in two main ways, viz., (1) it is ejected due to a kinetic energy (KE) boost in a close, super-elastic encounter (Heggie 1975) with a hard binary where no exchange of members has occurred (a flyby encounter) and (2) it is ejected from a hard binary being replaced by a more massive star (an exchange encounter). In both cases, the KE of ejection is of the order of the binary's binding energy, on average. The trend of the lower boundary can be understood by noting that a single star is likely to experience a flyby in a strong encounter with a hard binary (most of them have nearly equal mass components in our model) only if the binary components are more massive than the intruder, since otherwise the latter is likely to get trapped in the binary by exchanging with one of its (lighter) members. In other words, to recoil a mass M in a flyby encounter, one needs a binary with components of mass m_s ≳ M. Similarly, a star of mass M can preferentially be ejected from a binary only by being replaced by an intruder of mass m_s ≳ M. Since the ejected star, in both cases, will typically have a KE of the order of the binary's binding energy (Heggie 1975; Heggie et al. 1996), we have MV² ∼ m²_s/a, where a is the binary's semimajor axis. Hence, the condition m_s ≳ M implies V² ≳ M/a; the smallest possible ejection velocity increases with M, for a given a. Hence, the widest (hard) binaries in our model give rise to the lower boundary of the scatter in Figure 3 (lower panel). Indeed, the eye-estimated lower limit of the data for the lower mass ejecta can be represented well by a curve V∝M^1/2 shown by the pink line segment in Figure 3 (lower panel). This is also true for the massive end; their lower boundary can be represented by another curve (the green line segment) but with the same proportionality. The discontinuity is expected due to that in the binary period distribution (cf., Section 2.1; Figure 1).

On the other hand, the upper boundary of the plot is formed by the ejecta generated by binaries having the highest binding energies. Figure 1 shows that the hardest binaries are also the most massive ones, as a result of our chosen primordial binary distribution (see Section 2.1). They can therefore eject stars over the entire mass range with KE of the order of their binding energy forming the upper boundary that is independent of M. Note that as an order-of-magnitude estimate, the highest binding energies are ∼10⁶–10⁷ M_☉ pc² Myr⁻² which yield V ∼ 10^2.5–10³ km s⁻². This indeed agrees with the highest Vs in Figure 3 in terms of order of magnitude.

In Figure 3 (lower panel) it can be seen that the median of the ejection velocities (thick black line) tend toward V ≈ 50 km s⁻¹ for the most massive runaways and the corresponding quartiles (thin black lines) span between about 40–60 km s⁻¹. Taking into account the outliers with V smaller than the first quartile in the massive end of the M–V plot, these ejection velocities are similar to the inferred true velocity of VFTS 682 (Bestenlehner et al. 2011). This result, therefore, already strongly indicates that VFTS 682 is a typical runaway VMS from R136. The highest V reached in our computations, however, exceeds 300 km s⁻¹ as can be read from Figure 3.

We further inspect that all of our computations eject VMSs with kinematic properties agreeing fairly with those of VFTS 682, viz., V ≈ 40 km s⁻¹, R ≈ 30 pc, and M ≈ 150 M_☉. Such instances are shown in Table 1. These results dictate that the VFTS 682 is an expected, very likely runaway VMS from R136.

Table 1. A List of Runaway Single Stars with Properties Fairly Close to Those of VFTS 682, as Obtained from Our Computations

Model Number	Time t	Mass M	Distance R	Velocity V
	(Myr)	( M☉)	(pc)	(km s−1)
1	2.8	256.4	31.9	27.5
	3.2	135.9	26.6	34.8
2	2.6	126.4	27.7	45.7
	2.6	125.9	29.9	49.4
3	2.6	106.9	45.7	27.3
4	1.9	169.1	29.3	29.0
	1.9	116.9	35.2	32.8
VFTS 682	<3.0	≈150.0	≈30.0	≈40.0

Notes. The columns are as follows. Column 1: model number; Column 2: evolutionary time t at which the runaway is detected; Column 3: mass M of the runaway at t; Column 4: its distance R from the cluster's center of density; and Column 5: its three-dimensional velocity V. M, R, and V of these runaway stars agree fairly with those estimated observationally for VFTS 682 (last line; Bestenlehner et al. 2011).

Download table as:  ASCIITypeset image

As additional information to the reader and for convenience of comparison with observations, we provide the data corresponding to Figure 3 in the form of an online table (Table 2). In addition to the three-dimensional velocities and the (instantaneous) masses of all the runaways from all of our computations, we also provide the corresponding line of sight and tangential velocities and projected positions. For the aid of the reader in performing further analyses, we also provide the corresponding components of the positions and the velocities in an additional online table (Table 3).

Table 2. Data for All the Runaway Stars as Obtained from All of Our Computations at 1 Myr and 3 Myr

M	vt	vrad	V	r	Single/Binary	I
(M☉)	(km s−1)a	(km s−1)b	(km s−1)	(pc)
5.9	60.0	−25.9	65.3	15.3	1	1626
8.6	10.4	10.7	15.0	9.1	1	960
0.2	8.4	8.8	12.2	7.9	1	4491
0.4	15.8	44.8	47.5	5.8	1	6332
0.3	6.8	3.3	7.6	9.6	1	8042
40.5	6.5	−8.2	10.4	19.5	2	170968
94.0	50.3	209.9	215.8	83.7	1	10
236.3	31.3	−72.8	79.2	9.4	1	5

Notes. The columns are as follows. Column 1: mass M of the runaway at 1 Myr/3 Myr; Column 2: its tangential velocity v_t (all velocities are relative to the cluster's center of mass; see footnote "a"); Column 3: line of sight or radial velocity v_rad (see footnote "b"); Column 4: three-dimensional velocity V; Column 5: its projected distance r from the cluster's center of density; Column 6: 1⇒ single, 2⇒ binary; and Column 7: identity of the star I. ^aThe reference frame is chosen fixed and originated at the cluster's center of mass. Its axes are taken to be oriented arbitrarily as there is no preferred directionality due to the absence of an external field. All the projections are taken on the x–y plane. ^bThe velocity component normal to the plane of projection or the z velocity component is simply taken to be the radial velocity due to the arbitrariness in the choice of the coordinate axes. The corrections due to the spatial extent are negligible for the distance of R136.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 3. Additional Data for All the Runaway Stars as Obtained from All of Our Computations at 1 Myr and 3 Myr

M	I	x	y	z	vx	vy	vz	Single/Binary
(M☉)		(pc)a	(pc)a	(pc)a	(km s−1)	(km s−1)	(km s−1)
5.9	1626	−12.6	−8.7	−6.5	−49.7	−33.6	−25.9	1
8.6	960	−6.9	5.9	9.0	−8.0	6.7	10.7	1
0.2	4491	4.4	6.6	8.8	3.4	7.6	8.8	1
0.4	6332	−2.2	5.3	16.4	−6.5	14.5	44.8	1
0.3	8042	−6.8	−6.7	4.5	−4.7	−4.9	3.3	1
40.5	170968	−16.1	−11.1	−25.2	−5.4	−3.5	−8.2	2
107.2	170726	10.7	−55.2	47.6	6.3	−32.7	28.2	2
236.3	5	−8.8	−3.4	−94.9	−7.8	−30.3	−72.8	1

Notes. The columns are as follows. Column 1: mass M of the runaway at 1 Myr/3 Myr; Column 2: its identity I; Columns 3–5: components of its position relative to the cluster's center of density (see footnote "a"); Columns 6–8: components of its velocity relative to the cluster's center of mass; Column 9: 1⇒ single, 2⇒ binary. ^aThe reference frame is chosen fixed and originated at the cluster's center of mass. Its axes are taken to be oriented arbitrarily as there is no preferred directionality due to the absence of an external field.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

3.2. Another Runaway: 30 Dor 016

3.3. Mass Dependence of Runaway Stars

It is very interesting to study the mass (or spectral class) dependence of the runaway members as it directly points to the pure dynamical nature of the ejection process. Although the number of stars bound to the cluster decreases strongly with increasing stellar mass, more massive stars and binaries are more centrally concentrated due to mass segregation and hence interact with each other more efficiently and frequently. This causes the runaway fraction of stars (defined as the ratio of the number of stars, within a bin around mass M, moving away from the cluster at R > 10 pc to the total number of stars in the whole system, i.e., including all the bound and the ejected members, within the same mass bin) to increase with the stellar mass as shown in Figure 4, where the outcomes of all the computations are overplotted at t = 1 Myr and 3 Myr. In Figure 4, the ejected fraction g(M) increases considerably with M for M ≳ 5 M_☉, but remains nearly flat and small for lower M. This is probably an artifact of our initial binary population where only the stars with m_s > 5 M_☉ are in binaries and therefore are efficient in ejecting only those members which are more massive than 5 M_☉. Current technology does not allow us to perform direct N-body integrations of R136-type clusters in which all stars are initially in binaries. Such computations are available only for moderate mass clusters (Kroupa 1998; Kroupa et al. 2001).

Zoom In Zoom Out Reset image size

In the present work, we have computed the dynamical evolution of model clusters whose structure conforms with the observed global properties of R136—the central massive cluster in the 30 Dor complex of the LMC, using the direct N-body integration method. The evolution of the individual stars, chosen initially from the canonical IMF with the standard 150 M_☉ upper cutoff (Weidner & Kroupa 2004), has also been incorporated as well as the evolution of the individual binaries. We focus on the ejection of massive stars from our model clusters, which is a process that depends crucially on the properties of the primordial binaries in the cluster. For computational simplicity, we have taken only the stars with initial masses m_s > 5 M_☉ to be in binaries (see Section 2.1), which are the only ones that are efficient in ejecting the massive stars. Hence, the properties of the massive runaways are not expected to be affected significantly by the absence of lower mass binaries.

Recent observations of the R136 and the 30 Dor region have raised fundamental questions regarding massive star formation mechanisms. In particular, the apparently isolated and relatively slow moving single VMS VFTS 682 has raised suspicion that it might be an instance of isolated massive star formation (Bestenlehner et al. 2011; see Section 1). Additionally, the presence of single VMSs in R136 with inferred initial masses up to ≈320 M_☉ has questioned the canonical 150 M_☉ upper limit of the IMF (Crowther et al. 2010).

The most important outcome of our computations is the confirmation that a "slow runaway," with a three-dimensional velocity similar to that of VFTS 682, is in fact the most probable type of ejected VMS from an R136-like cluster (cf., Figure 3; Section 3.1). In fact, all of our computed models yield one or more runaways with kinematic properties agreeing fairly with those of VFTS 682 (cf., Table 1). Given such a likeliness of a VFTS 682-type runaway from R136, this apparently isolated star clearly does not imply isolated massive star formation and it is very likely a former member of R136.

Furthermore, as explained in Section 3, massive, close binaries are necessary to dynamically eject VMSs from star clusters, which, in turn, are susceptible to merge due to their hardening and/or eccentricity pumping, by the frequent close encounters that they receive. As our computations show (see Section 3), such massive binary mergers can easily produce single stars, within a few Myr, with masses well exceeding the 150 M_☉ upper limit, even if the cluster begins with the canonical upper limit. Our four computations have produced up to ≈250 M_☉ members and merger products up to ≈300 M_☉ are possible if the most massive binaries merge. Therefore, it does not seem to be a surprise that R136 has up to ≈320 M_☉ single-star members, given the large uncertainties in the stellar evolution models used to infer the masses (see Crowther et al. 2010, and references therein) and the presence of these super-canonical VMSs (or the presence of a super-saturated MF) is not an indication that R136 was born with an IMF having an upper limit that substantially exceeds the widely accepted 150 M_☉ limit. These VMSs can as well be massive binary mergers.

It is, of course, important to note that there are aspects in our assumptions and initial conditions that favor ejection of VMSs and formation of massive merger products. First, the assumption of no mass-loss during MS–MS mergers (see Section 2.2) is an idealization. Although SPH computations indicate that no mass is lost during an MS–MS collision (Sills et al. 2001), the studied collisions are typically between low-mass stars and such mass conservation is not necessarily true for massive MS–MS mergers. Furthermore, the merged MS star can spin rapidly, even beyond breakup, in which case the net mass gain would be much smaller or nil. Second, the assumption of a complete mixing (see Section 2.2) is also an idealization; it is possible that the denser helium-rich core(s) preferably sink to the center of the merger product, resulting only a partial mixing of unburned hydrogen with the helium core. This would lead to a partial refueling or rejuvenation of the merged MS star's core and hence a smaller lifetime. In other words, the treatment of complete mixing maximizes the lifetime of the merged star and hence its probability of getting dynamically ejected and/or being observationally detected. Although idealized, these conditions are those which can currently be adapted at best in a direct N-body calculation; other widely used direct N-body codes such as the "NBODY6++" (Spurzem 1999) and the "STARLAB" (Portegies Zwart et al. 2001) also utilize similar synthetic stellar evolution and merger recipes.

Third, the VMS ejection and the massive merger formation are also facilitated by our adapted initial or primordial mass segregation (see Section 2.1). While the motivation for choosing this condition is indeed to maximize these effects, primordial mass segregation is inferred to be true for several Galactic globular clusters (Baumgardt et al. 2008; Marks & Kroupa 2010; Strader et al. 2011) and open clusters (e.g., Bonnell & Davies 1998; Littlefair et al. 2003; Chen et al. 2007). An important outlook would be to study the effects of a varying degree of initial mass segregation.

Another ingredient that causes the ejection of runaway VMSs and the presence of m_s > 150 M_☉ VMSs in our computations likely is the presence of relatively close primordial binaries beyond m_s > 20 M_☉ (see Section 2.1). However, such a binary distribution is in accordance to the observational fact that O-stars are generally found in close binaries. Notably, the range of the O-star binaries' initial orbital periods in our models is taken to be similar to the observed one (Sana & Evans 2011). In this context, recall that we adapt a P-distribution that follows a single Öpik's law and also take the binary mass-ratios (q) close to unity (see Section 2.1). However, the P-distribution of O-stars in O-rich clusters is observed to follow a bi-uniform distribution in log₁₀P with the break at P ≈ 10 days; there is a significant overabundance of binaries (comprising 50%–60% of the O-star binary population) for P ≲ 10 days (see Equation 5.2 of Sana & Evans 2011). Furthermore, their mass ratios are found to be uniformly distributed between 0.2 ⩽ q ⩽ 1.0 (Sana & Evans 2011). A very important development would be to incorporate such more realistic distributions of O-star binaries' orbital periods and mass ratios. In fact, the Sana & Evans (2011) P-distribution would even more strongly favor the ejection of VMSs and the formation of super-saturated stars than the presently adapted distribution due to the overabundance of P ≲ 10 days binaries in the former distribution. This would, of course, be at least partially counterbalanced by the q-distribution; the smallest q gives a binary with a binding energy smaller by a factor of a few for the same primary mass and semimajor axis. It is, therefore, an essential outlook to study the net outcome of such effects. However, note that the gross conclusions, as obtained from the present models, are unlikely to get affected by these details since the range of P and hence the range of binary binding energies are indeed taken to be similar.

Finally, note that the chosen masses of the initial Plummer clusters are the upper mass limit of that of R136 (≈10⁵ M_☉; Crowther et al. 2010). A lower mass limit of R136 is ≈5 × 10⁴ M_☉ (Weidner & Kroupa 2004). This being only by a factor of two smaller than the upper limit, we do not expect that the results presented here would be affected significantly had we chosen the lower mass limit as the total initial cluster mass.

A corollary of our above two results is that runaway VMSs and super-canonical VMSs (with masses exceeding the canonical upper limit) would coexist, in general. Of course, as explained in Section 3, the latter type of VMSs can easily be ejected from their host clusters and may be found only as runaways. It would be worthwhile to search for other massive (and young enough) star clusters, in addition to R136, where such super-canonical VMSs exist and the clusters can be related to VMS runaways at the same time. Incidentally, the VMS ζ Pup (≈70 M_☉), whose kinematics indicate that it is a runaway from the Vela R2 association, has recently been inferred to be a merger product of at least two massive stars (Vanbeveren 2011).

The main drawback of this study is the incompleteness of the primordial binary population. As already stressed earlier, the absence of initial binaries for m_s < 5 M_☉ is anyway unlikely to significantly influence the ejected stellar population for M ≳ 5 M_☉, which justifies such a truncation for our purposes. The truncated primordial binary population, of course, has effects on less-massive ejecta which are pointed out in Section 3. The purpose of eliminating the lower mass binaries is, of course, computational ease. The presence of binaries substantially diminishes the speed of direct N-body computations (even on GPU platforms; see Section 2.2) and such computations with a full spectrum of primordial binaries (i.e., 100% primordial binary fraction) of models of the mass that we use in this study (i.e., N ≈ 170, 000; see Section 2.2) are currently prohibitive even with the fastest available GPUs and workstation processors. Nevertheless, we wish to emphasize again that our computations in this study are the largest-sized direct N-body computations of initially fully mass-segregated clusters in which all massive stars are in binaries reported so far, thanks to the remarkable boost of computing speed provided by the recent GPUs (see the NVIDIA GEFORCE Web site¹) and to the efficient numerical algorithms of NBODY6 (Aarseth 2003). These calculations represent the first-ever attempt to directly evolve a real-sized massive star cluster.

An immediate improvement over our work is, of course, to explore a full primordial binary spectrum and the effect of their varied distribution functions which, we expect, will become plausible in the near future given the continuing remarkable algorithmic improvements of NBODY6 and the availability of increasingly faster GPUs. A more accurate structural modeling of R136 is also pending once its dimensions are more clearly settled (see Section 2.1).

In conclusion, our calculations imply that the observed VMSs outside R136 and also its bound VMS members are fully consistent with the standard paradigm of massive star formation exclusively in dense environments and following the canonical IMF. Hence, the observations discussed here cannot be considered as instances violating these regimes.

We are thankful to Sverre Aarseth of the Institute of Astronomy, Cambridge, U.K., for his effort in the continuing remarkable improvements of NBODY6 and his gracious help during our computations. We are also thankful to Keigo Nitadori of RIKEN, Japan, for his developments that resulted in the exemplary hardware-acceleration of NBODY6. We thank the anonymous referee whose suggestions have improved the quality of the manuscript substantially.