Massive Stars in the Tarantula Nebula: A Rosetta Stone for Extragalactic Supergiant HII Regions

Abstract

A review of the properties of the Tarantula Nebula (30 Doradus) in the Large Magellanic Cloud is presented, primarily from the perspective of its massive star content. The proximity of the Tarantula and its accessibility to X-ray through radio observations permit it to serve as a Rosetta Stone amongst extragalactic supergiant HII regions since one can consider both its integrated characteristics and the individual properties of individual massive stars. Recent surveys of its high mass stellar content, notably the VLT FLAMES Tarantula Survey (VFTS), are reviewed, together with VLT/MUSE observations of the central ionizing region NGC 2070 and HST/STIS spectroscopy of the young dense cluster R136, provide a near complete Hertzsprung-Russell diagram of the region, and cumulative ionizing output. Several high mass binaries are highlighted, some of which have been identified from a recent X-ray survey. Brief comparisons with the stellar content of giant HII regions in the Milky Way (NGC 3372) and Small Magellanic Cloud (NGC 346) are also made, together with Green Pea galaxies and star forming knots in high-z galaxies. Finally, the prospect of studying massive stars in metal poor galaxies is evaluated.

1. Introduction

The Tarantula Nebula (alias 30 Doradus) in the Large Magellanic Cloud (LMC) is the brightest supergiant HII region in the Local Group of galaxies, and serves as a local analogue to metal-poor starburst knots in high redshift galaxies [1]. Its proximity (50 kpc), and high galactic latitude (and hence low extinction) have permitted a myriad of ground-based and space-based surveys across the electromagnetic spectrum, revealing an exceptional population of massive stars (≥$8M_{\odot }$), including the dense, young star cluster R136 that is home to some of the most massive stars known. The advent of modern highly multiplexed spectrographs coupled with large ground-based telescopes, has permitted multi-epoch optical spectroscopic surveys of the massive star population of the Tarantula Nebula for the first time. The VLT FLAMES Tarantula Survey, hereafter VFTS [2], has provided the multiplicity, rotational velocities and initial mass function of massive stars. In distant starburst regions it is not possible to resolve individual stars, such that studies rely on techniques based on their integrated properties, which themselves involve assumptions of binarity, stellar rotation and mass function.

The Tarantula Nebula is not the sole example of a supergiant H ii region within the Local Group. However, studies of the richest Milky Way star-forming regions are limited by dust extinction (e.g., Westerlund 1, NGC 3603). Counterparts in other Local Group galaxies suffer from a number of limitations, including a relatively modest stellar content (NGC 346 in the Small Magellanic Cloud, SMC) or much greater distance (NGC 604 in M 33), such that the Tarantula—whose metallicity is approximately half-solar [3]—serves as the only credible Rosetta Stone for rich extragalactic star-forming regions. This review will provide a brief overview of the structural properties of the Tarantula Nebula, but will largely focus on its massive star content, drawn from results from VFTS and other optical spectroscopic surveys, supplemented by the recent deep Chandra X-ray survey ‘Tarantula-Revealed by X-rays’ (T-ReX). Finally, comparisons with local and high-redshift star-forming regions will be provided to put the properties of the Tarantula into a broader context. Indeed, the nebular properties of the central NGC 2070 region of the Tarantula are strikingly similar to Green Pea galaxies, exhibiting intense [O iii] $\lambda$4959, 5007 emission, and its star-formation rate is comparable to intense star-forming clumps at high redshift.

2. Tarantula Nebula

The Tarantula Nebula is the most striking star-forming region in the LMC, whether viewed in the far ultraviolet (hot, luminous stars), H$\alpha$ (ionized gas) or mid-IR (warm dust), extending over several hundred parsec, owing to a massive stellar content producing an ionizing output which is a thousand-fold higher than the Orion Nebula [4].

NGC 2070 is the dominant ionized region within the Tarantula Nebula, powered by a large number of hot luminous stars from R136 at its heart plus many more within its vicinity. NGC 2060, located 6 arcmin (90 pc) to the southwest, is host to a more modest population of OB stars plus an X-ray pulsar PSR J0537-6910 and its supernova remnant (SNR) N157B. Hodge 301 is located 3 arcmin (45 pc) to the north west of R136, but does not possess significant nebulosity since previous supernovae are likely to have cleared this region of gas, and its stellar population (B-type stars and red supergiants) does not possess a significant Lyman continuum output. NGC 2070 is often referred to as a cluster in the literature but it extends over tens of parsecs whereas genuine star clusters are an order of magnitude smaller, such that the only rich star clusters within the Tarantula are R136, with an age of 1–2 Myr [5] and Hodge 301 with an age of 20–30 Myr [6] with a few additional lower mass young, compact clusters (e.g., TLD1, SL 639).

Table 1. Physical scales within the Tarantula Nebula, adapted from Walborn [7]. Ionizing outputs, N(LyC) are obtained from the present work.

Region	Angular	Physical	N(LyC)	Content	Reference
	Radius (${}^{″}$)	Radius (pc)	(${10}^{51}$ ph s${}^{-1}$)
R136a	0.8	0.2	2	R136a1 (WN5h), R136a2 (WN5h)	[9]
R136	4	1.0	4	R136b (O4 If/WN8), R136c (WN5h+)	[5]
NGC 2070	80	20.	9	R140a (WC4 + WN6+), Mk34 (WN5h + WN5h)	[10,11]
Tarantula	600	150.	12	Hodge 301, PSR J0537-6910 (pulsar), N157B (SNR)	[12]

