The emergence of super-canonical stars in R136-type starburst clusters

Abstract

Among the most remarkable features of the stellar population of R136, the central, young, massive star cluster in the 30 Doradus complex of the Large Magellanic Cloud, are the single stars whose masses substantially exceed the canonical stellar upper mass limit of 150 M_⊙. A recent study by Banerjee, Kroupa & Oh (hereafter Paper I), indicates that such ‘super-canonical’ (SC) stars can be formed out of a dense stellar population with a canonical initial mass function (IMF) through dynamically induced mergers of the most massive binaries. The above study consists of realistic N-body computations of fully mass-segregated star clusters mimicking R136 in which all the massive stars are in primordial binaries. In the present work, we study the formation of SC stars in the computed R136 models of Paper I in detail. Taking into consideration that extraneous SC stars form in the computed models of Paper I due to the primordial binaries’ initial eccentricities, we compute additional models where all primordial binaries are initially circular. We also take into account the evolution of the mass of the SC stars and the resulting lifetime in their SC phase using detailed stellar evolutionary models over the SC mass range that incorporate updated treatments of the stellar winds. In all these computations, we find that SC stars begin to form via dynamical mergers of massive binaries from ≈1 Myr cluster age. We obtain SC stars with initial masses up to ≈250 M_⊙ from these computations. Multiple SC stars are found to remain bound to the cluster simultaneously within an SC lifetime. However, we also note that SC stars can be formed at runaway velocities which escape the cluster at birth. These properties of the dynamically formed SC stars, as obtained from our computations, are consistent with the observed SC stellar population in R136. In fact, the evolutionary models of SC stars imply that had they formed primordially along with the rest of the R136 cluster, i.e. violating the canonical upper limit, they would have evolved below the canonical 150 M_⊙ limit by ≈3 Myr, the likely age of R136. Thus, according to the new stellar evolutionary models, primordially formed SC stars should not be observable at the present time in R136. This strongly supports the dynamical formation scenario of the observed SC stars in R136.

1 INTRODUCTION

The study of the functional form of the number distribution of stars in galaxies, with which they are born, or the stellar initial mass function (IMF) has always been of fundamental importance (Bastian, Covey & Meyer ; Kroupa et al. ). The high-mass end of the stellar IMF is in particular focus due to the feedback it gives to the star-forming gas in the form of radiation pressure, kinetic energy from the stellar wind and supernova ejecta as well as the chemical enrichment from the latter source. An upper limit of ≈60 M_⊙ to the mass of a star could be set by the ‘Eddington limit’ (Eddington ) or the point of balance of gravity by the radiation pressure of a star. However, in practice stellar masses can easily exceed the Eddington limit since massive stars are not fully radiative but contain convective cores (Kippenhan & Weigert ). The next upper limit is set by the possibility of the destruction of a star due to thermal pulsations (Schwarzschild & Härm ; Beech & Mitalas ). Stothers () determined this limit to be at m_max* ≈ 120–150 M_⊙ for [Fe/H] ≈ 0 and m_max* ≈ 90 M_⊙ for [Fe/H] ≈ −1.

Difficulties also arise if one considers the growth of a protostar via gas accretion which is crucial for massive star formation. Here, one again encounters the ≈60 M_⊙ Eddington limit. Stellar formation models lead to a mass limit near 40–100 M_⊙ imposed by feedback on a spherical accretion envelope (Kahn ; Wolfire & Cassinelli ). Some observations suggest that protostars may be accreting material in discs rather than spheres (e.g. Chini et al. ), in which case it may be possible to overcome the radiation pressure at the equator of the protostar. Studies of the formation of massive stars through disc accretion with high accretion rates, thereby allowing the radiation to escape preferentially along the poles (e.g. Jijina & Adams ), indeed allow formation of stars with larger masses.

The feedback-induced mass limit can also be avoided if massive stars can form through mergers (Bonnell, Bate & Zinnecker ; Zinnecker & Yorke ). In this scenario, massive stars form via coalescence of intermediate-mass protostars in the cores of dense stellar clusters that have undergone core contraction due to rapid accretion of gas with low specific angular momentum. The high central density required for the mergers (10⁸ M_⊙ pc⁻³) is still difficult to achieve but it should be noted that an observable young cluster is necessarily disposed of a substantial fraction of its natal cloud and is thus likely to be always observed in an expanding and hence diluted phase. Recently set limits on the radii of embedded clusters indeed suggest them to form very compact (Marks & Kroupa ). In this context, it is worthwhile to note that in the computations discussed in the following sections we do get binary coalescence events even in moderately dense stellar clusters. This happens due to the presence of a substantial number of binaries in the cluster core which have much larger encounter cross-sections than single stars due to their bigger geometrical sizes.

While the existence of an upper limit of the stellar mass is rather obvious from theory since several decades, such an upper limit has been established from observations only recently. There have been some earlier indications of the presence of a cut-off near m_max* ≈ 150 M_⊙ from observations of young massive clusters, particularly in R136 (Massey & Hunter ; Massey ) in the Large Magellanic Cloud (LMC) and in Arches cluster (Figer ) close to our Galactic Centre, but these results were not sufficiently convincing in the sense that no statistical significance was attached to them. The observed upper limit was considered to be a limitation due to sampling rather than a true limit (Massey ). Elmegreen () also noted that random sampling from an unlimited IMF for all star-forming regions in the Milky Way would lead to the prediction of stars with masses ≳1000 M_⊙, unless there is a sharp downturn in the IMF beyond several 100 M_⊙.

