The R136 star cluster hosts several stars whose individual masses greatly exceed the accepted 150 M⊙ stellar mass limit (original) (raw)

Journal Article

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

,

Search for other works by this author on:

Search for other works by this author on:

• PDF
• Split View
• • Article contents
  • Figures & tables
  • Video
  • Audio
  • Supplementary Data
• Cite

Cite

Paul A. Crowther, Olivier Schnurr, Raphael Hirschi, Norhasliza Yusof, Richard J. Parker, Simon P. Goodwin, Hasan Abu Kassim, The R136 star cluster hosts several stars whose individual masses greatly exceed the accepted 150 M⊙ stellar mass limit, Monthly Notices of the Royal Astronomical Society, Volume 408, Issue 2, October 2010, Pages 731–751, https://doi.org/10.1111/j.1365-2966.2010.17167.x
Close

Navbar Search Filter Mobile Enter search term Search

Abstract

Spectroscopic analyses of hydrogen-rich WN 5–6 stars within the young star clusters NGC 3603 and R136 are presented, using archival Hubble Space Telescope and Very Large Telescope spectroscopy, and high spatial resolution near-IR photometry, including Multi-Conjugate Adaptive Optics Demonstrator (MAD) imaging of R136. We derive high stellar temperatures for the WN stars in NGC 3603 (T*∼ 42 ± 2 kK) and R136 (T*∼ 53 ± 3 kK) plus clumping-corrected mass-loss rates of 2–5 × 10−5 M⊙ yr−1 which closely agree with theoretical predictions from Vink et al. These stars make a disproportionate contribution to the global ionizing and mechanical wind power budget of their host clusters. Indeed, R136a1 alone supplies ∼7 per cent of the ionizing flux of the entire 30 Doradus region. Comparisons with stellar models calculated for the main-sequence evolution of 85–500 M⊙ accounting for rotation suggest ages of ∼1.5 Myr and initial masses in the range 105–170 M⊙ for three systems in NGC 3603, plus 165–320 M⊙ for four stars in R136. Our high stellar masses are supported by consistent spectroscopic and dynamical mass determinations for the components of NGC 3603A1. We consider the predicted X-ray luminosity of the R136 stars if they were close, colliding wind binaries. R136c is consistent with a colliding wind binary system. However, short period, colliding wind systems are excluded for R136a WN stars if mass ratios are of order unity. Widely separated systems would have been expected to harden owing to early dynamical encounters with other massive stars within such a high-density environment. From simulated star clusters, whose constituents are randomly sampled from the Kroupa initial mass function, both NGC 3603 and R136 are consistent with an tentative upper mass limit of ∼300 M⊙. The Arches cluster is either too old to be used to diagnose the upper mass limit, exhibits a deficiency of very massive stars, or more likely stellar masses have been underestimated – initial masses for the most luminous stars in the Arches cluster approach 200 M⊙ according to contemporary stellar and photometric results. The potential for stars greatly exceeding 150 M⊙ within metal-poor galaxies suggests that such pair-instability supernovae could occur within the local universe, as has been claimed for SN 2007bi.

1 INTRODUCTION

Although the formation of very massive stars remains an unsolved problem of astrophysics (Zinnecker & Yorke 2007), the past decade has seen a shift from a belief that there is no observational upper stellar mass cut-off (e.g. Massey 2003) to the widespread acceptance of a limit close to 150 M⊙ (Figer 2005; Koen 2006). If stars above this limit were to exist, they would be exclusive to the youngest, highest mass star clusters, which are very compact (Figer 2005). However, spatially resolved imaging of such clusters is currently exclusive to the Milky Way and its satellite galaxies, where they are very rare. In addition, the accurate determination of stellar masses generally relies upon spectroscopic and evolutionary models, unless the star is a member of an eclipsing binary system (Moffat 2008), something which occurs extremely rarely.

Most, if not all, stars form in groups or clusters (Lada & Lada 2003). An average star forms with an initial mass of ∼0.5 M⊙ while the relative proportion of stars of higher and lower mass obeys an apparently universal initial mass function (IMF; Kroupa 2002). In addition, there appears to be a relationship between the mass of a cluster and its highest mass star (Weidner & Kroupa 2006; Weidner, Kroupa & Bonnell 2010). High-mass stars (>8 M⊙), ultimately leading to core-collapse supernovae (SNe), usually form in clusters exceeding 100 M⊙ while stars approaching 150 M⊙ have been detected in two 104 M⊙ Milky Way clusters, namely the Arches (Figer 2005; Martins et al. 2008) and NGC 3603 (Schnurr et al. 2008a). Is this 150 M⊙ limit statistical or physical?

The determination of stellar masses of very massive stars from colour–magnitude diagrams is highly unreliable, plus corrections to present mass estimates need to be applied to estimate initial masses. Studies are further hindered by severe spatial crowding within the cores of these star clusters (Maíz Apellániz 2008). In addition, sophisticated techniques are required to extract the physical parameters of early-type stars possessing powerful stellar winds from optical/infrared spectroscopy (Conti, Crowther & Leitherer 2008). In general, the incorporation of line blanketing has led to a reduction in derived temperatures for O stars from photospheric lines (Puls, Vink & Najarro 2008) while the reverse is true for emission lines in Wolf–Rayet (WR) stars (Crowther 2007).

Here we re-analyse the brightest members of the star cluster (HD 97950) responsible for NGC 3603, motivated in part by the case of A1 which is an eclipsing binary system, whose individual components have been measured by Schnurr et al. (2008a). We also re-analyse the brightest members of R136 (HD 38268) – the central ionizing cluster of the Tarantula nebula (30 Doradus) within the Large Magellanic Cloud (LMC). R136 is sufficiently young (Massey & Hunter 1998) and massive (≤5.5 × 104 M⊙; Hunter, Shaya & Holtzman 1995) within the Local Group of galaxies to investigate the possibility of a physical limit beyond 150 M⊙ (Massey & Hunter 1998; Selman et al. 1999; Oey & Clarke 2005).

As recently as the 1980s, component ‘a’ within R136 was believed to be a single star with a mass of several thousand solar masses (Cassinelli, Mathis & Savage 1981; Savage et al. 1983), although others favoured a dense star cluster (Moffat & Seggewiss 1983). The latter scenario was supported by speckle interferometric observations (Weigelt & Bauer 1985) and subsequently confirmed by Hubble Space Telescope (HST) imaging (Hunter et al. 1995).

Here, we present new analyses of the brightest sources of NGC 3603 and R136. In contrast to classical WR stars, these WN stars are believed to be young, main-sequence stars, albeit possessing very strong stellar winds as a result of their high stellar luminosities (Crowther 2007). Section 2 describes archival HST and Very Large Telescope (VLT) observations, while the spectroscopic analysis is presented in Section 3. Comparisons with contemporary evolutionary models are presented in Section 4, revealing spectroscopic masses in excellent agreement with dynamical masses of the components of NGC 3603A1. For the R136 stars, exceptionally high initial masses of 165–320 M⊙ are inferred. We consider the possibility that these stars are binaries in Section 5, from archival X-ray observations of R136. In Section 6 we simulate star clusters to re-evaluate the stellar upper mass limit, and find that both R136 and NGC 3603 are consistent with a tentative upper limit of ∼300 M⊙. Section 7 considers the contribution of the R136 stars to the global properties of both this cluster and 30 Doradus. Finally, we consider the broader significance of an increased mass limit for stars in Section 8.

2 OBSERVATIONS

Our analysis of WN stars in NGC 3603 and R136 is based upon archival ultraviolet (UV)/optical HST and near-IR VLT spectroscopy, summarized in Table 1, combined with archival high spatial resolution near-IR imaging. We prefer the latter to optical imaging due to reduced extinction, efficient correction for severe spatial crowding (Schnurr et al. 2008a,b) and consistency with studies of other young, high-mass clusters (e.g. Arches; Martins et al. 2008).

Table 1

Log of spectroscopic observations of R136 and NGC 3603 stars used in this study.

Table 1

Log of spectroscopic observations of R136 and NGC 3603 stars used in this study.

2.1 NGC 3603

Our spectroscopic analysis of the three WN 6h systems within NGC 3603 is based upon archival HST/Faint Object Spectroscopy (FOS) from Drissen et al. (1995) plus integral field VLT/Spectrograph for Integral Field Observation in the Near-Infrared (SINFONI) spectroscopy from Schnurr et al. (2008a)– see Table 1 for the log of observations. The spectral resolution of the SINFONI data sets is _R_∼ 3000, with adaptive optics (AO) used to observe A1, B and C, versus _R_∼ 1300 for FOS for which a circular aperture of diameter 0.26 arcsec was used.

Differential photometry from VLT/SINFONI integral field observations (Schnurr et al. 2008a), were tied to unpublished VLT/Infrared Spectrometer and Array Camera (ISAAC) _K_s-band acquisition images from 2002 June 16 (ESO Programme 69.D-0284(A), P.I. Crowther) using the relatively isolated star NGC 3603C. The VLT/ISAAC frames are calibrated against Two Micron All Sky Survey (2MASS) photometry (Skrutskie et al. 2006) using six stars in common (±0.05 mag). These are presented in Table 2 together with absolute magnitudes resulting from a distance of 7.6 ± 0.35 kpc (distance modulus 14.4 ± 0.1 mag) plus an extinction law from Melena et al. (2008) using

where _R_MW_V_= 3.1 for the Milky Way foreground component and _R_NGC 3603_V_= 4.3 for the internal NGC 3603 component. We obtain mag from spectral energy distribution fits, which are consistent with recent determinations of E(B_−_V) = 1.39 mag from Melena et al. (2008). In their analysis, Crowther & Dessart (1998) used a lower overall extinction of E(B_−_V) = 1.23 mag, albeit a higher distance of 10 kpc (distance modulus of 15.0 mag) to NGC 3603.

Table 2

Stars brighter than mag within 13 arcsec (0.5 pc) of NGC 3603A1, together with _V_-band photometry from Melena et al. (2008) in parenthesis.

a For a distance of 7.6 ± 0.35 kpc (distance modulus 14.4 ± 0.1 mag).

b A1 is a 3.77 d double-eclipsing system, while C is a 8.9-d SB1 (Schnurr et al. 2008a).

c An updated classification scheme for Of, Of/WN and WN stars (Walborn & Crowther, in preparation) favours O3 If*/WN 6 for NGC 3603C.

Table 2

Stars brighter than mag within 13 arcsec (0.5 pc) of NGC 3603A1, together with _V_-band photometry from Melena et al. (2008) in parenthesis.

a For a distance of 7.6 ± 0.35 kpc (distance modulus 14.4 ± 0.1 mag).

b A1 is a 3.77 d double-eclipsing system, while C is a 8.9-d SB1 (Schnurr et al. 2008a).

c An updated classification scheme for Of, Of/WN and WN stars (Walborn & Crowther, in preparation) favours O3 If*/WN 6 for NGC 3603C.

2.2 R136

Our UV/optical/infrared spectroscopic analysis of all four hydrogen-rich WN 5h WR stars in R136 (Crowther & Dessart 1998) is also based upon archival HST and VLT spectroscopy – see Table 1 for the log of observations.