compares various regions within the Tarantula Nebula, adapted from a previous review by Walborn [7]. Although the focus of the present review is on spectroscopic results for massive stars, Sabbi et al. [8] have undertaken a deep Hubble Space Telescope (HST) multi-colour photometric survey, known as the Hubble Tarantula Treasury Project (HTTP) which permits lower mass stars in the Tarantula Nebula to be studied.

A number of complementary studies of the star formation history of the Tarantula Nebula have been carried out, exploiting pre-main sequence low mass stars [13,14] and massive stars [15]. Significant star formation commenced ∼25 Myr ago, as witnessed by Hodge 301, and reached a peak several Myr ago, with the young massive star cluster R136 at its heart no more than 2 Myr old. It is apparent that star formation within the Tarantula Nebula has not been limited to specific parsec-scale star clusters, such as Hodge 301 or R136, but has been distributed across the entire region, akin to a super OB association. Wright et al. [16] have established from proper motion observations that star formation in the far smaller Milky Way Cygnus OB2 region did not originate in a star cluster, but involved individual sub-regions in virial equilibrium. Indeed, median ages of massive stars show little radial dependence on their projected distance from R136, with very massive stars (≥$100M_{\odot }$) identified throughout the region [15]. Infrared and radio observations of the Tarantula reveal ongoing regions of massive star formation, to which the reader is referred to the review by Walborn [17] and a more recent study of the brightest embedded sources based on Spitzer/IRAC imaging [18]. Atacama Large Millimetre Array (ALMA) has obtained high resolution observations of parsec-scale clumps within the Tarantula [19], with the rate of star formation in the Tarantula anticipated to decline in the future.

3. Massive Star Content

The integrated nebular properties of the Tarantula Nebula and other selected giant H ii regions in the Local Group is presented in

Table 2. Integrated nebular properties of nearby giant HII regions, adapted from Kennicutt [20] for an assumed O7V Lyman continuum output of 10${}^{49}$ ph s${}^{-1}$.

Region (Galaxy)	Diameter	L(H$\mathit{\alpha }$)	N(O7V)
	pc	10${}^{39}$ erg s${}^{-1}$
NGC 3372 (Milky Way)	200	0.8	0045
NGC 346 (SMC)	220	0.8	0045
NGC 3603 (Milky Way)	100	1.5	0110
NGC 604 (M33)	400	4.5	0320
Tarantula (LMC)	370	150	1100

. The ionizing output of the Tarantula Nebula corresponds to the equivalent of over a thousand O7V stars, each with 10${}^{49}$ ph s${}^{-1}$. In reality, the Tarantula hosts somewhat fewer O-type stars since the most extreme examples—early O stars and luminous Wolf-Rayet (WR) stars—each produce an order of magnitude more Lyman continuum photons. Nevertheless, this population represents an order of magnitude more O and WR stars than any Milky Way or SMC giant HII regions, and is not likely to be improved upon until extremely large telescopes are capable of resolving the massive stellar content of more extreme giant HII regions, such as NGC 5461, 5462 and 5471 in M 101 [1].

Historically, there have been several photometric and spectroscopic surveys of early-type stars in the Tarantula Nebula, each of which have employed contemporary (Galactic) temperature calibrations to produce Hertzsprung-Russell diagrams. Parker [21,22] obtained the first extensive study of the entire region. Subsequently, Massey and Hunter [23] obtained high-spatial resolution HST spectroscopy of early-type stars in the central, crowded R136 region, revealing an exceptional population of very early O stars, whilst Melnick and colleagues [24,25] obtained spectroscopy for a large number of early-type stars in the NGC 2070 region. The advent of efficient multi-object spectrographs on 8–10 m telescopes, has permitted the most comprehensive optical spectroscopic survey of massive stars in the Tarantula Nebula to date through the VLT FLAMES Tarantula Survey (VFTS) [2]. Multi-object, multi-epoch spectroscopy of ∼800 massive stars across the entire region, for which detailed spectroscopic analyses have been undertaken for over 500 O and early B stars, such that temperature calibrations are no longer necessary. Although this represents the most extensive of early-type stars in a single star-forming region to date, this survey is incomplete owing the fibre-placement limitations, sampling ∼70% of massive stars exterior to the dense R136 cluster from comparison with photometric surveys [12].

Two additional surveys have recently provided complete optical spectroscopic observations of all bright sources within the central crowded region of the Tarantula: (a) the central 4 arcsec (1 pc) of the R136 cluster exploiting the high spatial resolution of HST/STIS [5]; (b) the central 2 × 2 arcmin (30 × 30 pc) region of NGC 2070 using the VLT/MUSE integral field spectrograph [10]. Although these lack the multi-epoch capabilities of VFTS, they are complementary since the richest stellar populations of the Tarantula are found within NGC 2070/R136, and provide extended wavelength coverage (yellow and red for MUSE, ultraviolet for STIS), albeit at reduced spectral resolution. A summary of these three contemporary surveys is provided in

Table 3. Summary of stellar content of the Tarantula Nebula from recent spectroscopic surveys (excluding sources in common, although including individual components within SB2 binaries).