Weidner & Kroupa () for the first time gave a critical look at the question of the existence of a physical upper limit of ≈150 M_⊙ in the IMF, as observed in R136 (the central massive star cluster in the 30 Doradus region of the LMC). Assuming R136 has a mass in the range 5 × 10⁴ < M_cl < 2.5 × 10⁵ M_⊙ (Selman et al. ), they found that for a canonical IMF with m_max* = ∞, 10 < N(>150 M_⊙) < 40 stars are missing in R136 that have masses >150 M_⊙. The probability that no stars are observed among the 10 expected ones, assuming m_max* = ∞, is p = 4.5 × 10⁻⁵, i.e. the observed massive stellar content of R136 implies a physical stellar mass limit at m_max* ≈ 150 M_⊙. Similarly, Figer () found a dearth of N(>150 M_⊙) = 33 stars in the Arches cluster, where the shallowing of its stellar mass function due to its rapid tidal dissolution has been incorporated. The corresponding number of missing stars for a canonical IMF is N(>150 M_⊙) = 18, which gives a non-detection probability as small as p = 10⁻⁸. Thus, the Arches cluster also has m_max* ≈ 150 M_⊙ with a very high significance. Finally, Oey & Clarke () studied nine clusters and associations in the Milky Way, LMC and Small Magellanic Cloud to investigate the expected masses of the most massive stars in these for different upper mass limits (120, 150, 200, 1000 and 10 000 M_⊙). They concluded that the observed number of massive stars supports the existence of a general physical upper mass cut-off within the range 120 < m_max* < 200 M_⊙ with a high significance.

The massive stellar population of young massive star clusters therefore indicates a physical upper mass limit near m_max* ≈ 150 M_⊙ which is the canonical upper limit of the stellar IMF. One can regard this mass as the limit imposed by the process of star formation under physical conditions achievable in a star-forming core. The origin of this very limit is under investigation and the value should currently be taken as empirical only. As coined by Kroupa et al. (), a stellar population containing stars with their zero-age main-sequence (ZAMS) masses up to 150 M_⊙ is ‘saturated’.

In this paper, our concern is again related to the massive stellar population in R136. The particular aspect that we focus on is related to a recent study by Crowther et al. (). These authors reanalysed the massive stellar population of R136 in unprecedented detail using Hubble Space Telescope and Very Large Telescope spectroscopy and high spatial resolution near-IR photometry to find four stars, within the central 1 × 1 pc of R136, with masses 165–320 M_⊙, i.e. substantially above the canonical limit. One can call such a stellar population, containing single stars with masses substantially exceeding 150 M_⊙, as ‘super-saturated’ and the stars surpassing the canonical upper limit can be called ‘super-canonical’ (SC; Kroupa et al. ). Although the stellar population of R136 has been studied by earlier authors, the crucial turn in Crowther et al. () is due to these authors’ consideration that the observed SC stars in R136, in spite of exhibiting WN-type spectra and possessing strong winds, are actually young, main-sequence (MS) stars rather than classical Wolf–Rayet (WR) stars. To take this into account, Crowther et al. () incorporate detailed stellar evolutionary models of very massive MS stars that include state-of-the-art treatment of stellar winds which lead to the inference of the SC masses.

More recently, however, Banerjee, Kroupa & Oh (, hereafter Paper I) showed that single stars with masses m_s > 150 M_⊙ can indeed form in R136 through mergers of the members in massive binaries. These authors computed the evolution of model clusters that resemble R136 using the state-of-the-art direct N-body integration code ‘nbody6’ (Aarseth ). In these computationally challenging models, all the O-stars were taken to be in tight binaries as observations indicate (Sana & Evans ; see Section ) but the ZAMS mass of a single component never exceeds m_max* = 150 M_⊙. They found that, in course of their dynamical evolution, the massive binaries in the model clusters can often merge due to hardening and/or eccentricity enhancement due to dynamical encounters to produce single stars of m_s > 150 M_⊙. This implies that the presence of SC stars in R136 does not illustrate a violation of the canonical upper limit.

Of course, the value of the upper stellar mass limit is not necessarily exactly = 150 M_⊙ (Oey & Clarke ). However, for definiteness, we take the following hypothesis: there exists a fundamental upper mass limit of m_max* = 150 M_⊙ or the canonical limit of ZAMS stars formed in clusters and any observed more massive (‘SC’) star is created from dynamically induced stellar mergers within the dense young cluster. In this study, we aim to test the above hypothesis or answer the following question: can the observed number of SC stars in R136 be explained by the naturally occurring stellar dynamical and stellar evolution processes in R136? To that end, we utilize the four direct N-body computed massive cluster models of Paper I mimicking R136 (Sections and ) and additionally perform five more similar computations (Section ) to trace the formation of SC stars through dynamical means. These computed models initiate with properties that are consistent with the observed properties of young clusters, in particular primordial mass segregation (Littlefair et al. ; Chen, de Grijs & Zhao ) and massive stellar binaries with component mass ratio close to unity (Sana & Evans ). Furthermore, we utilize detailed rotating stellar evolution models in the SC mass range to estimate the mass evolution of the dynamically formed SC stars and the resulting lifetimes in their SC phases (Section ). We conclude the paper by summarizing our results and pointing out the limitations of the present study (Section ).