Goddard High Resolution Spectrograph (GHRS) ultraviolet observations of R136a stars have been described by de Koter, Heap & Hubeny (1997), achieving a spectral resolution of _R_∼ 2000–3000 using the SSA aperture of 0.22 × 0.22 arcsec. This prevented individual UV spectroscopy for R136a1 and R136a2 which are separated by ∼0.1 arcsec. For context, 0.1 arcsec subtends 5000 au at the distance of the LMC. R136a2 UV spectroscopy was attempted by S. Heap (Programme 5297), but was assessed to be unreliable by de Koter et al. (1997) and so is also excluded here. Note also that R136c has not been observed with GHRS. Visual FOS data sets, with a circular aperture of diameter 0.26 arcsec achieved _R_∼ 1300, although R136a1 and R136a2 once again suffer significant contamination from one another. It is solely at near-IR wavelengths (SINFONI) that R136a1 and R136a2 are spectrally separated, for which R136b served as an AO reference star (Schnurr et al. 2009). The spectral resolution of the SINFONI data sets is _R_∼ 3000.

We employ high spatial resolution _K_s-band photometry of R136. Differential _K_s photometry from AO assisted VLT/SINFONI integral field data sets are tied to identical spatial resolution wider field VLT Multi-Conjugate Adaptive Optics Demonstrator (MAD) imaging (Campbell et al. 2010) using the relatively isolated star R136b (WN 9h). Three overlapping fields were observed with VLT/MAD, for which Field 1 provided the highest quality in R136 [full width at half-maximum (FWHM) ∼ 0.1 arcsec], as shown in Fig. 1, itself calibrated using archival HAWK-I and 2MASS data sets (see Campbell et al. 2010). For 12 stars in common between MAD photometry and HST Near Infrared Camera and Multi Object Spectrometer (NICMOS) F205W imaging (Brandner et al. 2001) transformed into the Cerro Tololo Inter-American Observatory (CTIO) _K_-band system, Campbell et al. (2010) find mag.

Figure 1

VLT MAD _K_s-band 12 × 12 arcsec2 (3 × 3 pc2 for the LMC distance of 49 kpc) image of R136 (Campbell et al. 2010) together with a view of the central 4 × 4 arcsec2 (1 × 1 pc2) in which the very massive WN 5h stars discussed in this paper are labelled (component b is a lower mass WN 9h star). Relative photometry agrees closely with integral field SINFONI observations (Schnurr et al. 2009).

In Table 3 we present _K_s aperture photometry and inferred absolute magnitudes of stars brighter than mag within 20 arcsec (5 pc) of R136a1. Of these only R136a, R136b and R136c components lie within a projected distance of 1 pc from R136a1. Spectral types are taken from Crowther & Dessart (1998) or Walborn & Blades (1997), although other authors prefer alternative nomenclature, e.g. O2 If for both Mk 39 and R136a5 according to Massey et al. (2004, 2005).

Table 3

Stars brighter than mag within 20 arcsec (5 pc) of R136a1, together with _V_-band photometry from Hunter et al. (1995) in parenthesis.

a For a distance modulus of 18.45 ± 0.18 mag (49 ± 4 kpc).

b R136c is variable and X-ray bright, Mk 39 is a 92-d SB1 binary (Massey et al. 2002; Schnurr et al. 2008b) and Mk 34 is a binary according to unpublished Gemini observations (Schnurr et al., in preparation).

c An updated classification scheme for Of, Of/WN and WN stars (Walborn & Crowther, in preparation) favours O2 If* for Mk 42 and R136a5, O3.5 If for Mk 37a and O2 If/WN for Mk 39.

Table 3

Stars brighter than mag within 20 arcsec (5 pc) of R136a1, together with _V_-band photometry from Hunter et al. (1995) in parenthesis.

a For a distance modulus of 18.45 ± 0.18 mag (49 ± 4 kpc).

b R136c is variable and X-ray bright, Mk 39 is a 92-d SB1 binary (Massey et al. 2002; Schnurr et al. 2008b) and Mk 34 is a binary according to unpublished Gemini observations (Schnurr et al., in preparation).

c An updated classification scheme for Of, Of/WN and WN stars (Walborn & Crowther, in preparation) favours O2 If* for Mk 42 and R136a5, O3.5 If for Mk 37a and O2 If/WN for Mk 39.

Interstellar extinctions for the R136 WN 5h stars are derived from UV to near-IR spectral energy distribution fits (R136c lacks UV spectroscopy), adopting foreground Milky Way (LMC) extinctions of (0.06) mag, plus variable internal 30 Doradus nebular extinction. The adopted extinction law follows Fitzpatrick & Savage (1984) as follows:

where _R_MW_V_=_R_LMC_V_= 3.2 and _R_30 Dor_V_= 4.0. We derive mag for the R136a stars and 0.30 mag for R136c and estimate mag for other stars except that Mk 34 (WN 5h) mirrors the higher extinction of R136c. From Table 3, R136c is 0.5 mag fainter than R136a2 in the V band (Hunter et al. 1995) but is 0.06 mag brighter in the _K_s band, justifying the higher extinction. An analysis based solely upon optical photometry could potentially underestimate the bolometric magnitude for R136c with respect to R136a2. The main source of uncertainty in absolute magnitude results from the distance to the LMC. The mean of seven independent techniques (Gibson 2000) suggests an LMC distance modulus of 18.45 ± 0.06 mag. Here we adopt an uncertainty of ±0.18 mag owing to systematic inconsistencies between the various methods.

3 SPECTROSCOPIC ANALYSIS

For our spectroscopic study we employ the non-local thermodynamic equilibrium (LTE) atmosphere code cmfgen (Hillier & Miller 1998) which solves the radiative transfer equation in the comoving frame, under the additional constraints of statistical and radiative equilibrium. Since cmfgen does not solve the momentum equation, a density or velocity structure is required. For the supersonic part, the velocity is parametrized with an exponent of β= 0.8. This is connected to a hydrostatic density structure at depth, such that the velocity and velocity gradient match at this interface. The subsonic velocity structure is defined by a fully line-blanketed, plane-parallel TLUSTY model (Lanz & Hubeny 2003) whose gravity is closest to that obtained from stellar masses derived using evolutionary models, namely log _g_= 4.0 for R136 stars and log _g_= 3.75 for NGC 3603 stars. cmfgen incorporates line blanketing through a superlevel approximation, in which atomic levels of similar energies are grouped into a single superlevel which is used to compute the atmospheric structure.

Stellar temperatures, T*, correspond to a Rosseland optical depth 10, which is typically 1000 to 2000 K higher than effective temperatures _T_2/3 relating to optical depths of 2/3 in such stars.

Our model atom include the following ions: H i, He i–ii, C iii–iv, N iii–v, O iii–vi, Ne iv–v, Si iv, P iv–v, S iv–v, Ar v–vii, Fe iv–vii, Ni v–vii, totalling 1141 superlevels (29 032 lines). Other than H, He, CNO elements, we adopt solar abundances (Asplund et al. 2009) for NGC 3603, which are supported by nebular studies (Esteban et al. 2005; Lebouteiller et al. 2008). LMC nebular abundances (Russell & Dopita 1990) are adopted for R136, with other metals scaled to 0.4 Z⊙, also supported by nebular studies of 30 Doradus (e.g. Peimbert 2003; Lebouteiller et al. 2008). We have assumed a depth-independent Doppler profile for all lines when solving for the atmospheric structure in the comoving frame, while in the observer's frame, we have adopted a uniform turbulence of 50 km s−1. Incoherent electron scattering and Stark broadening for hydrogen and helium lines are adopted. With regard to wind clumping (Hillier 1991), this is incorporated using a radially dependent volume filling factor, f, with _f_∞= 0.1 at _v_∞, resulting in a reduction in mass-loss rate by a factor of ⁠.

3.1 NGC 3603

Fig. 2 compares spectral energy distributions for the WN 6h stars in NGC 3603 with reddened, theoretical models, including the individual components of A1. For this system we obtain mag, such that we adopt Δ_m_=_m_A1a−_m_A1b=−0.43 ± 0.3 mag for the individual components (Schnurr et al. 2008a). We have estimated individual luminosities in two ways. First, on the basis of similar mean molecular weights μ for individual components, stellar luminosities result from _L_∝_M_α, with α∼ 1.5 for zero-age main-sequence (ZAMS) models in excess of 85 M⊙. We derive L(A1a)/L(A1b) = 1.62+0.6−0.4 from the dynamical masses of 116 ± 31 and 89 ± 16 M⊙ obtained by Schnurr et al. (2008a) for A1a and A1b, respectively. Alternatively, we have obtained a luminosity ratio for components within A1 from a fit to the NICMOS light curve (Moffat et al. 2004; E. Antokhina, private communication), using the radii of the Roche lobes and light-curve-derived temperatures. We also obtain L(A1a)/L(A1b) = 1.62 from this approach.

Figure 2

Spectral energy distributions of NGC 3603 WN 6h stars from HST/FOS together with _K_s photometry from VLT/SINFONI and VLT/ISAAC plus reddened theoretical spectral energy distributions (red lines). Components A1a (green lines) and A1b (blue lines) are also shown separately.

Although we have matched synthetic spectra to near-IR photometry, we note that no major differences would be obtained with either HST Wide Field and Planetary Camera 2 (WFPC2) or Advanced Camera for Surveys (ACS) data sets. For example, our reddened spectral energy distribution for NGC 3603B implies _B_= 12.40 mag, which matches HST/ACS photometry to within 0.06 mag (Melena et al. 2008).

We have used diagnostic optical and near-IR lines, including N iii 4634–41, 2.103–2.115 μm, N iv 3478–83, 4058, together with He ii 4686, 2.189 μm plus Brγ. We are unable to employ helium temperature diagnostics since He i lines are extremely weak at optical and near-IR wavelengths. Spectroscopic comparisons with HST/FOS data sets are presented in Fig. 2. Overall, emission features are well reproduced, although absorption components of higher Balmer–Pickering lines are too strong, especially for NGC 3603C. Fits are similar to Crowther & Dessart (1998), in spite of an improved TLUSTY structure within the photosphere. Fortunately, our spectral diagnostics are not especially sensitive to the details of the photosphere since they assess the inner wind conditions. Spectroscopically, we use VLT/SINFONI spectroscopy of A1a and A1b obtained close to quadrature (Schnurr et al. 2008a) to derive mass-loss rates and hydrogen contents. The ratio of He ii 2.189 μm to Brγ provides an excellent diagnostic for the hydrogen content in WN stars, except for the latest subtypes (approximately WN 8 and later). As such, we are able to use primarily optical spectroscopy for the determination of stellar temperatures, with hydrogen content obtained from near-IR spectroscopy. Near-IR comparisons between synthetic spectra and observations are presented for each of the WN 6h stars within NGC 3603 in Fig. 3.

Figure 3

Rectified, spatially resolved near-IR (VLT/SINFONI) spectroscopy of NGC 3603 WN 6h stars (Schnurr et al. 2008a, black), including A1 close to quadrature (A1a blueshifted by 330 km s−1 and A1b redshifted by 433 km s−1), together with synthetic spectra (red, broadened by 50 km s−1). N iii 2.103–2.115 μm and He i 2.112–2.113 μm contribute to the emission feature blueward of Brγ.