Telescope/inst	Target	N(O-type)	N(B-type)	N(WR)	N(Of/WN)	N(A+)	Reference
VLT/FLAMES	30 Dor	369	436	9	6	35	[2]
HST/STIS	R136	57	..	3	2	..	[5]
VLT/MUSE	NGC 2070	115	79	..	..	1	[10]
Other	30 Dor	29	8	16	..	5	[22,23,24,26]
Total	30 Dor	570	523	28	8	41

, together with literature results. In total, approximately 1100 early-type massive stars have been spectroscopically observed. Our discussion of the massive star content of the Tarantula will largely draw upon results from VFTS, but include others where appropriate. It should be noted that analyses based on spectroscopic fibre-fed observations of early-type stars in regions of strong, highly variable nebulosity, is inherently problematic owing to the lack of local sky subtraction. Such issues do no affect long-slit STIS or integral field MUSE observations.

VFTS and other surveys have revealed that the Tarantula hosts extreme examples of stellar exotica, including the most massive stars known [9], a Luminous Blue Variable, R143 [27], a very massive runaway [28], the fastest rotating stars [29], a massive overcontact binary [30], and a supernova remnant N157B [31], with SN 1987A located ∼300 pc to the south west (Figure 1).

R136a, the central cluster, merits special consideration since it was considered by some to be a supermassive star as recently as the early 1980s [32]. Speckle interferometry resolved R136a into multiple sources [33], and it was subsequently established as a compact star cluster [34]. Massey and Hunter [23] established that dozens of the brightest sources within the central parsec were hot, early O stars. Spectroscopic studies of the brightest components R136a1, a2, a3, with nitrogen-sequence Wolf-Rayet spectral types, indicated masses of ∼$100M_{\odot }$ [23,35]. They established that these relatively weak-lined WN stars are luminous main-sequence stars close to their Eddington limits, rather than classical Wolf-Rayet stars. Subsequent analyses of the WN stars in R136 indicated significantly higher masses of 150–300 $M_{\odot }$ [9] as a result of increased spectroscopic luminosities, owing to higher stellar temperatures and IR photometry less affected by dust extinction. Indirectly inferred masses of massive stars are notoriously imprecise, and if binarity were established for individual stars their inferred luminosities and masses would be reduced. To date, faint companions to members of R136a have been detected with extreme adaptive optics imaging [36]. Melnick 34 is spectroscopically similar to the WN5-stars in R136a, and has recently been shown to be a colliding wind binary system comprising two WN5 components, with a total mass exceeding 250 $M_{\odot }$ [37]. Figure 2 shows that Melnick 34 is an order of magnitude brighter than R136a in X-rays, indicating that there are no colliding wind binaries comparable to Melnick 34 within R136a [38].

Overall, the Tarantula hosts a remarkable number of ∼50 early-type stars with bolometric luminosities exceeding ${10}^6{L}_{\odot }$. For reference, the Milky Way’s Carina Nebula (NGC 3372) hosts ≈5 massive stars with such extreme properties [39]. As such, the Tarantula Nebula represents our best opportunity to study the highest mass stars known, both individually and collectively. Schneider et al. [40] analysed VFTS spectroscopic results to establish an excess of massive stars with respect to a standard Salpeter Initial Mass Function (|MF), indicating 1/3 more stars with ≥$30M_{\odot }$ in 30 Doradus compared to expectations from a standard IMF. Finally, although we focus primarily on high mass early-type stars in the Tarantula Nebula, it also hosts red supergiants (RSG). Since RSG are the evolved descendants of moderately massive stars, and star formation in the Tarantula has peaked relatively recently, of order ∼10 RSG are known, most of which are associated with mature star clusters Hodge 301 and SL 639 [6].

4. Physical Properties

The determination of physical properties ($T_{\mathrm{eff}},logg$) of individual early-type stars ideally requires high S/N (≥50) intermediate resolution spectroscopy of suitable diagnostics, usually He i–ii lines (Si ii–iv) for temperatures of O-type (B-type) stars, plus Balmer lines for surface gravities, plus grids of model atmospheres obtained with modern codes, such as FASTWIND [41] for O stars and blue supergiants, CMFGEN [42] or PoWR [43] for emission-line stars, or TLUSTY [44] for low luminosity B stars. Analysis of late-type supergiants requires model atmosphere codes in which molecular opacities have been incorporated, such as MARCS [45]. Hot O stars require alternative temperature diagnostics to helium, with nitrogen commonly used since the blue visual spectrum of O stars includes lines of N iii–v. Wolf-Rayet stars are especially problematic since photospheres are masked by the dense wind, such that gravities cannot be directly measured and temperatures usually refer to deep layers, with an optical depth of $\tau \sim$ 10–20, rather than the effective temperature at $\tau =2/3$. If stellar distances are uncertain, comparisons with evolutionary models can be made using the so-called spectroscopic Hertzsprung-Russell (sHR) diagram [46], involving temperature and $\mathcal{L}=T_{\mathrm{eff}}^4/g$, the inverse of the flux-weighted gravity, where g is the surface gravity.