2 COMPUTATIONS

State-of-the-art calculations have been performed in Paper I to evolve model star clusters that mimic R136 in terms of observable parameters. The primary objective of these calculations was to study the ejected OB-stars from R136, in particular, whether massive runaway stars like VFTS 682 and 30 Dor 016 can indeed be ejected dynamically from R136 as suspected (Evans et al. ; Bestenlehner et al. ). Besides good agreement of the kinematics of the massive ejecta with those of these noted runaways, it is also noted in Paper I that the tight massive binaries which drive these runaways can merge due to the dynamical interactions leading to the formation of SC stars in each of the models, in agreement with the observed stellar population in R136 (Crowther et al. ). In this study, we continue to utilize the results of these computations. The initial conditions and the method of these calculations are accounted in detail in section of Paper I which we briefly recapitulate here.

2.1 Initial conditions

The above computations comprise direct N-body integrations of Plummer spheres (Kroupa ) with parameters conforming with those of R136. The initial mass of the Plummer spheres is M_cl(0) ≈ 10⁵ M_⊙, which is an upper limit for R136 (Crowther et al. ), and the half-mass radii are taken to be r_h(0) ≈ 0.8 pc. The clusters are made of stars with ZAMS masses drawn from a canonical IMF (Kroupa ) over the range 0.08 < m_s < 150 M_⊙ and are of metallicity Z = 0.5 Z_⊙.

As for the primordial binary population, stars with m_s > 5 M_⊙ are all in binaries, while all lighter stars are kept single. This termination of the binary population is due to computational ease; binaries bottleneck the calculation speed of direct N-body integration significantly so that adopting a full spectrum of primordial binaries (i.e. 100 per cent initial binary fraction) in models of the size that we compute becomes prohibitive. Disregarding the binary population for m_s < 5 M_⊙, of course, does not affect the ejection of massive stars and mergers of massive binaries significantly (see Paper I for details). Following the observed period distribution of O-star binaries (Sana & Evans ), the orbital periods of the binaries having primary masses m_s > 20 M_⊙ are chosen from a uniform distribution in the range 0.5 < log₁₀P < 4, where P is the orbital period in days. For primary masses 5 < m_s < 20 M_⊙, the ‘birth period distribution’ of Kroupa () is adopted. The eccentricity distribution is chosen to be thermal (Spitzer ; Kroupa ). These distributions, however, are not modified via ‘eigenevolution’ (Kroupa ) as the binary systems’ internal interactions between massive pre-MS stars are not yet quantified.

The initial configurations are fully mass segregated (Baumgardt, De Marchi & Kroupa ) so that the massive binaries concentrate within the clusters’ central region. In other words, it is assumed that the clusters are primordially segregated, which was the case for several Galactic globular clusters (Marks & Kroupa ) and young star clusters (Littlefair et al. ; Chen et al. ). In these mass-segregated initial models, the stars are redistributed over the cluster in such a way that the individual mass groups are in equipartition, whilst the overall cluster is in virial equilibrium as well. The algorithm for creating such systems is described in the appendix of Baumgardt et al. (). We emphasize that we are conforming with the observed states of massive, young clusters by starting our models binary rich and mass segregated, i.e. we are not choosing some particular fine-tuned initial conditions in order to enhance our results.

2.2 N-body integration, stellar evolution and stellar mergers

With such an initial set-up, these computations are the very first direct N-body calculations under such extreme conditions (all massive stars in binaries, fully mass segregated). To perform the direct N-body computations, the state-of-the-art ‘nbody6’ integrator (Aarseth ) has been used. In addition to integrating the particle orbits using the highly accurate fourth-order Hermite scheme and dealing with the diverging gravitational forces in close encounters through regularizations, nbody6 also employs the well-tested analytical stellar and binary evolution recipes of Hurley, Pols & Tout () and Hurley, Tout & Pols () (the SSE and the BSE schemes). These prescriptions are based on model stellar evolution tracks computed by Pols et al. (). The wind mass loss of the massive MS stars is adopted using the empirical Nieuwenhuijzen & de Jager () formula. The mass-loss rates are, of course, appropriately modified on the giant branches and other evolved phases for stars of all masses (Hurley et al. ). Admittedly, the Nieuwenhuijzen & de Jager () mass-loss rate is too simple for massive stars; it is likely to grossly underestimate the wind mass loss of massive stellar merger products, which we focus on here, as these objects usually become WR stars even in their MS phase (Glebbeek et al. ). The above stellar wind recipe is currently the best available one in a direct N-body evolution code; other widely used direct N-body codes, e.g. the ‘starlab’ (Portegies Zwart et al. ) and the ‘nbody6++’ (Spurzem ), utilize analogous wind treatments. Since the wind mass loss is crucial for very massive stellar entities (Glebbeek et al. ), we consider the effects of more detailed stellar winds separately in Section .