Physical properties for these stars in NGC 3603 are presented in Table 4, including evolutionary masses obtained from solar-metallicity non-rotating models. Errors quoted in the table account for both photometric and spectroscopic uncertainties, and represent the range of permitted values. However, for NGC 3603C we assume that the primary dominates the systemic light since the companion is not detected in SINFONI spectroscopy (Schnurr et al. 2008a). Should the companion make a non-negligible contribution to the integrated light, the derived properties set out in Table 4 would need to be corrected accordingly. The improved allowance for metal line blanketing implies ∼10–20 per cent higher stellar temperatures (T*∼ 40 000–44 000 K) and, in turn, larger bolometric corrections (_M_Bol−_MV_∼−3.7 ± 0.2 mag) for these stars than earlier studies (Crowther & Dessart 1998), who adopted a different combination of E(B_−_V) and distance for NGC 3603. We note that Schmutz & Drissen (1999) have previously derived T*∼ 46 000 K for NGC 3603B, resulting in a similar luminosity (log L/L⊙= 6.4) to that obtained here. Analysis of the NICMOS light curve (Moffat et al. 2004) yields T(A1a)/T(A1b) = 1.06, in support of the spectroscopically derived temperature ratio, albeit 7 per cent higher in absolute terms. In view of the underlying assumptions of the photometric and spectroscopic techniques, overall consistency is satisfactory.

Table 4

Physical properties of NGC 3603 WN 6h stars.

a Corresponds to the radius at a Rosseland optical depth of τRoss= 10.

mag for A1, for which we adopt Δ_m_=_m_A1a−_m_A1b=−0.43 ± 0.30 mag (Schnurr et al. 2008a). The ratio of their luminosities follows from their dynamical mass ratios together with _L_∝μ_M_1.5 (and is supported by NICMOS photometry from Moffat et al. 2004).

Table 4

Physical properties of NGC 3603 WN 6h stars.

a Corresponds to the radius at a Rosseland optical depth of τRoss= 10.

mag for A1, for which we adopt Δ_m_=_m_A1a−_m_A1b=−0.43 ± 0.30 mag (Schnurr et al. 2008a). The ratio of their luminosities follows from their dynamical mass ratios together with _L_∝μ_M_1.5 (and is supported by NICMOS photometry from Moffat et al. 2004).

For A1a and B, we adopt abundances of _X_C= 0.008 per cent, _X_N= 0.8 per cent and _X_O= 0.013 per cent by mass, as predicted by evolutionary models for _X_H∼ 60 per cent, versus _X_C= 0.005 per cent, _X_N= 0.6 per cent and _X_O= 0.25 per cent for A1b and C, for which _X_H∼ 70 per cent.

3.2 R136

Theoretical spectral energy distributions of all four WN 5h stars are presented in Fig. 4. Although we have matched synthetic spectra to VLT/MAD + SINFONI photometry, we note that no significant differences would be obtained from WFPC2 imaging (e.g. Hunter et al. 1995), providing these are sufficiently isolated. For example, our reddened spectral energy distribution for R136a3 implies _V_= 13.0 mag, in agreement with F555W photometry. VLT/MAD photometry is preferred to WFPC2 for the very crowded core of R136 (e.g. a1 and a2), although recent Wide Field Camera 3 (WFC3) imaging of R136 achieves a similar spatial resolution at visible wavelengths. Nevertheless, we adhere to VLT/MAD + SINFONI for our primary photometric reference since

Figure 4

Spectral energy distributions of R136 WN 5h stars from HST/FOS together using _K_s photometry from VLT/SINFONI calibrated with VLT/MAD imaging. Reddened theoretical spectral energy distributions are shown as red lines.

• solely VLT/SINFONI spectrally resolves R136a1 from R136a2 (Schnurr et al. 2009). We also note the consistent line to continuum ratios between optical and near-IR diagnostics, ruling out potential late-type contaminants for the latter;
• uncertainties in dust extinction at _K_s are typically 0.02 mag versus ∼0.2 mag in the V band, recalling ⁠. R136c is significantly fainter than R136a2 at optical wavelengths yet is brighter at _K_s and possesses a higher bolometric magnitude (contrast these results with the optical study of Rühling 2008);
• we seek to follow a consistent approach to recent studies of the Arches cluster which necessarily focused upon the _K_s band for photometry and spectroscopy (e.g. Martins et al. 2008).

We estimate terminal wind velocities from optical and near-IR helium lines. UV HST/GHRS spectroscopy has been obtained for the R136a stars (Heap et al. 1994; de Koter et al. 1997), although the severe crowding within this region results in the spectra representing blends of individual components for R136a1 and R136a2 (apparent in Fig. 4). As a result, we favour velocities from FOS and SINFONI spectroscopy for these cases.

In Fig. 5 we present UV/optical/near-IR spectral fits for the WN 5h star R136a3, which is sufficiently isolated that contamination from other members of R136a is minimal in each spectral region (in contrast to R136a1 and R136a2 at UV/optical wavelengths). Consistent optical and near-IR fits demonstrate that no underlying red sources contribute significantly to the _K_-band SINFONI data sets. Diagnostics include O v 1371, S v 1501, N iv 3478–83, 4058, N v 4603–20, 2.10 μm together with He ii 4686, 2.189 μm plus Hα, Hβ, Brγ. For T*≤ 50 000 K both O v 1371 and N v 2.11 μm are underestimated for the R136a WN 5h stars, while for T*≥ 56 000 K, N iv 3478–83 and 4057 become too weak, and S v 1501 becomes too strong. Therefore, we favour T*∼ 53000 ± 3000 K, except that T*∼ 51 000 ± 5000 K is preferred for R136c since UV spectroscopy is not available and N v 2.100 μm is weak/absent. Again, we are unable to employ helium temperature diagnostics since He i lines are extremely weak.

Figure 5

Rectified, ultraviolet (HST/GHRS), visual (HST/FOS) and near-IR (VLT/SINFONI) spectroscopy of the WN 5h star R136a3 together with synthetic UV, optical and near-IR spectra, for T*= 50 000 K (red) and T*= 55 000 K (blue). Instrumental broadening is accounted for, plus an additional rotational broadening of 200 km s−1.

Fig. 4 compares synthetic spectra with HST/FOS spectroscopy, for which emission features are again well matched, albeit with predicted Balmer–Pickering absorption components that are too strong. Fig. 6 presents near-IR spatially resolved spectroscopy of R136 WN stars (Schnurr et al. 2009) together with synthetic spectra, allowing for instrumental broadening (100 km s−1) plus additional rotational broadening of 200 km s−1 for R136a2, R136a3, R136c. As for NGC 3603 stars, the ratio of He ii 2.189 μm to Brγ provides an excellent diagnostic for the hydrogen content for the R136 stars. A summary of the resulting physical and chemical parameters is presented in Table 5, with errors once again accounting for both photometric and spectroscopic uncertainties, representing the range of permitted values, although the primary uncertainty involves the distance to the LMC.

Figure 6

Rectified, spatially resolved near-IR (VLT/SINFONI; Schnurr et al. 2009) spectroscopy of R136 WN 5h stars (black) together with synthetic spectra, accounting for instrumental broadening (100 km s−1) plus rotational broadening of 200 km s−1 for R136a2, R136a3 and R136c. Consistent hydrogen contents are obtained from the peak intensity ratio of Brγ/He ii 2.189 μm and optical (Pickering–Balmer series) diagnostics. Nitrogen emission includes N v 2.100 μm and N iii 2.103–2.115 μm.

Table 5

Physical properties of R136 WN 5h stars.

a Corresponds to the radius at a Rosseland optical depth of τRoss= 10.

b Component R136c is probably a colliding-wind massive binary. For a mass ratio of unity, initial (current) masses of each component would correspond to ∼160 M⊙ (∼130 M⊙).

c d_M_/d_t_Vink relates to Vink et al. (2001) mass-loss rates for _Z_= 0.43 Z⊙.

Table 5

Physical properties of R136 WN 5h stars.

a Corresponds to the radius at a Rosseland optical depth of τRoss= 10.

b Component R136c is probably a colliding-wind massive binary. For a mass ratio of unity, initial (current) masses of each component would correspond to ∼160 M⊙ (∼130 M⊙).

c d_M_/d_t_Vink relates to Vink et al. (2001) mass-loss rates for _Z_= 0.43 Z⊙.

With respect to earlier studies (Heap et al. 1994; de Koter et al. 1997; Crowther & Dessart 1998), the improved allowance for metal line blanketing also infers 20 per cent higher stellar temperatures (T*∼ 53 000 K) for these stars, and in turn, larger bolometric corrections (_M_Bol−_MV_∼−4.6 mag). Such differences, with respect to non-blanketed analyses, are typical of WR stars (Crowther 2007). Similar temperatures (T*= 50 000 or 56 000 K) to the present study were obtained for these WN stars by Rühling (2008) using a grid calculated from the Potsdam line-blanketed atmospheric code (for a summary see Ruehling, Gräfener & Hamann 2008), albeit with 0.1–0.2 dex lower luminosities (log L/L⊙= 6.4–6.7) using solely optically derived (lower) extinctions and absolute magnitudes.

To illustrate the sensitivity of mass upon luminosity, we use the example of R136a5 (O2–3 If/WN) for which we estimate T*∼ 50 000 K using the same diagnostics as the WN 5h stars, in good agreement with recent analyses of O2 stars by Walborn et al. (2004) and Evans et al. (2010b). This reveals a luminosity of log L/L⊙= 6.35, corresponding to an initial mass in excess of 100 M⊙, versus T*∼ 42 500 K, log L/L⊙= 5.9 and ∼65 M⊙ according to de Koter et al. (1997).

Regarding elemental abundances, we use the Hβ to He ii 5412 and Brγ to He ii 2.189 μm ratios to derive (consistent) hydrogen contents, with the latter serving as the primary diagnostic for consistency with NGC 3603 stars (current observations do not include He ii 5412). We adopt scaled solar abundances for all metals other than CNO elements. Nitrogen abundances of _X_N= 0.35 per cent by mass, as predicted by evolutionary models, are consistent with near-IR N v 2.100 μm recombination line observations, while we adopt carbon and oxygen abundances of _X_C=_X_O= 0.004 per cent by mass.

4 EVOLUTIONARY MODELS

We have calculated a grid of main-sequence models using the latest version of the Geneva stellar evolution code for 85, 120, 150, 200, 300 and 500 M⊙. The main-sequence evolution of such high-mass stars does not suffer from stability issues. Although a detailed description is provided elsewhere (Hirschi, Meynet & Maeder 2004), together with recent updates (Eggenberger et al. 2008), models include the physics of rotation and mass loss, which are both crucial to model the evolution of very massive stars. Although many details relating to the evolution of very massive stars remain uncertain, we solely consider the main-sequence evolution here using standard theoretical mass-loss prescription for O stars (Vink, de Koter & Lamers 2001) for which Mokiem et al. (2007) provide supporting empirical evidence. In the models, we consider the onset of the WR phase to take place when the surface hydrogen content _X_H < 30 per cent if _T_eff≥ 10 000 K, during which an empirical mass-loss calibration is followed (Nugis & Lamers 2000). The post-main-sequence evolution is beyond the scope of this study and will be discussed elsewhere.

All the main effects of rotation are included in the calculations: centrifugal support, mass-loss enhancement and especially mixing in radiative zones (Maeder 2009), although predictions for both non-rotating and rotating models are considered here. For rotating models we choose an initial ratio of the velocity to critical (maximum) rotation of _v_init/_v_crit= 0.4, which corresponds to surface equatorial velocities of around 350 km s−1 for the 85 M⊙ model and around 450 km s−1 for the 500 M⊙ case.