The determination of stellar luminosities requires comparisons between synthetic spectral energy distributions and photometry, taking account of interstellar extinctions and distance moduli, 18.5 mag in the case of the Large Magellanic Cloud. Visual extinctions of early-type stars in the Tarantula Nebula are usually modest, although near-IR photometric comparisons usually lead to more robust luminosities since typical dust extinctions are 0.1–0.2 mag in the K-band, versus 1–2 mag in the V-band, and the lack of sensitivity of K-band extinctions to any variations in the overall extinction law. Luminosities of RSG can also be reliably estimated by integrating observed spectral energy distributions from visual to mid-IR wavelengths [47,48]. Historically, stellar estimates of masses and ages from evolutionary models involved by-eye comparisons between their position on a conventional Hertszprung-Russell (HR) diagram and theoretical isochrones. However, additional physical information is often available, such as helium abundance or projected rotational velocities. Tools now exist which additionally take such information into account for the calculation of stellar ages and initial masses such as BONNSAI [49]. Significant discrepancies exist between current mass estimates from spectroscopic ($logg$) and evolutionary approaches for a subset of VFTS O dwarfs [50].

Figure 3 presents the HR diagram of the Tarantula Nebula, comprising single star results from VFTS [50,51,52,53,54], VLT/MUSE [11], HST/STIS [55] and literature results for other stars within 160 parsec of R136, including Wolf-Rayet stars [56]. Results for binary systems have been incorporated, primarily drawn from [57] for VFTS B-type binaries, Tarantula Massive Binary Monitoring (TMBM) for VFTS O-type binaries [58] and recent literature for WR stars [37,59]. Evolutionary tracks for non-rotating, LMC metallicity massive stars up to the onset of He-burning have been included for reference [60,61]. Over 1170 massive stars have been included, revealing a well populated main sequence population up to ∼$200M_{\odot }$, plus classical Wolf-Rayet stars to the left of the main sequence, and evolved blue supergiants up to $log(L/L_{\odot })\sim 6$, and cool supergiants, up to $log(L/L_{\odot })\sim$5.3 [48]. The addition of all luminous early-type stars from R136 and NGC 2070 fills in the extreme upper main sequence which is somewhat under populated from VFTS alone [40]. The overwhelming majority of the older massive stellar population—i.e., evolved stars with masses below $30M_{\odot }$—are spatially exterior to NGC 2070 (open symbols), although NGC 2070 is host to one luminous M supergiant, Melnick 9. Conversely, beyond NGC 2070, the main-sequence population at the highest stellar masses is relatively underpopulated, albeit with several WN5h stars (R146, R147) and early O stars (VFTS 16, BI 253) located 95 ± 25 parsec from R136.

5. Binaries, Rotation and Runaways

Until recently, the significance of close binary evolution for massive stars was not fully recognised, in spite of a few binary “champions” [62]. The high frequency of close binaries amongst O stars in young Galactic clusters obtained from radial velocity monitoring established that only a minority of massive stars follow single stellar evolution [63]. In contrast with the majority of previous spectroscopic surveys of early-type stars, VFTS comprised multiple epochs, such that [64] were able to establish that 53% of O stars in the Tarantula Nebula inhabit a binary system with a period below 1500 days, such that binary interaction will occur. In total, 18% of O-type binaries, those with very short periods, are anticipated to merge with a companion, 27% will be stripped of their envelopes (primaries, mass donors) and 8% are predicted to be spun up (secondaries, mass gainers), as summarised in Figure 4 together with counterparts in Milky Way clusters. Broadly similar results have been obtained for VFTS B-type stars [65].

The inferred rate of envelope stripping and spin-up in the Tarantula is rather lower than [63] obtained for O stars in Milky Way clusters (Figure 4), but it is probable that the true incidence is rather higher since a subset of the current O star sample is likely to have already undergone binary evolution. Although VFTS has established the binary frequency amongst massive stars in the Tarantula, binary orbits require follow-up surveys, notably the TMBM survey [58,66].

Consequences of close binary evolution include mass gaining secondaries being spun-up, and a subset of secondaries possessing high space velocities as a result of the disruption of the binary following the core-collapse supernova (ccSN) of the original primary. Other “observables" are more challenging, including the identification of stripped primaries in close binaries which will usually be masked at visual wavelengths by mass gaining secondaries [67] which should be more common in older star clusters such as Hodge 301.

Reliable measurements of projected rotation rates, $v_{e}sini$, are often problematic because strong hydrogen and helium lines are predominantly affected by pressure broadening, with an additional contribution from“macroturbulence”. A Fourier Transform approach applied to metallic lines offers the most robust results, albeit requiring high resolution, high S/N spectroscopy of suitable diagnostics [68]. The lack of a spectral features originating in the hydrostatic layers of Wolf-Rayet stars prevents a direct determination of their rotational velocities. Figure 5 compares projected rotational velocities for a large sample of VFTS O and B-type stars. Typical rotational velocities of single O stars are modest, with $v_{e}sini\sim$ 100 km s${}^{-1}$, albeit with 10% exceeding 300 km s${}^{-1}$ [69], of which some examples are rotating close to their critical rates [29]. This high velocity tail is suspected of being spun-up mass gainers in former close binaries. The lack of fast rotators amongst O stars in VFTS binary systems supports this interpretation [70]. Figure 5 illustrates that rotational velocities of VFTS B-type stars are higher, with $v_{e}sini\sim$ 200 km s${}^{-1}$ on average, and 20% exceeding 300 km s${}^{-1}$ [71].

Of particular interest is the rotational velocity distribution of massive stars within the young R136 star cluster, whose severe crowding prevented inclusion in VFTS. Bestenlehner et al. [55] utilised HST/STIS spectroscopy to reveal 150 km s${}^{-1}$ on average for a sample of 55 massive stars within the central parsec of R136, with no examples exceeding 250 km s${}^{-1}$, although the low S/N of these datasets prevented distinguishing between rotational and macroturbulence, adding to the interpretation that rapid rotators originate from close binary evolution. Wolff et al. [72] obtained somewhat higher rotational velocities for OB stars in the periphery of R136, where contamination from the field population is significant.