Mergers among MS stars are treated in nbody6 by adopting the schemes of Hurley et al. (, ). When two MS stars collide (a collision between two single MS stars or between two MS components in a binary due to eccentricity induced by close encounters and/or due to encounter hardening), the two main assumptions in the above scheme are (i) the merged object is an MS star with the stellar material completely mixed and (ii) no mass is lost from the system during the hydrodynamical process leading to the merger. This no mass-loss assumption is based on the results of smoothed particle hydrodynamics simulations of MS–MS mergers (e.g. Sills et al. , who show that the mass loss is up to 5 per cent for low-mass MS stars) but such studies also show that the mixing is only partial. The age of the merged MS star is assigned based on the amount of unburnt hydrogen fuel gained by the hydrogen-burning core as a result of the mixing. In the case of a mass transfer across an MS–MS binary, the cases of the original accretor MS star having a radiative and a convective core are treated separately (included in the BSE algorithm). For a convective core, the core grows with the gain of mass and mixes with the unburnt hydrogen fuel so that the accreting MS star appears younger. For the case of a radiative core, the fraction of the hydrogen burnt in the hydrogen-burning core remains nearly unaffected by the gain of mass so that the effective age of the MS star decreases. However, in any case of a merger between two MS stars, an appropriate amount of mass is removed from the final product if the kinetic energy of their final approach exceeds the binding energy of the merged MS star. Admittedly, there exists some arbitrariness in applying the treatments, appropriate for mergers of low-mass stars, to the cases involving mergers of ≈100 M_⊙ stars which we consider in the following sections. The structure of these stars is very different so that it is unclear whether the above collisional mass-loss and mixing criteria will remain valid for such massive stellar mergers. Unfortunately, the outcomes of massive stellar collisions are currently unclear from theoretical studies. A relatively simple but very fast algorithm for treating collisions of massive stars (the MMAMS scheme) by Gaburov, Lombardi & Portegies Zwart () shows <10 per cent mass loss and partial mixing in encounters of ∼10 M_⊙ stars. In the absence of a clear understanding, we continue to utilize the scheme for low-mass stellar mergers for massive mergers also.

Note that the computed models 1–4 of Paper I constitute exactly the same cluster model of four different randomized (discrete) realizations. Each of these models are evolved up to ≈3 Myr, which is the widely used age of R136 as inferred from its low-mass stellar population (Andersen et al. ). The above age estimate is perhaps the most robust one for the R136's stellar population.

2.3 Additional model computations

In the four model calculations presented in Paper I (models 1–4), there are typically one to two massive binary mergers per cluster at the very beginning of the computations that form SC stars. These mergers are found to happen mostly due to high initial eccentricities of the corresponding massive primordial binaries. Although these merger products are highly SC, we consider such immediate mergers as ‘accidental’ since the parent binaries would have already merged during their pre-MS evolution. However, we adhere to the canonical upper limit of m_max* = 150 M_⊙ for R136, which is a plausible assumption as the limit is deduced from observations of young massive star clusters. Hence, we consider the SC stars that appear due to a collision between the two members of a primordial binary at the very beginning of the cluster evolution (cf. Fig. ⁰⁰⁰¹) as spurious and unrelated to the formation of the initial ZAMS stellar population. These are artefacts of the chosen eccentricity distribution. The SC stars that are formed later in the course of the cluster evolution are, of course, considered genuine.

Notably, the presence of one to two spurious (single) SC stars is not expected to affect the dynamical evolution of our model clusters in any significant way as, for a massive primordially mass-segregated system like ours, there are anyway a large number of massive stars in the cluster's central region. Therefore, the inclusion of a few more massive single stars is unlikely to provide a substantial additional effect. It is, of course, not impossible that two massive stars happen to collide to form a single SC star at the birth of a cluster, but it seems unphysical that it could be an eccentricity effect which is always the case for the above-mentioned accidental SC stars, which is why we consider them spurious. The conditions under which such ‘primordial mergers’ can happen and their corresponding probabilities are beyond the scope of the present study. In any case, apart from their negligible contribution to the dynamics of their host clusters, it is unlikely to find any such primordially merged star in its SC state in R136 at its current age, due to the short SC lifetime, so that the currently observed ones should all be dynamically formed. We elaborate the latter point in Section .

To see whether the chosen thermal eccentricity distribution in the computed models of Paper I does have a systematic effect on the dynamical formation of SC stars, we perform five additional computations with the same initial configuration (but different realizations) as in Paper I (see Section ) except that all the primordial binaries are initially taken to be circular. Of course, the most realistic approach would be to ‘eigenevolve’ (Kroupa ), i.e. to incorporate the pre-MS evolution of each of the primordial binaries to determine their eccentricities and mass ratios at the beginning of the computations. However, the nature of eigenevolution is not yet known for progenitors of massive stars. Typically, the tidal interaction between the proto-massive-star members would circularize the tight binary orbits but wide binaries would still remain eccentric. Currently, it is unclear in the literature what would be the exact outcome of eigenevolution of massive protostellar binaries. As we shall see in Section , if only the dynamically formed SC stars are considered, their population is similar for the calculations in Paper I and in these new computations, i.e. the primordial binaries’ eccentricity distribution does not play a significant role.