We have calculated both solar (_Z_= 1.4 per cent by mass) and LMC (_Z_= 0.6 per cent by mass) metallicities. The evolution of 120 M⊙ models in the Hertzsprung–Russell diagram is presented in Fig. 7. The effects of rotation are significant. Because of additional mixing, helium is mixed out of the core and thus the opacity in the outer layers decreases. This allows rotating stars to stay much hotter than non-rotating stars. Indeed, the effective temperature of the rotating models stay as high as 45 000–55 000 K, whereas the effective temperature of non-rotating models decreases to 20 000–25 000 K. Therefore, rapidly rotating stars progress directly to the classical WR phase, while slow rotators are expected to become η Car-like luminous blue variables (see Meynet & Maeder 2005).

Figure 7

Comparison between main sequence evolutionary predictions for both rotating (_Z_= 0.006: H120z06S400) and non-rotating (_Z_= 0.006: H120z06S000, _Z_= 0.014: H120z14S000) 120 M⊙ models. Horizontal lines in the top left-hand panel correspond to the transition to the WR phase.

Rotating models can reach higher luminosities, especially at very low metallicity (see Langer et al. 2007). This is explained by additional mixing above the convective core. Finally, by comparing models at solar metallicity and LMC metallicity, we can see that lower metallicity models reach higher luminosities. This is due to weaker mass loss at lower metallicity.

Langer et al. (2007) provide 150 M⊙ tracks at 0.2 and 0.05 Z⊙ with high initial rotation velocities of 500 km s−1. From these it is apparent that rapidly rotating very metal-deficient models can achieve high stellar luminosities. According to Langer et al. (2007) a 150 M⊙ model reaches log L/L⊙∼ 6.75 at 0.05 Z⊙ and 6.5 at 0.2 Z⊙. Could the R136 stars represent rapidly rotating stars of initial mass 150 M⊙? We have calculated a Small Magellanic Cloud (SMC) metallicity (0.14 Z⊙) model for a 150 M⊙ star initially rotating at _v_init= 450 km s−1, which achieves log L/L⊙≈ 6.5 after 1.5 Myr and 6.67 after 2.5 Myr. Our models are thus compatible with Langer et al. (2007) described above. We nevertheless exclude the possibility that the R136 WN 5 stars possess initial masses below 150 M⊙ since

• the metallicity of 30 Doradus is a factor of 3 higher than the SMC (e.g. Peimbert 2003; Lebouteiller et al. 2008);
• R136 has an age of less than 2 Myr (de Koter, Heap & Hubeny 1998; Massey & Hunter 1998), since its massive stellar population is analogous to the young star clusters in Car OB1 (1–2 Myr; Walborn 2010) and NGC 3603 (this study);
• clumping-corrected mass-loss rates for the R136 WN stars agree well with LMC metallicity predictions (Vink et al. 2001), which are also supported by studies of O stars in the Milky Way, LMC and SMC (Mokiem et al. 2007).

The surface abundances of a majority of high-mass stars can be well reproduced by models of single stars including the effects of rotation. However, the VLT FLAMES survey has highlighted some discrepancies between models and observations (Hunter et al. 2008), questioning the efficiency of rotation-induced mixing. We are currently investigating this matter by comparing the composition of light elements like boron and nitrogen between models and observations and we find that models including rotation-induced mixing reproduce the abundances of most stars well (Frischknecht et al. 2010).

We can nevertheless consider how a less efficient rotation-induced mixing (or absence of mixing) would affect the conclusions of this paper. If rotation-induced mixing was less efficient than our models predict, the masses derived here would remain at least as high, and would usually be higher. Indeed, less efficient mixing prevents the luminosity to increase as much since less helium is mixed up to the surface and the mean molecular weight, μ, remains lower (_L_≈μ_M_1.5). In particular, less efficient mixing would prevent stars with initial masses around or below 150 M⊙ from reaching such high luminosities as predicted in Langer et al. (2007, see discussion in previous paragraph) and would bring further support to our claim that these stars are more massive than 150 M⊙. The challenge would then be to explain the high effective temperature (_T_eff∼ 50 kK) observed for the R136 stars, which is best explained by rotation induced mixing. Note that less efficient mixing would also make impossible the quasi-chemical evolution of fast rotating (single or binary) stars (Yoon, Langer & Norman 2006) which is currently one of the best scenarios for long/soft gamma-ray bursts progenitors.

4.1 NGC 3603

Fig. 8 compares various observational properties of NGC 3603 WN stars with solar metallicity evolutionary predictions. Non-rotating models imply current masses of 120+26−17 and 92+16−15 M⊙ for A1a and A1b, respectively, at an age of ∼1.5 ± 0.1 Myr. These are in excellent agreement with dynamical mass determinations of 116 ± 31 and 89 ± 16 M⊙ for the primary and secondary A1 components (Schnurr et al. 2008a). Initial and current stellar mass estimates are shown in Table 4, and include a (high) initial mass of 166 ± 20 M⊙ for NGC 3603B. Independent age estimates using pre-main sequence isochrones of low-mass stars also favour low 1 ± 1 Myr ages (Sung & Bessell 2004), while Crowther et al. (2006) estimated 1.3 ± 0.3 Myr for NGC 3603 from a comparison between massive O stars and theoretical isochrones (Lejeune & Schaerer 2001).

Figure 8

Comparison between solar-metallicity (_Z_= 1.4 per cent) models calculated for the main-sequence evolution of 85–200 M⊙ stars [initially rotating at _V_init/_v_crit= 0.4 (dotted) and 0 (solid)] and the physical properties derived from spectroscopic analysis of NGC 3603 WN 6h stars. We obtain excellent agreement with dynamical masses for A1a and A1b for initially non-rotating models at ages of 1.5 ± 0.1 Myr. Current mass-loss rates match solar-metallicity theoretical predictions (Vink et al. 2001) to within 0.2 dex.

4.2 R136

Fig. 9 compares the derived properties of R136a1, R136a2, R136a3 and R136c with LMC metallicity evolutionary predictions, under the assumption that these stars are single. Initial stellar masses in the range 165–320 M⊙ are implied, at ages of 1.7 ± 0.2 Myr, plus high initial rotational rates, in order that the observed surface hydrogen contents of 30–40 per cent by mass are reproduced. Initial and current stellar mass estimates are included in Table 5. Differences in age estimates reflect variations in initial rotation rates. Nevertheless, equatorial rotation rates of _v_e∼ 200 (300) km s−1 are predicted after ∼1.75 (2.75) Myr for a 300 (150) M⊙ star.

Figure 9

Comparison between LMC-metallicity models calculated for the main-sequence evolution of 85–500 M⊙ stars [initially rotating at _v_init/_v_crit= 0.4 (dotted) or 0 (solid)] and the physical properties derived from our spectroscopic analysis. We obtain excellent agreement for initially rapidly rotating, 165–320 M⊙ stars at ages of ∼1.7 ± 0.2 Myr. Current mass-loss rates match LMC-metallicity theoretical predictions (Vink et al. 2001) to within 0.2 dex.

In the absence of photospheric absorption features, the N v 2.100 μm feature provides the best diagnostic of rotational broadening for the R136 stars. This recombination line is intrinsically narrow at the SINFONI spectral resolution (FWHM ∼ 15 Å) because it is formed extremely deep in the stellar wind. Indeed, FWHM ∼ 15 Å is observed for R136a1, as expected either for a non-rotating star, or one viewed pole-on (see Fig. 10). In contrast, R136a2 and R136a3 reveal FWHM ∼ 40 Å, corresponding to _v_e sin _i_≈ 200 km s−1. We are unable to quantify FWHM for R136c since the N v feature is very weak, although it too is consistent with a large rotation rate. Therefore, R136a2, R136a3 and probably R136c show spectroscopic evidence for rapid rotation as shown in Fig. 10, while R136a1 could either be a very slow rotator, or a rapid rotator viewed close to pole-on.

Figure 10

Spectral comparison between VLT/SINFONI spectroscopy of N v 2.10 μm in R136 WN 5 stars (Schnurr et al. 2009) and synthetic spectra (red), allowing for instrumental broadening plus rotational broadening of 200 km s−1 for R136a2, R136a3 and R136c.

For comparison, Wolff et al. (2008) found that R136 lacks slow rotators among lower mass 6–30 M⊙ stars. Wolff et al. derived _v_e sin _i_= 189 ± 23 km s−1 for eleven 15–30 M⊙ stars within R136 versus _v_e sin _i_= 129 ± 13 km s−1 from equivalent mass field stars within the LMC. A much more extensive study of rotational velocities for O stars in 30 Doradus will be provided by the VLT-FLAMES Tarantula Survey (Evans et al. 2010a).

Finally, although we defer the possibility that the WN 5 stars in R136 are (equal mass) binaries until Section 5, it is necessary to remark upon the possibility of chance superpositions within the observed cluster. Line-of-sight effects should not play a significant role in our interpretation of bright systems such as R136a1 (Maíz Apellániz 2008). Chance superposition of other stars has been calculated to contribute at most ∼10 to 20 per cent of the (visible) light (J. Maíz Apellániz, private communication). For example, a 0.2 mag decrease in the absolute _K_-band magnitude of R136a1 arising from the contribution of lower mass stars along this sightline would lead to a 10 per cent reduction in its initial (current) stellar mass, i.e. to 285 M⊙ (235 M⊙).

4.3 Synthetic spectra

We have calculated synthetic spectra for each of the LMC metallicity models, both at the ZAMS and ages corresponding to the surface hydrogen compositions reaching _X_H= 30 per cent for the rotating models. These are presented in Fig. 11 in which mass-loss rates follow the theoretical mass-loss recipes from Vink et al. (2001), both for the case of radially dependent clumped winds with a volume filling factor of _f_∞= 0.1 at _v_∞ (solid) and smooth (_f_∞= 1, dotted) winds. For clumped winds, these result in ZAMS synthetic O supergiant spectra which are equivalent to O3 If (for 85 M⊙) and O2–3 If/WN 5–6 (for 200 M⊙) subtypes. Weaker He iiλ4686 emission would naturally be predicted if we were to instead adopt smooth winds, including O3 V for the 85 M⊙ case and O2 III–If for the 200 M⊙ case at the ZAMS.

Figure 11

Synthetic visual spectra calculated for ZAMS stars (black) of initial mass 85 and 200 M⊙, initially rotating at _v_init/_v_crit= 0 or 0.4, assuming that winds are clumped (_f_∞= 0.1, solid lines) and follow Vink et al. (2001). Synthetic spectra are also shown at the time at which _X_H= 30 per cent by mass is predicted at the surface of the initially rotating model, i.e. 1.75–2.75 Myr (red). Synthetic spectra are broadened by 50 km s−1 (_v_init/_v_crit= 0), _v_e sin _i_= 250 km s−1 (_v_init/_v_crit= 0.4, 0.01 Myr), _v_e sin _i_= 150 km s−1 (_v_init/_v_crit= 0.4, 1.75–2.75 Myr). He iiλ4686 emission is predicted for all phases for these clumpy winds – most prominently during the later phases of initially rapidly rotating very massive stars. If we were to adopt non-clumped winds (_f_∞= 1, dotted lines), He iiλ4686 absorption would be predicted for the ZAMS 85 M⊙ models.