The Tarantula hosts a number of candidate early-type runaway stars from their measured (radial and/or tangential) velocities with respect to the average for their environment, although radial velocity outliers may be unresolved binaries. Runaways can originate either from disrupted secondaries following the core collapse of primaries in close binaries, or following the dynamical ejection of stars from young star-forming regions. Platais et al. [73] have investigated high proper motion stars in the Tarantula from HST imaging obtained 3 years apart, revealing a number of potential stars ejected from R136, while Lennon et al. [28] have exploited Gaia DR2 proper motions to conclude that VFTS 16 (O2 III) [74] was likely ejected from R136 during its formation 1–2 Myr ago. Renzo et al. [75] discuss the origin of the candidate ‘walkaway’ very massive star VFTS 682 (WN5h).

The Tarantula hosts several notable massive binary systems, whose physical and orbital properties have been obtained from spectroscopic monitoring, with searches for massive binaries also greatly benefitting from the recent Chandra T-ReX survey (PI Leisa Townsley) which monitored the Tarantula in X-rays for almost 2 years with a total integration time of 2 Ms. Single hot, luminous stars tend to produce (thermal) X-rays due to shocks in the winds, but these are generally soft X-ray emitters with $L_X/L_{\mathrm{Bol}}\sim {10}^{-7}$ [76]. Massive stars in binary systems may lead to excess X-ray emission arising from wind-wind collisions, usually relatively hard, providing the separations are not too small (low wind velocities) or too large (low wind densities) [77]. A close binary comprising an early-type star and compact remnant (neutron star or black hole) will be extremely X-ray bright if the accretion disk of the remnant is being fed by the wind of the massive star or via Roche Lobe overflow.

A number of eclipsing binaries in the proximity of R136 have been identified [78], including # 38 from [34] comprising an O3 V + O6 V in a circular 3.4 day orbit, with component masses 57 and 23 $M_{\odot }$. This represented the first robust stellar mass determination for an O3 star in the LMC. VFTS revealed a large number of binaries within the Tarantula, many of which have been followed-up with TMBM. Most notably R139 has been established as an eccentric system comprising a pair of mid O supergiants in a 154 day orbit [79] with lower limits of ∼$66+78M_{\odot }$ for individual component, recently revised downward to $54+69M_{\odot }$ [58]. R139 is amongst the brightest X-ray sources in the Tarantula in T-ReX with $L_{\mathrm{X},\mathrm{corr}}\sim 5\times {10}^{33}$ erg s${}^{-1}$ and an enhanced $L_{\mathrm{X},\mathrm{corr}}/L_{\mathrm{Bol}}\sim 9\times {10}^{-7}$ based on TMBM bolometric luminosities [58].

The most remarkable X-ray source in the Tarantula is VFTS 399 with $L_{\mathrm{X},\mathrm{corr}}\sim 5\times {10}^{34}$ erg s${}^{-1}$ despite being associated with a low luminosity O9 giant, implying $L_{\mathrm{X},\mathrm{corr}}/L_{\mathrm{Bol}}\sim 2\times {10}^{-4}$. Clark et al. [80] conclude that VFTS 399 is a high-mass X-ray binary hosting a neutron star remnant, with the O giant known to be a rapid rotator ($v_{e}sini=324$ km s${}^{-1}$, according to [53]), as one would expect for a mass gainer in a close binary system.

Two point sources in the Tarantula are even brighter in X-rays than VFTS 399. Of these, R140a is a compact group of stars including two WR stars, so likely hosts one or more colliding wind binaries, while X-ray variability for Melnick 34 reveals a 155 day period, peaking at $L_{\mathrm{X},\mathrm{corr}}=3.2\times {10}^{35}$ erg s${}^{-1}$, exceeding $\eta$ Carina at X-ray maximum [38]. Melnick 34 has been confirmed to be a colliding wind binary system in a 155 day eccentric orbit, with minimum masses of 60–65 $M_{\odot }$ for each of the WN5h components [37]. Individual masses of 130–140 $M_{\odot }$ are favoured from spectroscopic analysis, such that Melnick 34 is likely to be the most massive binary known to date, with $L_{\mathrm{X},\mathrm{corr}}/L_{\mathrm{bol}}=1.7\times {10}^{-5}$ at X-ray maximum. This arises from the collision of dense, fast moving winds at a minimum separation of ∼1.2 AU or 13 stellar radii for an assumed orbital inclination of $i={50}^{\circ }$. Almeida et al. [30] identify the overcontact binary VFTS 352 as a prototype of systems which, at low metallicity, are plausible black hole-black hole merger progenitors.

6. Wind Properties

The underlying theory responsible for outflows by hot luminous stars has been known for several decades [81], with a metallicity dependence primarily arising from the variation in iron-peak elemental abundances [82,83]. Individual wind properties of O stars or blue supergiants usually rely on spectroscopic fits to H$\alpha$, as parameterised by the wind strength parameter, Q [41], from which the mass-loss rate requires knowledge of the physical radius (from $T_{\mathrm{eff}},logL$) and measured or adopted wind velocity. Wind velocities of OB stars cannot be measured from optical spectroscopy, so usually spectral type calibrations are adopted based on measured velocities from UV P Cygni profiles of C iv, N v or Si iv [84]. Until recently, high S/N, high quality UV spectroscopy of OB stars in the Large Magellanic Cloud has been in short supply, but the situation has improved via the HST Large Program METAL (GO 14675, PI Julia Roman-Duval) [85] and upcoming ULLYSES initiative 1.