We discuss the formation of SC stars from the computations of Paper I and from these new computations in Section .

3 RESULTS: FORMATION OF SUPER-CANONICAL STARS

Fig. ⁰⁰⁰¹ (top panel) shows the masses of the SC stars, M_SC, with the cluster evolution time for the computed ‘model 3’ of Paper I. Each appearance of a line in this figure corresponds to the appearance of an SC star and its slow decline is due to the wind mass loss of the SC star as provided by the currently available stellar wind prescriptions in nbody6 (see Section ). As pointed out in Section , this stellar wind scheme is inappropriate for very massive stars such as the SC ones, as it substantially underestimates the mass loss and enhances the lifetime in the SC phase. We address this issue in Section .

The SC stars represented by the grey lines that appear from the very beginning of the evolution are formed due to immediate mergers of highly eccentric massive primordial binaries which are unintended and considered spurious as explained in Section . The rest of them are of course genuine SC stars that are formed through dynamically induced mergers of massive binaries (see Section ). This particular model formed a rather marked number of SC stars. Notably, an M_SC ≈ 240 M_⊙ SC star is formed for the first time as early as T₀ = 0.7 Myr.

Fig. ⁰⁰⁰¹ (bottom panel) similarly shows the formation of SC stars for ‘model 4’ of Paper I. Here, only one spurious SC star is formed which is only marginally more massive than the canonical upper limit and hence disappears soon due to wind mass loss. The subsequent cluster contains only dynamically formed SC stars. In this model also, the first appearance of an SC star (M_SC ≈ 170 M_⊙) occurs quite early (T₀ = 1.2 Myr).

Fig. ⁰⁰⁰² shows the formation of SC stars in the additional models ‘C2’, ‘C5’ and ‘C10’ with all primordial binaries initially circular as explained in Section . While the total number of SC stars formed in these models is typically similar to that in the models of Paper I (cf. Table ⁰⁰⁰¹), none of the stars is formed artificially as can be seen in Fig. ⁰⁰⁰². The SC stars start to appear from typically T₀ ≈ 1 Myr for these models also. The remaining two of the additionally computed models produce SC stars only with masses marginally above 150 M_⊙.

3.1 Dynamical formation of SC stars

In the above computations, we find several dynamical channels for the formation of the SC stars. The most frequent way is found to be an abrupt eccentricity enhancement of a massive binary due to a close encounter with a single star. A collision and subsequent merger occur between the components when their periastron separation becomes smaller than the sum of their radii. The single stars in the cluster's central region arise from dissociation of wider binaries by encountering with harder binaries or by direct ionization (i.e. detachment of a binary due to its encounter with a single star of kinetic energy larger than its binding energy). This channel is most fruitful for the models 1–4 as they generally contain more eccentric binaries.

The next important channel is the eccentricity augmentation of the inner massive binary due to Kozai cycles with an outer companion. Such hierarchical triples can often form in the cluster centre as a result of binary–binary encounters. Notably, a close encounter of a single star with a hard binary can also form a temporary triple system (resonance encounters; see Heggie & Hut ). The orbital changes due to tidal interactions within the inner massive binary and also that of the outer member with the inner members as well as the merger process of the inner binary can often make the triple unstable. This causes the inner binary to become isolated from the outer companion in the course of its eccentricity growth. In our computations, a few newly formed SC stars are still found bound in wide orbits with the outer member that soon gets ionized (within a few orbital times). In a few cases, an SC star is also found to form when the outer member of a triple collided with one of the inner members.

SC stars also form due to encounter hardening (Heggie ) of the orbits of appropriate binaries resulting in a Roche lobe overflow. The approach to the Roche lobe contact is generally assisted by the evolutionary expansion (during its MS) of the more massive binary member which therefore fills the Roche lobe first (cf. Hurley et al. ). The subsequent unstable mass transfer from the more massive member to the less massive one coalesces the binary.

Some of the SC stars are found at high velocities at their birth. These are the ones whose progenitor binaries suffered substantial recoils in the close encounters that caused them to merge and they escape from the cluster right at their birth. As examples, the radial distance, R_SC, versus t plots for the SC stars in Fig. ⁰⁰⁰³ show such runaways in the models ‘3’ and ‘C10’ exhibiting the monotonic outgoing trajectories. In fact, such born escapers constitute most of the runaway stars in the SC mass range. This can be expected since, once formed, it is unlikely to eject an SC star later by a super-elastic close encounter (Heggie ) with a hard primordial binary unlike the non-SC stars for which such an ejection channel is usual. The SC star, being the most massive participant in such an encounter (the primordial binary members are <150 M_⊙), is likely to get trapped in the binary by exchanging with one of the companions and hence receiving less recoil.

The rest of the SC stars initially form and remain bound to the cluster. Because they are among the most massive members of the cluster, they continue to reside close to the cluster's centre, occasionally several of them simultaneously (cf. Fig. ⁰⁰⁰³). However, a few of them can subsequently be ejected by dynamical recoils (only two such SC ejecta occurred from all of our computations). Notably, the ejection or runaway fraction of stars in the SC mass range is >0.5 as determined in Paper I (see fig. 4 of Paper I).