Although specific details depend upon the degree of wind clumping in O star winds and validity of the smooth Vink et al. predictions, stars whose masses exceed a given threshold are expected to display a supergiant signature from the outset. Overall, these results suggest a paradigm shift in our understanding of early O dwarfs, which have hitherto been considered to represent ZAMS stars for the highest mass stars (Walborn et al. 2002).

For the NGC 3603 and R136 WN stars studied here, wind clumping is required to reproduce the electron scattering wings of He iiλ4686 (Hillier 1991). Overall, we find ⁠, in close agreement with Martins et al. (2008) who found that spectroscopic mass-loss rates for the Arches WN 7–9 stars exceeded predictions by +0.2 dex.

Very powerful winds naturally result from the dependence of the mass-loss rate upon the ratio of radiation pressure to gravity, Γe∝L/M. Since _L_∝_M_1.7 for ZAMS stars in the range 30–300 M⊙ based upon evolutionary calculations discussed above, Γe∝_M_0.7. A 30 M⊙ main-sequence star possesses an Eddington parameter of Γe∼ 0.12 (Conti et al. 2008), therefore, a 300 M⊙ star will possess Γe∼ 0.5. The latter achieves Γe∼ 0.7 after 1.5 Myr, typical of the cluster age under investigation here. Γe increases from 0.4 to 0.55 for stars of initial mass 150 M⊙ after 1.5 Myr.

Therefore, the very highest mass stars in very young systems may never exhibit a normal O-type absorption line spectrum. Indeed, O2–3 V stars may rather be limited to somewhat lower ZAMS stars, with a higher mass threshold at lower metallicity. Recall the relatively low dynamical mass of 57 M⊙ for the O3 V primary in R136–38 (Massey, Penny & Vukovich 2002), which currently represents the most massive O2–3 V star dynamically weighed, versus 83 M⊙ (WR20a; Bonanos et al. 2004; Rauw et al. 2005), ≥87 M⊙ (WR21a; Niemela et al. 2008) and 116 M⊙ (NGC 3603A1a; Schnurr et al. 2008a) for some of the most massive H-rich WN stars. O If* or even WN-type emission line spectra may be expected for the highest mass stars within the very youngest giant H ii regions, as is the case for R136 and NGC 3603 here, plus Car OB1 (Smith 2006) and W43 (Blum, Damineli & Conti 1999).

5 BINARITY

NGC 3603A1 and NGC 3603C are confirmed binaries (Schnurr et al. 2008a), and R136c is a probable binary (Schnurr et al. 2009) while component B and R136a1, R136a2, R136a3 are presumed to be single. Of course, we cannot unambiguously confirm that these are genuinely single but if their mass ratios differ greatly from unity the derived physical properties reflect those of the primary. Therefore, here we shall focus upon the possibility that the R136 stars are massive binaries, whose ratio is close to unity (Moffat 2008).

If we adopt similar mean molecular weights, μ, for individual components, current (evolutionary) masses of 150 + 150, 200 + 100 and 220 + 55 M⊙ would be required to match the current properties of R136a1 for mass ratios of 1, 0.5 and 0.25 since _L_∝μ_M_1.5. Of these possibilities, only near-equal mass binaries would contradict the high stellar masses inferred here. If the separation between putative binary components were small, radial velocity variations would be expected. There is no unambiguous evidence for binarity among the stars under investigation, although R136c does exhibit marginal radial velocity variability (Schnurr et al. 2009).

In contrast, the WN 6h system A1 in the young Milky Way cluster NGC 3603 is a short-period massive binary (Schnurr et al. 2008a). However, longer period systems with periods of months to years could have easily avoided spectroscopic detection. We shall employ anticipated properties of colliding wind systems and dynamical interactions to assess whether the R136a stars are realistically long-period binaries. This approach is especially sensitive to equal wind momenta (equal mass) systems since the X-ray emission is maximal for the case of winds whose momenta are equal. A1 in NGC 3603 is relatively faint in X-rays because the components are in such close proximity that their winds collide at significantly below maximum velocity.

5.1 X-rays

If both components in a binary system possessed similar masses, once outflow velocities of their dense stellar winds achieve asymptotic values, collisions would produce stronger X-ray emission than would be expected from a single star. Analytical estimates of X-ray emission from colliding wind systems are approximate (Stevens, Blondin & Pollock 1992). Nevertheless, this approach does provide constraints upon orbital separations to which other techniques are currently insensitive and is especially sensitive to systems whose wind strengths are equal. According to Pittard & Stevens (2002) up to 17 per cent of the wind power of the primary (and 17 per cent of the secondary) can be radiated in X-rays for equal wind momenta systems, versus 0.4 per cent of the primary wind power (56 per cent of the secondary) for (unequal mass) systems whose wind momentum ratio is η= 0.01.

To illustrate the diagnostic potential for X-ray observations, let us consider the expected X-ray emission from R136c under the assumption that it is single. The intrinsic X-ray luminosities of single stars can be approximated by _L_X/_L_Bol∼ 10−7 (Chlebowski, Harnden & Sciortino 1989). From our spectroscopic analysis we would expect _L_X≈ 2 × 1033 erg s−1, yet Chandra imaging reveals an intrinsic X-ray luminosity which is a factor of 30–50 times higher (Portegies Zwart, Pooley & Lewin 2002; Townsley et al. 2006; Guerrero & Chu 2008), arguing for a colliding wind system in this instance. Analytical models favour equal-mass components separated by ∼100 au.

In contrast, the expected X-ray emission from the sum of R136a1, R136a2, R136a3 (and R136a5) – unresolved at Chandra resolution – is _L_X≈ 7 × 1033 erg s−1 under the assumption that they are single. In this case, the intrinsic X-ray emission from R136a is observed to be only a factor of ∼3 higher (Portegies Zwart et al. 2002; Guerrero & Chu 2008) arguing against short-period colliding wind systems from any R136a components. Multiple wind interactions (outside of the binaries but within the cluster) will also produce shocks and X-ray emission (see e.g. Reyes-Iturbide et al. 2009).

Empirically, colliding winds within O-type binaries typically exhibit _L_X/_L_Bol∼ 10−6 (Rauw et al. 2002), although binaries comprising stars with more powerful winds (i.e. WR stars) often possess stronger X-ray emission. Indeed, NGC 3603C has an X-ray luminosity of _L_X≳ 4 × 1034 erg s−1 (Moffat et al. 2002), corresponding to _L_X/_L_Bol≳ 5 × 10−6 based upon our spectroscopically derived luminosity (similar results are obtained for R136c). If we were to assume this ratio for the brightest components of R136a (a1, a2, a3 and a5) under the assumption that each were close binaries – i.e. disregarding predictions from colliding wind theory – we would expect an X-ray luminosity that is a factor of 15 times higher than the observed value. Therefore, the WN stars in R136a appear to be single, possess relatively low-mass companions or have wide separations.

If we make the reasonable assumption that equal-mass stars would possess similar wind properties, we can consider the effect of wind collisions upon the production of X-rays. Let us consider a 150 + 150 M⊙ binary system with a period of 100 d in a circular orbit with separation 3 au, whose individual components each possess mass-loss rates of M⊙ yr−1 and wind velocities of _v_∞= 2600 km s−1. A pair of stars with such properties could match the appearance of R136a1, with individual properties fairly representative of R136a3.

Since the ratio of their wind momenta, η is unity, the fraction, Ξ, of the wind kinetic power processed in the shock from each star is maximal (Ξ∼ 1/6; Pittard & Stevens 2002). One must calculate the conversion efficiency, χ, of kinetic wind power into radiation for each star. From Stevens et al. (1992),

where _v_8 is the wind velocity in units of 108 cm s−1, _d_12 is the distance to the interaction region in units of 1012 cm (i.e. half the separation for equal wind momenta) and is the mass-loss rate in units of 10−7 M⊙ yr−1. We set a lower limit of χ= 1 for cases in which the system radiates all of the collision energy. This allows an estimate of the intrinsic X-ray luminosity (J. Pittard, private communication):

For a 150 + 150 M⊙ system with separation of 3 au, we derive χ∼ 3.7 and predict _L_X= 2.7 × 1036 erg s−1 for each component. The total exceeds the measured X-ray luminosity of R136a by a factor of ∼200 (Townsley et al. 2006; Guerrero & Chu 2008), as shown in Table 6. Recall that if R136a1 were to be a massive binary with a mass ratio of 0.25, the secondary would possess a mass of ∼55 M⊙. This would contribute less than 10 per cent of the observed spectrum, with a secondary to primary wind momentum ratio of η∼ 0.05. Predictions for this scenario are also included in Table 6, revealing X-ray luminosities a factor of ∼6 lower than the equal wind (equal mass) case.

Table 6

Predicted X-ray luminosities from R136a1 for scenarios in which it is either a colliding wind binary using analytical predictions from Stevens et al. (1992) and Pittard & Stevens (2002), or single/multiple for various empirical _L_X/_L_Bol values. Analytical calculations adopt component separations of _d_= 3, 30, 300 au, for a mass ratio (i) _q_= 1 (150 + 150 M⊙) and equal wind momenta, yr−1, _v_∞,1= 2600 km s−1); (ii) _q_= 0.25 (220 + 55 M⊙) and yr−1, _v_∞,1= 2600 km s yr−1, _v_∞,2= 2000 km s−1). The intrinsic X-ray luminosity of R136a is taken from Guerrero & Chu (2008).

Table 6

Predicted X-ray luminosities from R136a1 for scenarios in which it is either a colliding wind binary using analytical predictions from Stevens et al. (1992) and Pittard & Stevens (2002), or single/multiple for various empirical _L_X/_L_Bol values. Analytical calculations adopt component separations of _d_= 3, 30, 300 au, for a mass ratio (i) _q_= 1 (150 + 150 M⊙) and equal wind momenta, yr−1, _v_∞,1= 2600 km s−1); (ii) _q_= 0.25 (220 + 55 M⊙) and yr−1, _v_∞,1= 2600 km s yr−1, _v_∞,2= 2000 km s−1). The intrinsic X-ray luminosity of R136a is taken from Guerrero & Chu (2008).

Stevens et al. (1992) also provide an expression for the characteristic intrinsic column density in colliding wind binaries, (their equation 11), from which would be inferred for R136a1 is it were an equal mass binary system, separated by 3 au. This is an order of magnitude higher than the estimate for R136a (Townsley et al. 2006; Guerrero & Chu 2008), arising from a combination of foreground and internal components.

Of course, a significant fraction of the shock energy could produce relativistic particles rather than X-ray emission (e.g. Pittard & Dougherty 2006). Formally, an equal mass binary system whose components are separated by 600 au is predicted to produce an X-ray luminosity of _L_X= 2.7 × 1034 erg s−1, which is comparable to the total intrinsic X-ray luminosity of R136a (Guerrero & Chu 2008). If 90 per cent of the shock energy were to contribute to relativistic particles, the same X-ray luminosity could arise from a system whose components were separated by 60 au. Still, the relatively low X-ray luminosity of R136a suggests that if any of the R136a WN 5 stars are composed of equal mass binaries, their separations would need to be in excess of 100–200 au. How wide could very massive binaries be within such a dense cluster? To address this question we now consider dynamical interactions.