Specifically for the Tarantula, low resolution UV spectroscopy of the R136 star cluster has added a significant number of wind measurements for early O stars [5]. Wind velocities exceed 3000 km s${}^{-1}$ at the earliest subtypes (O2–3), reducing to ∼1500 km s${}^{-1}$ for late O-types. Mass-loss rates of emission line stars typically rely on an alternative wind scaling relation, namely the transformed radius, $R_t$ [86], which also requires knowledge of wind velocities, although these can be estimated from optical spectroscopy for Wolf-Rayet stars. An added complication arises because radiatively-driven winds are known to be inherently unstable, leading to clumped winds. Wolf-Rayet winds have been known to be clumped for 30 years [87,88], but the degree of clumping for O stars via the H$\alpha$ diagnostic remains unclear, with conflicting results obtained from UV resonance lines [89] unless the wind comprises a mixture of optically thin and thick clumps within a much lower density inter-clump medium [90].

Figure 6 presents unclumped mass-loss rates of O, Of/WN and Wolf-Rayet stars in the Tarantula Nebula, obtained from VFTS [50,51,53], HST/STIS [55] and other literature results. Uncertainties have been included wherever possible. Since the primary wind diagnostic in the majority of instances presented here is H$\alpha$, it is apparent that uncertainties are large for those stars with weak stellar winds. In addition, mass-loss rates for Of/WN and Wolf-Rayet stars are anticipated to be reduced significantly owing to wind clumping, as indicated with downward arrows. If volume filling factors are ∼10%, mass-loss rates will be reduced by a factor of $\sqrt{10}$.

Theoretical mass-loss rates [83] for zero-age main sequence stars at LMC composition [60,61] are included in Figure 6. At face value it would appear that the theoretical mass-loss rates of LMC O stars are supported by theory. However, the following should be borne in mind. It is not clear how significantly wind clumping affects the inferred mass-loss rates of normal O stars, although H$\alpha$ results for supergiants are likely to be sensitive to wind clumping. In addition, the vast majority of mass-loss rates of VFTS O stars shown here have been inferred by adopting wind velocities from an assumed scaling relation involving escape velocities [91,92], which are themselves dependent upon spectroscopic gravities. Exceptions are HST/STIS results for early-type stars in R136 which are based on measured UV wind velocities, and span dwarfs, giants, supergiants and main sequence WN stars. In order to verify predictions for lower luminosity ($logL/L_{\odot }<5.5$) O stars, more sensitive diagnostics would need to be employed, such as UV P Cygni lines, providing complications such as porosity are accounted for [90].

It is clear from Figure 6 that rates for the highest luminosity main-sequence Of/WN and WN stars significantly exceed theoretical predictions. This discrepancy is partially addressed through wind clumping, but very massive stars close to their Eddington limits are observed to exhibit enhanced mass-loss rates which are not taken into account in standard theoretical predictions [93,94]. Unsurprisingly, classical Wolf-Rayet stars with $log(L/L_{\odot })=5.5-6$ possess the strongest winds amongst early-type stars in 30 Doradus, with clumping-corrected wind densities an order of magnitude higher than O stars with similar luminosities. It is well known that the wind momenta of WR stars, $\dot{M}{v}_{\infty }$, exceeds the momentum provided by their radiation field, $L/c$, owing to multiple photon absorption and re-emission within their optically thick winds, permitting $\dot{M}{v}_{\infty }/(L/c)>1$ [95].

7. Fate of Massive Stars in the Tarantula

The conventional picture of massive star evolution close to solar metallicity is that those with initial masses of 8–25 $M_{\odot }$ will end their lives as red supergiants (RSG), undergo a H-rich core collapse supernova, leaving behind a neutron star remnant, while higher mass counterparts will either circumvent the RSG phase or subsequently proceed to a Wolf-Rayet stage prior to undergoing core collapse, leading to a H-deficient supernova (neutron star remnant) or faint/failed supernova (black hole remnant) [96]. If one considers the global WR vs. RSG population in the LMC, the lower boundary to the luminosity of WR stars is $log(L/L_{\odot })$ = 5.3, while the upper luminosity of RSG is $log(L/L_{\odot })$ = 5.5 [48], supporting a transition from RSG to WR for higher mass progenitors at ∼25 $M_{\odot }$.

Close binary evolution severely complicates this scenario, since primaries below 25 $M_{\odot }$ can be stripped of their hydrogen envelope, leading to a type IIb or Ib/c instead of a H-rich supernova, while secondaries will be rejuvenated, spun-up, with the potential for a core-collapse supernova for secondaries whose initial masses fall below 8 $M_{\odot }$. To date, there are no unambiguous cases of pre-supernova close binaries in the Tarantula hosting stripped (Wolf-Rayet or helium) stars, although it has been suggested that the WN3 binary BAT99-49 elsewhere in the LMC is the product of close binary evolution [59]. Rapid rotation of the bright O giant component of the high mass X-ray binary VFTS 399 is consistent with this evolutionary scenario. The absence of low luminosity Wolf-Rayet stars in the Tarantula does not exclude the binary channel since low-mass stripped stars would be unlikely to exhibit a Wolf-Rayet spectral appearance [67].