It is to be noted that the above dynamical channels are also applicable for formation of merger products with m_s < 150 M_⊙.

4 EVOLUTION OF SUPER-CANONICAL STARS: THE EFFECT OF STELLAR WINDS

As already raised in Section , the treatment of stellar wind mass loss in nbody6 is inappropriate for massive (MS) entities like the SC stars which grossly underestimates their wind mass loss, thereby enhancing their lives in the m_s > 150 M_⊙ state. In Fig. ⁰⁰⁰¹, it can be seen that all the SC stars formed in ‘model 3’ remain SC until the end of the computation at t ≈ 3 Myr. This is also true for all the substantially SC merger products in the other computations.

To take into account an appropriate mass loss for the SC stars, we consider the stellar evolutionary models of rotating MS stars by Köhler & Langer (in preparation), which encompass the SC mass range and are computed for the LMC metallicity. The above authors use the one-dimensional hydrodynamic binary evolution code of Heger, Langer & Woosley () which includes direct treatments of stellar rotation, magnetic field and detailed stellar winds. The code has been recently improved as in Petrovic et al. () and Yoon, Langer & Norman ().

Stellar wind is the most important factor for the evolution of the mass of an isolated massive star during its MS. In the computations by Köhler & Langer (in preparation), the MS stellar wind of Vink et al. () is applied that includes the bi-stability jump. During the bi-stability jump, a linear interpolation is applied as in Brott et al. (). The mass-loss rate is switched from Vink et al. () to that of Nieuwenhuijzen & de Jager () when the latter becomes stronger (at effective temperatures <22 000 K). When the surface helium abundance becomes Y_s ≥ 0.7, the star is considered to enter the WR phase in which case the much enhanced Hamann, Koesterke & Wessolowski () wind, divided by a factor of 10 (Yoon et al. ), is applied.

Fig. ⁰⁰⁰⁴ shows the resulting mass evolution of SC stars with initial surface rotation velocities of V_rot(0) = 50 and 500 km s⁻¹. The SC stars undergo a steep mass loss as they become WR stars. Because of the large mass loss in the WR phase, these stars become less massive than 150 M_⊙, i.e. cease to be SC, in τ_SC ≈ 1.4–1.8 Myr from their zero age. The lifetime in the SC phase therefore becomes much shorter with the more detailed mass-loss prescription in a hydrodynamical stellar evolution calculation. We consider this more realistic lifetime as the lifetime of an SC star rather than that obtained from within nbody6. We note that the SC stars formed in our computations do not undergo further coalescence or Roche lobe overflow with a secondary (although these are in principle possible) so that the individually obtained mass-loss rates can be applied to them (see also Section ). It can be noted from Fig. ⁰⁰⁰⁴ that the mass-loss rate is only moderately affected by the stellar spin. The SC stars being binary merger products are indeed likely to form with high surface rotation velocities.

Notably, in a dynamically active environment, an SC star can undergo further mass evolution due to subsequent merger of the SC product with other stars or its acquirement of a binary companion causing Roche lobe overflow. However, in our computations, we hardly find subsequent mergers of the SC stars or them being overflown as confirmed by their continuous and smooth mass depletion (cf. Figs ⁰⁰⁰¹ and ⁰⁰⁰²), justifying the application of isolated stellar mass loss to them.

4.1 Can the dynamically formed SC stars in R136 be observed?

Given that the lifetime of the SC phase (for the LMC metallicity) is τ_SC ≈ 1.4–1.8 Myr, how likely is it to detect the dynamically formed SC stars in R136 at the present day? Clearly, the SC stars must form at most τ_SC time earlier to be presently visible in their SC phases. From Figs ^0001- 0003 and Table ⁰⁰⁰¹, it is clear that (a) SC stars are very likely to form at any time within 1–3 Myr cluster age and (b) multiple SC stars are likely to appear and remain bound to the cluster within an SC lifetime of τ_SC = 1.5 Myr, the latter value being taken representatively.

These conclusions imply that in spite of the severe wind mass loss of the SC stars, it is still feasible that they can be detected at the present day in multiple numbers near the centre of an R136-type cluster. Such SC stars are the ones that are formed via dynamical mergers of massive binaries late enough that they remain SC until now. Our calculations indicate that the satisfaction of this condition is quite feasible for an R136-type cluster. For example, at the cluster age of 2.5 Myr, ‘model 3’ would retain three SC stars (cf. Figs ⁰⁰⁰¹ and ⁰⁰⁰³) and C5 and C10 would retain three and two SC stars, respectively (cf. Figs ⁰⁰⁰² and ⁰⁰⁰³).

By the same argument, the above detailed evolutionary models of the SC stars imply that those observed in R136 would have evolved below the canonical limit by the present day had they been primordial, i.e. had they been formed at the cluster's birth, thereby violating the canonical upper limit, assuming the widely used ≈3 Myr age of R136 (Andersen et al. ). This strongly supports our dynamical formation scenario of the observed SC stars in R136 which causes them to form sufficiently later so as to make them visible during their SC phases at the present day. This is well supported by our computations.