5.2 Dynamical interactions

The binding energy of a ‘hard’100 + 100 M⊙ binary with a separation of 100 au is comparable to the binding energy of a massive cluster (1041 J). Such a system will have frequent encounters with other massive stars or binaries. The collision rate for a system of separation, a, in a cluster of number density, n, and velocity dispersion, σ, is

where θ is the Safranov number, indicating the importance of gravitational focusing (Binney & Tremaine 1987). Obviously, an encounter with a low-mass star will not affect a very massive binary, so we assume only that encounters with stars in excess of 50 M⊙ will significantly affect a very massive binary.

For a 100 + 100 M⊙ binary with 100 au separation in a cluster with a velocity dispersion of σ= 5 km s−1, initially containing a number density of 500 pc−3 (_n_= 1.5 × 10−47 m−3) then θ= 35. Each very massive binary should have an encounter every 1.8 Myr. For binaries wider than 100 au the encounter rate will be even higher, since the encounter rate scales with the square of the separation, i.e. encounters every ∼0.2 Myr for separations of 300 au or 0.05 Myr for separations of 600 au.

Such close encounters will create an unstable multiple system that will rapidly decay by ejection of the lowest mass star (Anosova 1986). This would not necessarily destroy the binary. Instead, it will harden the system in order to gain the energy required to eject the other star. Thus very massive systems separated by more than ∼100 au will reduce their separations through dynamical interactions with other high-mass stars within the dense core of R136a.

Should one of the high-mass binaries in R136a elude hardening in the way outlined above, would it still escape detection if it remained at a large separation? One binary with components of 150 + 150 M⊙ separated by ∼300 au for which ∼30 per cent of the shock energy contributes to the X-ray luminosity, together with intrinsic X-ray luminosities from single R136a members, obeying _L_X/_L_Bol∼ 10−7 (Chlebowski et al. 1989) would indeed mimic the observed R136a X-ray luminosity (Guerrero & Chu 2008).

In summary, if we adopt similar ratios of X-ray to bolometric luminosities as for R136c and NGC 3603C, we would expect an X-ray luminosity from R136a that is a factor of 15 times higher than the observed value if R136a1, R136a2 and R136a3 were each colliding wind systems. Alternatively, we have followed the colliding wind theory of Stevens et al. (1992) and Pittard & Stevens (2002), and conservatively assume 30 per cent of the shock energy contributes to the X-ray luminosity. If all the very massive stars within the Chandra field-of-view (R136a1, R136a2, R136a3 and R136a5) were members of equal mass binaries with separations of 300 au, we would expect an X-ray luminosity that is over a factor of 2 higher than observed. In view of dynamical effects, at most one of the WN components of R136a might be a long period, large separation (≥300 au) equal mass binary. We cannot rule out short-period, highly unequal-mass binary systems of course, but such cases would have little bearing upon our derived stellar masses.

6 CLUSTER SIMULATIONS

We now compare the highest mass stars in R136 and NGC 3603 with cluster predictions, randomly sampled from the stellar IMF, for a range of stellar mass limits. For R136, the stellar mass within a radius ∼4.7 pc from R136a1 has been inferred for ≥2.8 M⊙ at 2.0 × 104 M⊙ (Hunter et al. 1995). If the standard IMF (Kroupa 2002) is adopted for <2.8 M⊙ a total cluster mass of ≤5.5 × 104 M⊙ is inferred. Andersen et al. (2009) estimated a factor of 2 higher mass by extrapolating their measured mass down to 2.1 M⊙ within 7 pc with a Salpeter slope to 0.5 M⊙. These are upper limits to the mass of R136 itself since the mass function in the very centre of R136 cannot be measured (Maíz Apellániz 2008), such that they include the associated halo of older (≥3 Myr) early-type stars viewed in projection – indeed Andersen et al. (2009) obtained an age of 3 Myr for this larger region. Examples of such stars include R136b (WN 9h), R134 (WN 6(h)) and Melnick 33Sb (WC 5).

For NGC 3603, we adopt a cluster mass of 1.0 × 104 M⊙, the lower limit obtained by Harayama, Eisenhauer & Martins (2008) from high-resolution VLT NAOS/CONICA imaging. These are presented in Table 7, together with the Arches cluster, for which Figer (2008) estimated a cluster mass of 2 × 104 M⊙ based on the mass function of Kim et al. (2006).

Table 7

Comparison between the properties of the NGC 3603, Arches and R136 star clusters and their most massive stars, _M_Max.

Table 7

Comparison between the properties of the NGC 3603, Arches and R136 star clusters and their most massive stars, _M_Max.

6.1 Cluster populations

We simulate a population of clusters and stars by randomly sampling first from a power-law cluster mass function (CMF), and then populating each cluster with stars drawn randomly from a stellar IMF (Parker & Goodwin 2007). Cluster masses are selected from a CMF of the form N(M) ∝_M_−β, with standard slope β= 2 (Lada & Lada 2003) between cluster mass limits 50 and 2 × 105 M⊙. These limits enable the full range of cluster masses to be sampled, including those with masses similar to that of R136. The total mass of clusters is set to 109 M⊙ to fully sample the range of cluster masses. Each cluster is populated with stars drawn from a three-part IMF (Kroupa 2002) of the form

where _m_0= 0.02 M⊙, _m_1= 0.1 M⊙ and _m_2= 0.5 M⊙. We use three different values for _m_3 in the simulations: 150, 300 and 1000 M⊙. Stellar mass is added to the cluster until the total mass is within 2 per cent of the cluster mass. If the final star added to the cluster exceeds this tolerance, then the cluster is entirely repopulated (Goodwin & Pagel 2005).

Our random sampling of the IMF allows low-mass clusters to be composed of one massive star and little other stellar material. This contravenes the proposed fundamental cluster mass–maximum stellar mass (CMMSM) relation (Weidner & Kroupa 2006; Weidner et al. 2010). However, the average maximum stellar mass for a given cluster mass closely follows the CMMSM relation (Parker & Goodwin 2007; Maschberger & Clarke 2008). The results for our three Monte Carlo runs are shown in Fig 12. They reveal that the average relation between the maximum stellar mass and cluster mass (in the 102 and 104 M⊙ interval) is recovered (Weidner et al. 2010), without the constraint that the cluster mass governs the maximum possible stellar mass (Parker & Goodwin 2007). We also indicate 25 and 75 per cent quartiles for the instances of a 104 M⊙ (NGC 3603), 2 × 104 M⊙ (Arches) and 5 × 104 M⊙ (R136) cluster.

6.2 R136

If we assume that R136a1, R136a2 and R136c are either single or the primary dominates the optical/IR light, we obtain an average of 260 M⊙ for their initial masses. In reality this value will be an upper limit due to binarity and/or line-of-sight effects. Nevertheless, from Fig. 12, we would expect the average of the three most massive stars in a cluster of mass 5 × 104 M⊙ to be 150, 230 and 500 M⊙ for an adopted upper mass limit of _m_3= 150, 300 and 1000 M⊙, respectively. Therefore, an upper limit close to 300 M⊙ is reasonably consistent with the stellar masses derived (see also Oey & Clarke 2005).

Figure 12

Average mass of the highest mass star (solid black line) and 25 and 75 per cent quartiles for 104, 2 × 104 and 5 × 104 M⊙ (vertical error bars), plus the average of the three highest mass stars (dotted red line) versus cluster mass. Clusters were randomly populated by sampling a three-part IMF (Kroupa 2002) using a Monte Carlo simulation, under the assumption that the upper mass limit is (left) 150 M⊙; (middle) 300 M⊙; (right) 1000 M⊙. In order to minimize sampling problems we take an ‘average’ cluster as the mean of 1000 realizations of a cluster of a particular mass.

In the lower panel of Fig. 13 we therefore present typical mass functions for a cluster of mass 5 × 104 M⊙ for an adopted upper limit of _m_3= 300 M⊙. From a total of 1.05–1.10 × 105 stars we would expect ∼14 initially more massive than 100 M⊙ to have formed within 5 pc of R136a1, since this relates to the radius used by Hunter et al. (1995) to derive the cluster mass. Indeed, of the 12 brightest near-IR sources within 5 pc of R136a1, we infer initial evolutionary masses in excess of 100 M⊙ for 10 cases (all entries in Table 3 except for R136b and R134). Of course, this comparison should be tempered by (a) contributions of putative secondaries to the light of the primary in case of binarity; (b) dynamical ejection during the formation process (e.g. Brandl et al. 2007). Indeed, Evans et al. (2010b) propose that an O2 III-If* star from the Tarantula survey is a potential high velocity runaway from R136, in spite of its high stellar mass (∼90 M⊙).

Figure 13

Theoretical mass functions for clusters with stellar masses of (upper) 104 M⊙ and (lower) 5 × 104 M⊙, for the scenario with an upper mass limit of _m_3= 300 M⊙. Three and 14 stars with initial masses in excess of 100 M⊙ are anticipated, respectively. Our Monte Carlo approach (Parker & Goodwin 2007) breaks down for very massive stars in the upper panel due to small number statistics.

Overall, R136 favours a factor of ≈2 higher stellar mass limit than is currently accepted. We now turn to NGC 3603 to assess whether it supports or contradicts this result.

6.3 NGC 3603

On average, the highest initial stellar mass expected in a 104 M⊙ cluster would be ∼120 M⊙ with an upper mass limit of _m_3= 150 M⊙, ∼200 M⊙ for _m_3= 300 M⊙ and ∼400 M⊙ for M3= 103 M⊙. Since we infer an initial mass of 166+20−20 M⊙ for component B, a stellar limit intermediate between 150 and 300 M⊙ might be expected, if it is either single or dominated by a single component.1 However, no stars initially more massive than 150 M⊙ would be expected in 25 per cent of 104 M⊙ clusters for the case of _m_3= 300 M⊙ (recall middle panel of Fig. 12).

In the upper panel of Fig. 13 we present typical mass functions for a cluster of mass 1 × 104 M⊙ for an adopted upper mass limit of _m_3= 300 M⊙. We would expect three stars initially more massive than 100 M⊙, versus four observed in NGC 3603 (Table 2). If we were to extend the upper limit of the IMF to _m_3= 1000 M⊙, then the average of the three most massive stars of a 104 M⊙ cluster would be 230 M⊙ which is not supported by NGC 3603 (150 M⊙).

Therefore, the very massive stars inferred here for both R136 and NGC 3603 are fully consistent with a mass limit close to _m_3= 300 M⊙, both in terms of the number of stars initially exceeding 100 M⊙ and the most massive star itself. Before we conclude this section, let us now consider the Arches cluster.

6.4 Arches

The upper mass limit of ∼150 M⊙ from Figer (2005) was obtained for the Arches cluster from photometry. This was confirmed from spectroscopic analysis by Martins et al. (2008) who estimated initial masses of 120–150 M⊙ range for the most luminous (late WN) stars. How can these observations be reconciled with our results for NGC 3603 and R136?

Let us consider each of the following scenarios:

• the Arches cluster is sufficiently old (2.5 ± 0.5 Myr; Martins et al. 2008) that the highest mass stars have already undergone core collapse;
• the Arches is unusually deficient in very massive stars for such a high mass cluster; or
• previous studies have underestimated the masses of the highest mass stars in the Arches cluster.