Initially very close binaries may merge on the main sequence, prior to following a relatively conventional evolution, albeit with unusually high rotation rates, which would lead to increased luminosities and potentially evolve blueward off the main sequence [60,61]. Extremely rapid rotation in some VFTS OB stars favours close binary evolution or stellar mergers. Very massive stars in the Tarantula up to ∼$300M_{\odot }$ are expected to lead to 30–50 $M_{\odot }$ CO cores and black hole fates, unlikely to produce any associated supernova [97], such that a subset of binary VMS are plausible progenitors of LIGO black hole binary mergers, although their exact fate crucially depends on their mass-loss properties, which remain uncertain.

8. Integrated Properties and Comparison With Star-Forming Regions, Near And Far

The Tarantula Nebula would subtend little more than one arcsec if it were located at a distance of 50 Mpc, and so provides us with a unique opportunity to compare the individual spatially-resolved properties of an intensively star-forming region and its aggregate characteristics. Using integrated H$\alpha$ observations of the Tarantula Nebula [98], an age of ∼3.5 Myr would be inferred from a comparison between the inferred H$\alpha$ equivalent width of 1100 Å and population synthesis models for a coeval population at LMC metallicity [12]. This is in reasonable agreement with the typical age of massive stars, albeit failing to reflect the complexity in its star-formation history (recall Figure 3).

An analysis of the integrated UV spectrum of NGC 2070 supports a young (≤3 Myr) starburst episode [99], while the high spatial resolution of HST/STIS has permitted a comparison between the individual and integrated UV spectroscopic appearance of the central R136 cluster [5]. Very massive stars contributed a significant fraction of its far UV continuum flux, and completely dominate the strong, broad He ii $\lambda$1640 emission. The integrated UV spectroscopic appearance of R136 closely resembles some star clusters in star-forming galaxies at Mpc distances, such as NGC5253-5 [100], suggesting the presence of VMS in these other young massive clusters. From comparison with the predictions of standard population synthesis models, both Starburst99 [101] and BPASS [102] models fail to predict any significant emission prior to the conventional Wolf-Rayet phase (Figure 7), owing to the use of inadequate wind theory for VMS [93,94]. Consequently, neither Starburst99 nor BPASS accounts for the powerful winds of very massive stars in R136, and the adopted mass function follows a Salpeter slope, rather than the top heavy IMF identified by [40] for the Tarantula region as a whole.

An estimate of the cumulative ionizing and mechanical feedback from massive stars within the Tarantula has revealed a major contribution from VMS towards the collective ionizing output and a dominant role from WR stars to the mechanical feedback [12]. However this analysis relied on calibrations and estimates of spectral types for a significant subset of the massive star content, so we are able to provide updates from recent spectroscopic observations (e.g., VLT/MUSE, HST/STIS) and analyses. We present the updated cumulative ionizing output from 1170 massive stars in the Tarantula Nebula in Figure 8, indicating a total Lyman continuum ionizing output of 1.2 $\times {10}^{52}$ ph s${}^{-1}$ within 150 pc of R136a. A quarter of the total ionizing radiation originates from the R136a cluster, while members of NGC 2070 produce three quarters of the global feedback. Ten systems alone, listed in

Table 4. Top ten stellar systems contributing to the Lyman continuum output of the Tarantula Nebula, comprising very massive early O stars and WN5 stars, and classical Wolf-Rayet stars, updated from [12].

Star (Alias)	Sp Type	$log\mathit{L}/{\mathit{L}}_{\odot }$	log N(LyC)	Reference
R136a1 (BAT99-108)	WN5h	6.8	50.6	[55]
Mk 34 (BAT99-116)	WN5h + WN5h	6.4 + 6.4	50.6	[37]
R136a2 (BAT99-109)	WN5h	6.7	50.5	[55]
R144 (BAT99-118)	WN5-6 + WN6-7	6.7	50.5	[56,64]
R136a3 (BAT99-106)	WN5h	6.6	50.5	[55]
Mk 49 (BAT99-98)	WN6(h)	6.7	50.5	[56]
R145 (BAT99-119)	WN6h + O3.5If/WN7	6.3+6.3	50.4	[59]
Mk 42 (BAT99-105)	O2 If	6.6	50.4	[51]
VFTS 682	WN5h	6.5	50.4	[51]
R136c (BAT99-112)	WN5h +?	6.6	50.4	[51]

, collectively contribute a quarter of the ionizing budget of the Tarantula Nebula, comprising main sequence and classical WR stars, plus early O supergiants. Indeed, half of the global ionizing output originates from 40 early O stars, main sequence WN stars and classical WR stars, with 1130 stars contributing the remaining 50%. The collective bolometric luminosity of these stars is ${10}^{8.4}{L}_{\odot }$, of which the 40 UV-bright stars contribute ${10}^{8.0}{L}_{\odot }$.