Admittedly, there are several drawbacks and incompleteness in the above line of arguments. First, the detailed stellar evolutionary models of Köhler & Langer (in preparation) have been considered independently of the dynamics of the model clusters. The newer stellar evolutionary models would result in stronger wind mass loss from the cluster core and would contribute to the dilution of the cluster's central potential to a larger extent than that from the nbody6's native stellar evolutionary scheme (Hurley et al. , ). However, in our computations the cluster's core, where the binary mergers occur, is populated by a large number of hard massive binaries due to our initial conditions (Section ), which are all in accordance with observations of young star clusters (see Section and references therein). In that case, the dynamics of the cluster's central region will be dominated by the energy generated in the super-elastic binary–single and binary–binary encounters (Heggie ; Heggie & Hut ) and the stellar mass loss would not affect the dynamical encounter rates in the cluster's centre significantly.

It is true that the mass loss from the individual companions of a binary would result in expansion of the binary orbit and hence dilute its hardness. However, as shown by Fujii & Portegies Zwart (), this effect is counteracted by the dynamical hardening (Heggie ) of the binaries so that they remain hard throughout. In the N-body calculations of Fujii & Portegies Zwart (), even just a single binary survives to remain hard in their model cluster's core and dominates the energy generation in the cluster's core (until its ejection). With a large number of binaries, as in our case, the encounter hardening and the corresponding energy generation would be substantially stronger as the binaries provide much larger encounter cross-sections than single stars due to their much bigger geometrical size (see also Fujii, Saitoh & Portegies Zwart ). Notably, even after including stellar wind mass loss, the computed single-star-only models of Fujii & Portegies Zwart () undergo deep core collapse within ≈3 Myr, i.e. the effects of dynamical relaxation and negative specific heat (Spitzer ) themselves dominate over the wind mass loss. Our models, of course, do not core collapse due to the substantial energy release from the central primordial binaries.

Another drawback of the above arguments is that they are based on stellar evolution models that begin from ZAMS with normal composition. The SC stars are merger products of O-stars that are evolved from their zero age and therefore would be He enriched and substantially rotating. While we do consider rotating models, the initial He fraction is taken to be that appropriate for the LMC, viz. Y = 0.2562 (Köhler & Langer, in preparation). To obtain a basic estimate of the He richness of the merged product, let us consider two 100 M_⊙ stars that merge to form a typical 200 M_⊙ SC star at a typical cluster age of 1.5 Myr. Each of the initial stars (before merger) would possess an H-burning core of M_c ≈ 80 M_⊙ and at the above merging age their core He fraction Y_c ≈ 0.5 (i.e. ≈40 M_⊙ He per star) according to the above newer models (Köhler, private communication). Assuming an ideal complete mixing and no mass loss during the merger process, the new He abundance of the merged product would be Y ≈ 0.4 which would also be the rejuvenated surface value Y_s; for a partial mixing the surface abundance would be lower. While this Y_s is large enough to drive a stronger wind from the newly formed SC star compared to its ZAMS value Y_s ≈ 0.25, it is still moderately above the ZAMS value and less than the abundance at which the WR winds turns on (Y_s ≥ 0.7; see Section ). Therefore, a drastically large wind from the newly formed SC star is not expected.

The above discussions imply that none of the shortcomings or uncertainties of our approach dictates the impossibility of observations of SC stars in R136; rather such observations are still feasible.

5 DISCUSSION AND SUMMARY

Table ⁰⁰⁰¹ summarizes the SC star formation in the computed models of Paper I and in the new computations. Our general conclusions from the computations are as follows.

• Formation of SC stars due to dynamically induced mergers of massive binaries is common in an R136-type cluster. In most of the computations, we find multiple SC stars formed within <3 Myr which is appropriate for R136 (Andersen et al. ; see Figs ⁰⁰⁰¹ and ⁰⁰⁰²). Most SC stars are formed single, while a few of them initially form in wide binaries.

• The SC stars typically begin to appear from T₀ ≈ 1 Myr cluster age or even earlier (cf. Table ⁰⁰⁰¹). They tend to form with equal likeliness over a cluster age of 1–3 Myr.

• The most massive SC star formed in a given computed model is typically close to M_max ≈ 200 M_⊙ and the most massive one formed is ≈250 M_⊙ (‘model 3’; cf. Table ⁰⁰⁰¹). The cluster age T_max corresponding to the formation of the most massive SC star is typically well within <3 Myr (cf. Table ⁰⁰⁰¹).

• Multiple SC stars are occasionally found to exist simultaneously near the cluster's centre which are bound to the cluster, over a representative SC phase lifetime of τ_SC = 1.5 Myr, within <3 Myr (cf. Table ⁰⁰⁰¹) cluster age. This, along with the second point above, implies that it is quite plausible that R136 harbours multiple SC stars at the present day.

• Some SC stars are formed with runaway velocities and escape (cf. Section ; Fig. ⁰⁰⁰³).

These conclusions conform with the observations by Crowther et al. (), who found four SC stars in R136 in the initial mass range of 165–320 M_⊙. In our computations, SC stars of up to ≈250 M_⊙ are formed, which is consistent with the observations taking into account the large uncertainties in the stellar evolution models in this mass range. Also, only one of the models (‘model 3’) has up to four SC stars simultaneously bound to the cluster (over a 1.5 Myr period) although other models also host multiple SC stars close to the cluster's centre. Some models, however, contain only one bound SC star within 3 Myr (cf. Table ⁰⁰⁰¹).