We have compared the derived properties for the two most luminous WN stars in the Arches cluster (F6 and F9) from Martins et al. (2008) to the solar metallicity grids from Section 4. These are presented in Table 8, revealing initial masses of ≥115+30−25 M⊙ from comparison with non-rotating models. Uncertainties from solely _K_-band spectroscopy are significantly higher than our (UV)/optical/_K_-band analysis of WN stars in NGC 3603 and R136. Masses represent lower limits since Martins et al. (2008) indicated that the metallicity of the Arches cluster is moderately supersolar.

Table 8

Comparison between stellar properties (and inferred masses) of the most luminous Arches stars from Martins et al. (2008) obtained using different photometry, foreground interstellar extinctions and distances to the Galactic Centre.

*Includes –0.25 mag offset to bolometric correction relative to F. Martins (private communication), as reported in Clark, Crowther & Mikles (2009).

‡NICMOS F205W and F160W filters.

Table 8

Comparison between stellar properties (and inferred masses) of the most luminous Arches stars from Martins et al. (2008) obtained using different photometry, foreground interstellar extinctions and distances to the Galactic Centre.

*Includes –0.25 mag offset to bolometric correction relative to F. Martins (private communication), as reported in Clark, Crowther & Mikles (2009).

‡NICMOS F205W and F160W filters.

Core hydrogen exhaustion would be predicted to occur after ∼2.5 Myr for the non-rotating 150 M⊙ solar metallicity models, with core-collapse SN anticipated prior to an age of 3 Myr. However, comparison between the properties of F6 and F9 – the most luminous WN stars in the Arches cluster – with our solar metallicity evolutionary models suggest ages of ∼2 Myr. On this basis, it would appear that scenario (i) is highly unlikely.

Let us therefore turn to scenario (ii), on the basis that the initial mass of the most massive star within the Arches cluster was 150 M⊙. On average, the highest mass star in a 2 × 104 M⊙ cluster would be expected to possess an initial mass in excess of 200 M⊙ for an upper mass limit of _m_3= 300 M⊙. From our simulations, the highest mass star within such a cluster spans a fairly wide range (recall Fig. 12). However, an upper limit of 120 (150) M⊙ would be anticipated in only 1 (5) per cent of cases if we were to adopt _m_3= 300 M⊙. Could the Arches cluster be such a statistical oddity?

Our simulations indicate that a 2 × 104 M⊙ cluster would be expected to host six stars initially more massive than 100 M⊙ for _m_3= 300 M⊙. Our solar metallicity evolutionary models combined with spectroscopic results from Martins et al. (2008) suggest that the five most luminous WN stars (F1, F4, F6, F7, F9), with log L/L⊙≥ 6.15, are consistent with initial masses of ≥100 M⊙, providing they are single or one component dominates their near-IR appearance.

Finally, let us turn to scenario (iii), namely that the stellar masses of the Arches stars have been underestimated to date. Martins et al. (2008) based their analysis upon HST/NICMOS photometry from Figer et al. (2002) together with a (low) Galactic Centre distance of 7.6 kpc from Eisenhauer et al. (2005) plus a (low) foreground extinction of mag from Stolte et al. (2002). Let us also consider the resulting stellar parameters and mass estimates on the basis of the standard Galactic Centre distance of 8 kpc (Reid 1993), more recent (higher) foreground extinction of mag from Kim et al. (2006), plus VLT/NACO photometric results from Espinoza, Selman & Melnick (2009).2 These yield absolute _K_-band magnitudes that are 0.55–0.7 mag brighter than Figer et al. (2002), corresponding to 0.22–0.28 dex higher bolometric luminosities. With respect to the hitherto 150 M⊙ stellar mass limit identified by Figer (2005), these absolute magnitude revisions would conspire to increasing the limit to ≥200 M⊙.

In Table 8 we provide the inferred properties of F6 and F9 on the basis of these photometric properties plus the stellar temperatures and _K_-band bolometric corrections derived by Martins et al. (2008). Our solar metallicity non-rotating models indicate initial masses of ≥185+75−45 and ≥165+60−40 M⊙ for F6 and F9, respectively. These are lower limits to the actual initial masses, in view of the supersolar metallicity of the Arches cluster (Martins et al. 2008). In total, ∼five stars are consistent with initial masses in excess of ≃150 M⊙ (those listed above), plus a further five stars for which initial masses exceed ∼100 M⊙ (F3, F8, F12, F14, F15).

Overall, ten stars initially more massive than ∼100 would suggest either that the mass of the Arches cluster approaches 3 × 104 M⊙ or several cases are near-equal mass binaries (Lang et al. 2005; Wang, Dong & Lang 2006). As such, we would no longer require that the Arches cluster is a statistical oddity, since the highest mass star would not be expected to exceed 200 M⊙ in 25 per cent of 2 × 104 M⊙ clusters.

In conclusion, we have attempted to reconcile the properties of the highest mass stars in Arches cluster with an upper mass limit of _m_3= 300 M⊙. Based upon the Martins et al. (2008) study, comparison with evolutionary models suggests that the initial masses of the most massive stars are ≥115+30−25 M⊙ with ages of ∼2 Myr. We find that an upper mass of 120 M⊙ would be expected in only 1 per cent of 2 × 104 M⊙ clusters for which _m_3= 300 M⊙. However, use of contemporary near-IR photometry and foreground extinctions towards the Arches cluster together with the standard 8 kpc Galactic Centre distance reveal an initial mass of ≥185+75−45 M⊙ for the most massive star based on solar metallicity evolutionary models, with four–five stars consistent with initial masses ≥150 M⊙. An upper mass of 200 M⊙ would be expected in 25 per cent of cases. Robust inferences probably await direct dynamical mass determinations of its brightest members, should they be multiple.

7 GLOBAL PROPERTIES OF R136

We will re-evaluate the ionizing and mechanical wind power resulting from all hot, luminous stars in R136 and NGC 3603 elsewhere (Doran et al., in preparation). Here, we shall consider the role played by the very massive WN stars in R136. We have updated the properties of early-type stars brighter than _MV_=−4.5 mag within a radius of ∼5 pc from R136a1 (Crowther & Dessart 1998) to take account of contemporary _T_eff–spectral type calibrations for Galactic stars (Conti et al. 2008), and theoretical mass-loss rates (Vink et al. 2001) for OB stars. Until a census of the early-type stars within R136 is complete, we adopt O3 subtypes for those stars lacking spectroscopy (Crowther & Dessart 1998). Aside from OB stars and the WN 5h stars discussed here, we have included the contribution of other emission-line stars, namely O2–3 If/WN stars for which we adopt identical temperatures to those of the WN 5 stars, one other WN 5h star (Melnick 34), a WN 6(h) star (R134) for which we estimate T*≈ 42 000 K and a WN 9h star (R136b) for which we estimate T*≈ 35 000 K from its HST/FOS and VLT/SINFONI spectroscopy. Finally, a WC 5 star (Melnick 33Sb) lies at a projected distance of 2.9 pc for which we adopt similar wind properties and ionizing fluxes to single WC 4 stars (Crowther et al. 2002).

Fig. 14 shows the integrated Lyman continuum ionizing fluxes and wind power from all early-type stars in R136, together with the explicit contribution of R136a1, R136a2, R136a3 and R136c, amounting to, respectively, 46 and 35 per cent of the cluster total. If we were to adopt a 2000 K systematically higher temperature calibration for O-type stars in the LMC, as suggested by some recent results (Mokiem et al. 2007; Massey et al. 2009), the contribution of WN 5h stars to the integrated Lyman continuum flux and wind power is reduced to 43 and 34 per cent, respectively.

Figure 14

Upper panel: incremental Lyman continuum ionizing radiation from early-type stars within 5 pc from R136a1 (Crowther & Dessart 1998) either following the Galactic O subtype temperature calibration (solid black), or systematically increasing the temperature calibration by 2000 K for LMC stars (dotted line). The four very luminous WN 5h stars discussed here (dashed red line) contribute 43–46 per cent of the total. Lower panel: as above, except for the mechanical wind power from early-type stars within 5 pc of R136a1, based upon theoretical wind prescriptions (Vink et al. 2001) to which the four WN 5h stars contribute 34 per cent of the total.

We have compared our empirical results with population synthesis predictions (Leitherer et al. 1999) for a cluster of mass 5.5 × 104 M⊙ calculated using a standard IMF (Kroupa 2002) and evolutionary models up to a maximum limit of 120 M⊙. These stars alone approach the mechanical power and ionizing flux predicted for a R136-like cluster at an age of ∼1.7 Myr (Leitherer et al. 1999). Indeed, R136a1 alone provides the Lyman continuum output of seventy O7 dwarf stars (Conti et al. 2008), supplying 7 per cent of the radio-derived ∼1052 photon s−1 ionizing flux from the entire 30 Doradus region (Mills, Turtle & Watkinson 1978; Israel & Koornneef 1979). Improved agreement with synthesis models would be expected if evolutionary models allowing for rotational mixing were used (Vazquez et al. 2007), together with contemporary mass-loss prescriptions for main-sequence stars.

8 DISCUSSION AND CONCLUSIONS

If very massive stars – exceeding the currently accepted 150 M⊙ limit – were to exist in the local universe, they would be

• located in high mass (≥104 M⊙), very young (≤2 Myr) star clusters;
• visually the brightest stars in their host cluster, since _L_∝_M_1.5 for ZAMS stars above 85 M⊙,3 and surface temperatures remain approximately constant for the first 1.5 Myr;
• possess very powerful stellar winds, as a result of the mass dependence of the Eddington parameter, Γe∝L/_M_≈_M_0.5 for such stars. 150–300 M⊙ ZAMS stars, for which Γe≈ 0.4–0.55, would likely possess an O supergiant morphology, while a WR appearance would be likely to develop within the first 1–2 Myr, albeit with significant residual surface hydrogen.

In this study we have presented spectroscopic analyses of bright WN stars located within R136 in the LMC and NGC 3603 in the Milky Way that perfectly match such anticipated characteristics. The combination of line blanketed spectroscopic tools and contemporary evolutionary models reveals excellent agreement with dynamical mass determinations for the components of A1 in NGC 3603, with component B possessing a higher initial mass of ∼170 M⊙, under the assumption that it is single. Application to the higher temperature (see also Rühling 2008) brighter members of R136 suggests still higher initial masses of 165–320 M⊙. Owing to their dense stellar winds, they have already lost up to 20 per cent of their initial mass with the first ∼1.5 Myr of their main-sequence evolution (see also de Koter et al. 1998). The R136 WN 5 stars are moderately hydrogen depleted, to which rotating models provide the best agreement, whereas the NGC 3603 WN 6 stars possess normal hydrogen contents, as predicted by non-rotating models at early phases. These differences may arise from different formation mechanisms (e.g. stellar mergers). Wolff et al. (2008) found a deficit of slow rotators within a sample of R136 early-type stars, although statistically significant results await rotational velocities from the VLT-FLAMES Tarantula Survey (Evans et al. 2010a).