Recalling Table 2, the highest mass stars and evolved high mass stars in other giant H ii regions in the Local Group dominate their radiative and mechanical feedback, including NGC 3372 (Carina Nebula) in the Milky Way, N206 in the LMC and NGC 346 in the SMC. By way of example, Smith [39] established that only a handful of early O-type stars and H-rich WN stars contribute the majority of the Lyman continuum flux of the Carina Nebula, while $\eta$ Car, four Wolf-Rayet stars and two early O supergiants completely dominate the stellar mechanical luminosity. The central ionizing cluster of the Galactic NGC 3603 star-forming region is host to a stellar content analogous to R136a, including a number of early O stars, nitrogen-sequence Wolf-Rayet stars [103,104]. Weak main-sequence wind properties of metal-poor massive stars conspire to even fewer massive stars (HD 5980, Sk 80) dominating the cumulative stellar feedback in NGC 346. Similar conclusions were reached by Ramachandran et al. [105] for the supergiant shell in the wing of the SMC.

Although the Tarantula Nebula is the most extreme giant H ii region in the Local Group, how does it rank against star-forming regions of galaxies in the near universe or knots at high redshift? Figure 9 compares the star-formation rate versus size of regions spanning z = 0 to 3.4, adapted from [106], indicating that the Tarantula (red square) is forming stars more vigorously than typical low-redshift counterparts, resembling some star-forming regions at high redshift. Indeed, 80% of the cumulative ionizing radiation originates from NGC 2070, such that this region corresponds closely with typical clumps in the lensed galaxy SDSS J1110 + 6459 at z = 2.5 (green circles).

In addition to its unusually high star formation rate, the Tarantula also possesses high ionization parameter nebular properties with respect to star-forming galaxies in the local universe. Figure 10 presents a Baldwin, Philipps and Terlevich [107] (BPT) diagnostic diagram of Sloan Digital Sky Survey (SDSS) galaxies, in which the Tarantula (red square) has been indicated, along with Green Pea galaxies from Micheva et al. [108] which are low-metallicity, intensively star-forming galaxies exhibiting unusually strong [O iii] $\lambda$5007 emission. Steidel et al. [109] showed that z = 2–3 star forming galaxies share similar extreme nebular properties, and a subset of Green Pea galaxies have been established as Lyman continuum leakers [110,111]. Focusing again on NGC 2070, this sits amongst the extreme Green Pea galaxies in the BPT diagram.

9. Summary and Outlook

Recent comprehensive spectroscopic and imaging surveys have revealed that the Tarantula Nebula hosts the most exceptional massive star population within the Local Group of galaxies, including the most massive stars identified to date, the fastest rotating early-type stars, and the X-ray brightest colliding wind system. In particular, the VFTS survey has revealed an excess of massive stars with respect to a Salpeter IMF [40] and added support from previous results for the importance of close binary evolution in the evolution of massive stars [64]. As such, the Tarantula Nebula is the closest analogue to the Hubble Deep Field for the community interested in the evolution of massive stars since its richness provides us with a huge breadth of extreme stars.

The integrated appearance of R136 resembles young extragalactic star clusters, while the integrated nebular properties of NGC 2070 is analogous to extreme Green Pea galaxies at low redshift and star-forming knots in high-redshift galaxies. Typical metallicities of Green Pea galaxies are lower than the LMC, as measured from oxygen nebular lines, while the oxygen content of high-z star forming galaxies tends to be similar to those of the Magellanic Clouds. However, $\alpha$/Fe abundances of high redshift galaxies are likely to be higher than young populations in the Milky Way, such that winds from OB stars at high redshift are anticipated to be weaker than in the LMC or SMC [112].

The LMC metallicity is only a factor of two below that of the Solar neighbourhood [3], so ideally we would like to supplement the extensive survey of the Tarantula Nebula with counterparts at significantly lower metallicity. The SMC (1/5 solar) represents our best opportunity to study the formation and evolution of massive stars at a metallicity significantly below that of the LMC. Alas, it does not host as rich a massive star-forming region as the Tarantula, but cumulatively does host a substantial number of O stars so is key towards our improved understanding of massive stars at low metallicity, especially as it can be studied in exquisite detail with current instrumentation.

Since we need to look beyond the SMC to study metal-poor counterparts to the Tarantula Nebula, Figure 11 compares the star-formation rate, metallicity, and distance modulus of Local Group dwarf galaxies. Rates of star formation in metal-poor (≤20% of solar oxygen content) galaxies are significantly lower than the Magellanic Clouds, so there are no rich metal-poor massive star populations elsewhere in the Local Group, and those few that are present are much more distant. The absence of nearby metal-poor counterparts to the Tarantula implies that we currently have to rely on the interpretation of integrated populations in order to understand massive stellar evolution at low metallicity, notably extremely metal-poor dwarf star-forming galaxies I Zw 18 and SBS-0335.

In order to test predictions of the metallicity dependence of massive star winds, it is necessary to measure mass-loss properties across a wide range of metallicities. Although theory has been qualitatively supported from the observed wind properties of Milky Way, LMC and SMC early-type stars [82], some issues remain, including weak winds in low luminosity OB stars. Our only opportunity to study individual massive stars below 1/10th of the solar oxygen content is to observe O stars in Sextans A and B with 7% solar [114] or the Sagittarius Dwarf Irregular Galaxy (SgrDIR) with 5% solar [115]. Stellar winds from early-type stars at such low metallicities are anticipated to be much weaker than in metal-rich populations, which has been confirmed by UV spectroscopy. By way of example, Figure 12 compares the far UV spectrum of $\xi$ Per (O7.5 iii) with a counterpart in Sextans A [116], revealing negligible wind signatures in the latter (e.g., C iv $\lambda$1550).