At this point, it is worthwhile to note that our above conclusions depend somewhat on the age of the R136 cluster which is, of course, not fully settled yet. While the bulk of R136 is ≈3 Myr old (Andersen et al. ), the high-mass stellar population might be younger (de Koter, Heap & Hubeny ; Massey & Hunter ). In our computations, the SC stars typically begin to appear from ≈1 Myr cluster age and several of them appear by ≈2.5 Myr. This time frame is quite consistent with the possible ≲2 Myr age of R136's massive stellar population given the wide uncertainties in the age of such a young cluster particularly for the massive end of the IMF. There has so far been no direct estimate of the age of the massive end of the stellar population of R136 and the above age limit is based only on comparisons with the stellar populations of Car OB1 and NGC 3603 (Crowther et al. ). Given that the stellar IMF of R136 continues to maintain the canonical law (Kroupa ) from the low to the high mass range (1.1–120 M_⊙; Massey & Hunter ; Andersen et al. ), we find it more natural to consider that the whole stellar population of the R136 cluster has formed in a single starburst event without a significant age spread. The issue can, of course, be more resolved with better understanding of the evolution of very massive stars and their winds. Finally, a recent study by Chini et al. () shows that the O-star binary distribution can, in fact, be even harder (i.e. more bound or tighter) than that considered in our computations (see Section ). This would lead to the formation of SC stars even earlier, which constitutes an important future study.

It is important to remember that the formation of SC stars at early ages in the above computations is facilitated by the adopted initial complete mass segregation. This condition subjects the massive binaries to strong dynamical encounters from the beginning of the cluster evolution. Note, however, that primordial mass segregation is inferred to have been true for several Galactic globular clusters (Baumgardt, De Marchi & Kroupa ; Marks & Kroupa ; Strader, Cladwell & Seth ) and open clusters (Littlefair et al. ; Chen et al. ). Notably, for a completely unsegregated model, it would take ≈10 Myr for the massive binaries to segregate to the cluster's centre which makes the early formation of SC stars unlikely. However, the mass segregation time-scale shortens substantially with increasing compactness of the cluster. In this context, an important outlook would be to study the formation of SC stars with varying initial compactness and degree of primordial mass segregation.

The drawbacks of the computed R136 models, as discussed in Paper I (see section of Paper I), naturally carry over to the present analyses. These limitations, however, do not crucially affect SC formation. In particular, the truncation of the binary distribution at m_s ≈ 5 M_⊙ does not influence the mergers of the most massive binaries, the latter being much more centrally concentrated due to mass segregation. Also, the exact period distribution of the O-star binaries is not instrumental for the formation of the SC stars as long as O-stars are largely found in tight binaries (Sana & Evans ). Notably, our adopted range of the orbital period for the O-star binaries is similar to that reported by Sana & Evans () (see Paper I). The formation of SC stars is generally similar for models with initially thermally distributed eccentricities (Paper I) and for the newer ones with initially circular binaries (except for the spurious SC stars in the former models). This implies that the initial eccentricity distribution of the massive binaries does not crucially influence the formation of SC stars.

The most consistent and accurate way to address the above-discussed concerns regarding the effects of a stronger wind mass loss due to the newer stellar evolutionary models is to incorporate the latter in a direct N-body code, which is well beyond the scope of this paper. No such direct N-body code exists at the present time. Furthermore, there are substantial uncertainties in the physics of mergers of massive stars at the present time. Therefore, the conclusions in this work are the best one can draw given the current technical limitations. Perhaps, it would be possible to do such work in future with a distributively computing, highly modular N-body calculation framework such as ‘muse’ (Portegies Zwart et al. ). In spite of all these uncertainties, it is worthwhile to note that if one focuses to the particular case of R136, then (a) the existence of SC stars is strongly supported by observations and (b) considering the widely used ≈3 Myr age of R136, the SC stars must have formed later than the birth of R136. Our computations are consistent with the formation of SC stars through dynamically induced mergers of massive binaries. While there are technologically limited drawbacks in our present analysis, as elaborated in Section , none of these shortcomings dictates the non-observation of SC stars in R136. Therefore, we can say that we have justified the hypothesis of Section through our detailed N-body modelling of R136 with initial conditions chosen in accordance with observations of young star clusters. In other words, from our realistically modelled computations of R136-like star clusters, and from our analyses as presented above, it can now be more definitely concluded that the observed super-saturated stellar population in R136 does not imply a violation of the canonical stellar upper mass limit near m_max* = 150 M_⊙.

We are very thankful to Karen Köhler of the Argelander-Institut für Astronomie, University of Bonn, Germany, for providing us with the massive stellar evolutionary models and discussing them in detail, which have been instrumental for this work. We thank Sverre Aarsteh of the Institute of Astronomy, Cambridge, UK, for his effort in the continuing improvements of nbody6 and his help during our computations. We wish to thank Paul Crowther of the University of Sheffield, UK, for encouraging discussions. We thank the anonymous referee for guiding criticisms that led to a substantial improvement of the manuscript.

REFERENCES