We have assessed the potential that each of the R136 stars represents unresolved equal-mass binary systems. SINFONI spectroscopic observations argue against close, short-period binaries (Schnurr et al. 2009). R136c possesses an X-ray luminosity which is a factor of ∼100 times greater than that which would be expected from a single star and is very likely a massive binary. In contrast, X-ray emission from R136a exceeds that expected from single stars by only a factor of ∼3, to which multiple wind interactions within the cluster will also contribute. If the R136a stars possessed similar X-ray properties to R136c and NGC 3603C, its X-ray emission would be at least a factor of 15 times higher than the observed luminosity. At most, only one of the WN sources within R136a might be a very long period, large separation (∼300 au), equal-mass binary system if as little as 30 per cent of the shock energy contributes to the X-ray luminosity. Dynamical effects would harden such systems on ≪Myr time-scales. We cannot rule out shorter period, unequal-mass binary systems of course, but such cases would have little bearing upon our mass limit inferences.

In Table 9 we present a compilation of stars in R136 (30 Doradus), NGC 3603 and the Arches cluster whose initial masses challenge the currently accepted upper mass limit of ∼150 M⊙. Although the formation of high-mass stars remains an unsolved problem in astrophysics (Zinnecker & Yorke 2007), there are no theoretical arguments in favour of such a limit at 150 M⊙ (e.g. Klapp, Langer & Fricke 1987) – indeed Massey & Hunter (1998) argued against an upper stellar limit based on R136 itself. Observations of the Arches cluster provide the primary evidence for such a sharp mass cut-off (Figer 2005). However, as Table 9 illustrates, contemporary photometry and foreground extinction towards the cluster coupled with the spectroscopy results from Martins et al. (2008) suggest four–five stars initially exceed ≈150 M⊙, with an estimate of ≥185 M⊙ for the most luminous star. Recall Martins et al. (2008) used identical spectroscopic tools to those employed in the current study, plus near-IR spectroscopic observations which also form a central component of our study. On this basis the Arches cluster would no longer be a statistical oddity.

Table 9

Compilation of stars within R136/30 Dor, NGC 3603 and the Arches cluster whose initial masses exceed ≈150 M⊙ according to evolutionary models presented here. Photometry, extinctions and bolometric corrections are presented in this study, with the exception of the Arches cluster for which we follow Espinoza et al. (2009), Kim et al. (2006) and Martins et al. (2008), respectively. Known or suspected binaries are marked with *.

Table 9

Compilation of stars within R136/30 Dor, NGC 3603 and the Arches cluster whose initial masses exceed ≈150 M⊙ according to evolutionary models presented here. Photometry, extinctions and bolometric corrections are presented in this study, with the exception of the Arches cluster for which we follow Espinoza et al. (2009), Kim et al. (2006) and Martins et al. (2008), respectively. Known or suspected binaries are marked with *.

Monte Carlo simulations for various upper stellar mass limits permit estimates of the revised threshold from which ∼300 M⊙ is obtained. It may be significant that both NGC 3603 and R136 are consistent with an identical upper limit, in spite of their different metallicities. Oey & Clarke (2005) obtained an upper limit of ≪500 M⊙, primarily from a maximum stellar mass of 120–200 M⊙ inferred for stars within R136a itself at that time.

Would there be any impact of an upper mass limit of order ∼300 M⊙ on astrophysics? Population synthesis studies are widely applied to star-forming regions within galaxies, for which an upper mass limit of ∼120 M⊙ is widely adopted (Leitherer et al. 1999). A number of properties are obtained from such studies, including star formation rates, enrichment of the local interstellar medium (ISM) through mechanical energy through winds, chemical enrichment and ionizing fluxes. A higher stellar mass limit would increase their global output, especially for situations in which rapid rotation leads to stars remaining at high stellar temperatures.

Could the presence of very high mass stars, such as those discussed here, be detected in spatially unresolved star clusters? The high luminosities of such stars inherently lead to the development of very powerful stellar winds at early evolutionary phases (1–2 Myr). Therefore the presence of such stars ensures that the integrated appearance of R136a (and the core of NGC 3603) exhibits broad He iiλ4686 emission, which at such early phases would not be the case for a lower stellar limit. Other high-mass clusters, witnessed at a sufficiently early age, would also betray the presence of very massive stars through the presence of He ii emission lines at λ1640 and λ4686. Hitherto, such clusters may have been mistaken for older clusters exhibiting broad helium emission from classical WR stars.

Finally, how might such exceptionally massive stars end their life? This question has been addressed from a theoretical perspective (Heger et al. 2003), suggesting neutron star remnants following Type Ib/c core-collapse SNe close to solar metallicities, with weak SNe and black hole remnants at LMC compositions. Rapid rotation causes evolution to proceed directly to the classical WR phase, whereas slow rotators, such as the NGC 3603 WN 6h stars, will likely produce a η Car-like luminous blue variable phase. Extremely metal-deficient stars exceeding ∼140 M⊙ are thought to end their lives prior to core-collapse (Bond, Arnett & Carr 1984). They would undergo an electron–positron pair-instability explosion during the advanced stages that would trigger the complete disruption of the star (Heger & Woosley 2002).

Until recently, such events have been expected from solely metal-free (Population III) stars, since models involving the collapse of primordial gas clouds suggest preferentially high characteristic masses as large as several hundred M⊙ (Bromm & Larson 2004). Langer et al. (2007) have demonstrated that local pair-instability SNe could be produced either by slow rotating moderately metal-poor (≤1/3 Z⊙) yellow hypergiants with thick hydrogen-rich envelopes – resembling SN 2006gy (Ofek et al. 2007; Smith et al. 2007) – or rapidly rotating very metal-deficient (≤10−3 Z⊙) WR stars. Therefore, it is unlikely any of the stars considered here are candidate pair-instability SNe. Nevertheless, the potential for stars initially exceeding 140 M⊙ within metal-poor galaxies suggests that such pair-instability SNe could occur within the local universe, as has recently been claimed for SN 2007bi (Gal-Yam et al. 2009, see also Langer 2009).

Finally, close agreement between our spectroscopically derived mass-loss rates and theoretical predictions allows us to synthesize the ZAMS appearance of very massive stars. We find that the very highest mass progenitors will possess an emission-line appearance at the beginning of their main-sequence evolution due to their proximity to the Eddington limit. As such, spectroscopic dwarf O2–3 stars are not anticipated with the very highest masses. Of course the only direct means of establishing stellar masses is via close binaries. For the LMC metallicity, Massey et al. (2002) have obtained dynamical masses of a few systems, although spectroscopic and photometric searches for other high-mass, eclipsing binaries in 30 Doradus are in progress through the VLT-FLAMES Tarantula Survey (Evans et al. 2010a) and other studies (Schnurr et al., in preparation). Nevertheless, it is unlikely that any other system within 30 Doradus or indeed the entire Local Group of galaxies will compete in mass with the brightest components of R136 discussed here.

1

We do not attempt to claim that NGC 3603B is single, even though it is photometrically and spectroscopically stable and is not X-ray bright (Moffat et al. 2004; Schnurr et al. 2008a).

2

Indeed, if we adopt (_H_−_K_s)0=−0.11 mag (Crowther et al. 2006) for F6 and F9, mag is implied from the Espinoza et al. (2009) photometry.

3

The exponent flattens further at the highest masses, such that _L_∝_M_1.3 for ZAMS stars between 300 and 500 M⊙.

We are grateful to C. J. Evans for providing VLT MAD images of R136 and photometric zero-points in advance of publication, and E. Antokhina, J. Maíz Apellániz, J. M. Pittard and I. R. Stevens for helpful discussions. We thank C. J. Evans, N. R. Walborn and an anonymous referee for their critical reading of the manuscript. Based on observations made with ESO Telescopes at the Paranal Observatory during MAD Science Demonstration runs SD1 and SD2, plus programme IDs 69.D-0284 (ISAAC), 075.D-0577 (SINFONI), 076-D.0563 (SINFONI). Additional observations were taken with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Institute. STScI is operated by the association of Universities for Research in Astronomy, Inc., under the NASA contract NAS 5-26555. Financial support was provided to OS and RJP by the Science and Technology Facilities Council. NY and HAK gratefully acknowledge the University of Malaya and Ministry of Higher Education, Malaysia for financial support, while RH acknowledges support from the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. A visit by NY to Keele University was supported by UNESCO Fellowships Programme in Support of Programme Priorities 2008–2009.

REFERENCES

,

1987

,

Galactic Dynamics

.

Princeton Univ. Press

Princeton, NJ

et al.,

2004

,

ApJ

,

611

,

L33

,

2007

, in , eds, ASP Conf. Ser. Vol. 367,

Massive Stars in Interacting Binaries

.

Astro. Soc. Pac.

, San Francisco, p.

629

,

2008

,

Cambridge Astrophys. Ser. Vol. 45, From Luminous Hot Stars to Starburst Galaxies

.

Cambridge Univ. Press

, Cambridge

et al.,

2005

,

ApJ

,

628

,

246

et al.,

2010

, in , eds,

Proc. IAU Symp. 266

,

Star Clusters: Basic Galactic Building Blocks Throughout Time and Space

.

Cambridge Univ. Press

, Cambridge, p.

35

et al.,

2010

,

ApJ

,

715

,

L74

,

2008

, in , eds,

Proc. IAU Symp. 250

,

Massive Stars as Cosmic Engines

.

Cambridge Univ. Press

, Cambridge, p.

247

et al.,

2002

,

ApJ

,

581

,

258

et al.,

2009

,

Nat

,

462

,

624

,

2000

,

Mem. Soc. Astron. Ital.

,

71

,

693

et al.,

2008

,

ApJ

,

676

,

L29

,

1987

,

Revista Mexicana Astron. Astrofisica

,

14

,

265

et al.,

1999

,

ApJS

,

123

,

3

,

2009

,

Physics, Formation and Evolution of Rotating Stars

.

Springer-Verlag

, Berlin

,

2008

, in , eds,

Proc. IAU Symp. 250

,

Massive Stars as Cosmic Engines

.

Cambridge Univ. Press

, Cambridge, p.

119

et al.,

2002

,

ApJ

,

573

,

191

et al.,

2007

,

A&A

,

465

,

1003

et al.,

2007

,

ApJ

,

659

,

L13

,

2002

, in , eds, ASP Conf. Ser. Vol. 260,

Interacting Winds from Massive Stars

.

Astron. Soc. Pac.

, San Francisco, p.

449

et al.,

2005

,

A&A

,

432

,

985

,

2008

, in , eds,

Proc. IAU Symp. 250

,

Massive Stars as Cosmic Engines

.

Cambridge Univ. Press

, Cambridge, p.

536

,

2008

,

Diplomarbeit

, Universität Potsdam, Germany

,

1999

,

Revista Mexicana Astron. Astrofisica

,

8

,

41

et al.,

2006

,

AJ

,

131

,

1163

et al.,

2007

,

ApJ

,

666

,

1116

,

2010

, in , Bennett P. D., Morris P. W., Van Loon J. Th., eds, ASP Conf. Ser. Vol. 425, Hot and Cool: Bridging Gaps in Massive Star Evolution.

Astron. Soc. Pac.

, San Francisco, p. 45

et al.,

2002

,

AJ

,

123

,

2754

© 2010 The Authors. Journal compilation © 2010 RAS

Citations

Views

Altmetric

Metrics

Total Views 4,817

2,879 Pageviews

1,938 PDF Downloads

Since 1/1/2017

Citations

407 Web of Science

×

Email alerts

Citing articles via

• Latest

• Most Read

• Most Cited

More from Oxford Academic