An excess of massive stars in the local 30 Doradus starburst

See allHide authors and affiliations

Science  05 Jan 2018:
Vol. 359, Issue 6371, pp. 69-71
DOI: 10.1126/science.aan0106

Observing more massive stars

The number of stars that form at each mass is known as the initial mass function (IMF). For most masses, the IMF follows a power-law distribution, first determined by Edwin Salpeter in 1955. Schneider et al. used observations of the nearby star-forming region 30 Doradus (also known as the Tarantula Nebula) and combined these with stellar modeling to determine its IMF. They found more stars above 30 solar masses than predicted by the Salpeter distribution. Because the most massive stars also have the biggest influence on their surroundings—for instance, through ultraviolet radiation, stellar winds, supernova explosions, and production of heavy elements—this excess will have wide-ranging implications.

Science, this issue p. 69


The 30 Doradus star-forming region in the Large Magellanic Cloud is a nearby analog of large star-formation events in the distant universe. We determined the recent formation history and the initial mass function (IMF) of massive stars in 30 Doradus on the basis of spectroscopic observations of 247 stars more massive than 15 solar masses (Embedded Image). The main episode of massive star formation began about 8 million years (My) ago, and the star-formation rate seems to have declined in the last 1 My. The IMF is densely sampled up to 200 Embedded Image and contains 32 ± 12% more stars above 30 Embedded Image than predicted by a standard Salpeter IMF. In the mass range of 15 to 200 Embedded Image, the IMF power-law exponent is Embedded Image, shallower than the Salpeter value of 2.35.

Starbursts are large star-formation events whose feedback affects the dynamical and chemical evolution of star-forming galaxies throughout cosmic history (13). They are found at low and high redshift, with the earliest starburst galaxies contributing to the reionization of the universe (2, 4). In such starbursts, massive stars [≥10 solar masses (Embedded Image)] dominate the feedback through intense ionizing radiation, stellar outflows, and supernova explosions. Because of large distances to most starbursts, analyses have so far been restricted to either photometric observations or composite spectra of entire stellar populations. In the former case, the high surface temperature of massive stars precludes the determination of accurate physical parameters because their colors are too similar (5), and, in the latter case, physical parameters of individual stars cannot be determined (6). Greater understanding can be obtained by spectroscopically examining individual stars within star-forming regions.

The stellar initial mass function (IMF) influences many areas of astrophysics because it determines the relative fraction of massive stars; that is, those that undergo supernova explosions and drive the evolution of star-forming galaxies. Much effort has therefore gone into understanding whether the IMF is universal or varies with local environmental properties (7, 8). Over the past few decades, evidence has accumulated that the IMF slope may be flatter than that of a Salpeter IMF (9)—in other words, there are more high-mass stars than expected—in regions of intense star formation (1012). However, these studies are based on integrated properties of stellar populations, hampering the ability to infer IMFs.

The star-forming region 30 Doradus (30 Dor) lies within the Large Magellanic Cloud, a satellite galaxy of the Milky Way, and has a metallicity (total abundance of all elements heavier than helium) of ~40% of the solar value (13). At a distance of 50 kpc (14), 30 Dor is a nearby analog of distant starbursts and one of the brightest hydrogen-ionization (H ii) regions in the local universe (15). With a diameter of ~200 pc, 30 Dor hosts several star clusters and associations and is similar in size to luminous H ii complexes in more distant galaxies (16).

Through the use of the Fibre Large Array Multi Element Spectrograph (FLAMES) (17) on the Very Large Telescope (VLT), the VLT-FLAMES Tarantula Survey (VFTS) (18) has obtained optical spectra of ~800 massive stars in 30 Dor, avoiding the core region of the dense star cluster R136 because of difficulties with crowding (18). Repeated observations at multiple epochs allow determination of the orbital motion of potentially binary objects. For a sample of 452 apparently single stars, robust stellar parameters—such as effective temperatures, luminosities, surface gravities, and projected rotational velocities—are determined by modeling the observed spectra (19). Composite spectra of visual multiple systems and spectroscopic binaries are not considered here because their parameters cannot be reliably inferred from the VFTS data.

To match the derived atmospheric parameters of the apparently single VFTS stars to stellar evolutionary models, we use the Bayesian code Bonnsai, which has been successfully tested with high-precision observations of Galactic eclipsing binary stars (20). Bonnsai accounts for uncertainties in the atmospheric parameters and determines full posterior probability distributions of stellar properties, including the ages and initial masses of the VFTS stars (19). By summing these full posterior probability distributions of individual stars, we obtain the overall distributions of stellar ages and initial masses of massive stars currently present in 30 Dor (Fig. 1). These distributions are missing those stars that already ended their nuclear burning. However, given that we know both the present-day age and mass distributions, we can correct for these missing stars and derive the star-formation history (SFH) and IMF of massive stars in 30 Dor (19), allowing us to fully characterize this prototype starburst.

Fig. 1 Distributions of stellar ages and initial masses of massive stars in 30 Dor.

(A) Age and (B) initial mass (Mini) distribution of the VFTS sample stars more massive than 15 Embedded Image (black lines). Uncertainties are calculated by bootstrapping (19), and the 1σ regions are shaded blue. The best-fitting star-formation history (SFH) (A) and present-day distribution of initial masses (B) are plotted in red. For comparison, the expected present-day distribution of initial masses, assuming a Salpeter initial mass function (IMF), is also provided in (B) (note that these modeled mass distributions are no longer single power-law functions). About 140 stars above 15 Embedded Image are inferred to have ended their nuclear burning during the last ≈10 My, and their contribution to the SFH is shown by the red shaded region in (A). The peak star-formation rate (SFR) extrapolated to the entire 30 Dor region is ~Embedded Image (on the order of Embedded Image, depending on the exact size of 30 Dor). The age and mass distributions are number density functions (dN/dt and dN/dM) with respect to age t and mass M. γ is the power-law slope of the IMF. (C) Ratio of modeled to observed present-day mass functions, illustrating that the Salpeter IMF model underpredicts the number of massive stars in our sample, particularly above 30 Embedded Image.

When determining the SFH and IMF, it is necessary to account for selection biases. The VFTS target selection implemented a magnitude cut, observing only stars brighter than the 17th magnitude in the V band (18). Compared with a full photometric census of massive stars in 30 Dor (21), the VFTS sample is ~73% complete. Although, owing to the magnitude limit, the VFTS is incomplete for stars Embedded Image, the completeness shows no correlation with the V-band magnitude of stars more massive than 15 Embedded Image (19). Of the 452 stars with robust stellar parameters, 247 are more massive than 15 Embedded Image and form the basis of our determination of the SFH and high-mass end of the IMF. Incompleteness corrections are applied to account for our selection process (19). We assume the high-mass IMF is a power-law function, Embedded Image, where M is the mass and γ the slope, and compute the SFH and corresponding prediction of the distribution of initial masses for different IMF slopes until we best match (i) the number of stars above a given mass and (ii) the observed initial-mass distribution (19).

We find that the observed distribution of initial masses of stars in 30 Dor is densely sampled up to ~200 Embedded Image. It is shallower than that predicted by a Salpeter IMF with γ = 2.35, and the discrepancy increases with mass (Fig. 1C). Relative to the Salpeter index, we find an excess of Embedded Image (Embedded Image) stars more massive than 30 Embedded Image and Embedded Image (Embedded Image) stars more massive than 60 Embedded Image (Fig. 2 and fig. S5; unless stated otherwise, uncertainties are 68.3% confidence intervals). The hypothesis that a Salpeter IMF can explain the large number of stars more massive than 30 Embedded Image in our sample can thus be rejected with >99% confidence (19). The number of stars more massive than 30 Embedded Image is best reproduced by an IMF slope of Embedded Image (Fig. 2). Using our second diagnostic, a least-squares fit to the observed distribution of initial masses over the full mass range of 15 to 200 Embedded Image, our best fit is Embedded Image (Figs. 1 and 3), in agreement with our first estimate based on the number of massive stars ≥30 Embedded Image. Our high-mass IMF slope is shallower than the slope inferred by other studies for stars below ≈20 Embedded Image in the vicinity of R136 (22, 23).

Fig. 2 Number of massive stars in our sample.

Expected number of massive stars in our sample initially more massive than (A) 30 Embedded Image and (B) 60 Embedded Image as a function of the IMF slope γ (black solid lines). The blue and red shaded areas indicate the 68 and 95% confidence intervals (CIs), respectively, of the observed number of stars (see fig. S5). The IMF slopes best reproducing the observed number of stars and the associated 68% intervals are indicated by the vertical dashed lines and gray shaded regions and correspond to Embedded Image and Embedded Image for stars more massive than 30 Embedded Image and 60 Embedded Image, respectively.

Fig. 3 Probability density function of the inferred IMF slope in 30 Dor.

Results are based on a χ2 power-law fitting over the mass range of 15 to 200 Embedded Image. The shaded areas represent 1σ, 2σ, and 3σ confidence regions, and the slope of the Salpeter IMF is indicated by the vertical dashed line. Our inferred IMF is shallower than that of Salpeter (γ = 2.35), with 83% confidence.

The limitation of our sample to stars ≥15 Embedded Image means that we can reconstruct the SFH of 30 Dor over the past ≈12 million years (My). When also considering the 1- to 2-My-old stars in R136 that were not observed within VFTS (24), we find that the star-formation rate in 30 Dor sharply increased ~8 My ago and seems to have dropped ~1 My ago (Fig. 1A). If the currently observed drop continues for another million years, the duration of the main star-forming event will be shorter than ~10 My. This result complements a recent study (23) that found a similar time-dependence of star formation around the central R136 star cluster in 30 Dor on the basis of photometric data of low- and intermediate-mass stars. We therefore conclude that star formation in the 30 Dor starburst is synchronized across a wide mass range.

Our results challenge the suggested 150 Embedded Image limit (25) for the maximum birth mass of stars. The most massive star in our sample, VFTS 1025 (also known as R136c), has an initial mass of Embedded Image (19). From stochastic sampling experiments (19), we exclude maximum stellar birth masses of more than 500 Embedded Image in 30 Dor with 90% confidence because we would otherwise expect to find at least one star above 250 Embedded Image in our sample. Our observations are thus consistent with the claim that stars with initial masses of up to 300 Embedded Image exist in the core of R136 (26).

Approximately 15 to 40% of our sample stars are expected to be products of mass transfer in binary star systems (27). Binary mass transfer in a stellar population produces a net surplus of massive stars and rejuvenates stars such that they appear younger than they really are (28). Mass accretion alone biases the inferred IMF slope to flatter values, whereas rejuvenation steepens it. Taken together, we calculate that these two effects roughly cancel out in our case and, thus, binary mass transfer cannot explain the difference between our inferred IMF and that of Salpeter (19). Also, our final sample of stars contains unrecognized binaries, but they do not affect our conclusions (19).

The core of the R136 star cluster is excluded from the VFTS, but stars ejected from R136 (so-called runaway stars) may enter our sample. Runaway stars are biased toward high masses (29) and thus flatten the upper IMF. However, star clusters such as R136 typically eject about 5 to 10 stars above 15 Embedded Image (30, 31), which is insufficient to explain the expected excess of 25 to 50 stars above 30 Embedded Image in 30 Dor, after correcting for the completeness of our sample and that of the VFTS (19).

We conclude that the 30 Dor starburst has produced stars up to very high masses (Embedded Image), with a statistically significant excess of stars above 30 Embedded Image and an IMF shallower above 15 Embedded Image than a Salpeter IMF. Measuring the IMF slope above 30 to 60 Embedded Image has proven difficult (7), and, in general, large uncertainties in the high-mass IMF slope remain (32). This raises the question of whether star formation in 30 Dor proceeded differently. It has been suggested that starburst regions themselves provide conditions for forming relatively more massive stars by the heating of natal clouds from nearby and previous generations of stars (33). Alternatively, a lower metallicity may lead to the formation of more massive stars because of weaker gas cooling during star formation. An IMF slope shallower than the Salpeter value may then be expected at high redshift when the universe was hotter and the metallicity lower (33, 34).

Because massive-star feedback increases steeply with stellar mass, it is strongly affected by the IMF slope. Comparing an IMF slope of Embedded Image to the Salpeter slope, we expect Embedded Image more core-collapse supernovae and an increase of supernova metal-yields and hydrogen ionizing radiation by factors of Embedded Image and Embedded Image, respectively (19). The formation rate of black holes increases by a factor of Embedded Image (19), directly affecting the expected rate of black hole mergers detected by their gravitational wave signals. We also expect an increase in the predicted number of exotic transients that are preferentially found in starbursting, metal-poor dwarf galaxies such as long-duration gamma-ray bursts (35) and hydrogen-poor superluminous supernovae (36). Many population synthesis models and large-scale cosmological simulations assume an IMF that is truncated at 100 Embedded Image. Compared with such models, the various factors estimated above are even larger (19).

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S10

Tables S1 to S3

References (37116)

Data S1 to S3

References and Notes

  1. Materials and methods are available as supplementary materials.
  2. The Bonnsai web service is available at
  3. For the determination of atmospheric parameters of B-type stars, see supplementary materials section S3.2.
  4. For the determination of atmospheric parameters of A-type and later-type stars, see supplementary materials section S3.3.
  5. For the determination of atmospheric parameters of classical Wolf-Rayet stars, see supplementary materials section S3.1.
Acknowledgments: We thank the referees for constructive feedback that helped to improve this work. Our findings are based on observations collected at the European Southern Observatory (ESO) under program ID 182.D-0222. This work was supported by the Oxford Hintze Centre for Astrophysical Surveys, which is funded through generous support from the Hintze Family Charitable Foundation. H.S. acknowledges support from the FWO-Odysseus program under project G0F8H6N. G.G. acknowledges financial support from Deutsche Forschungsgemeinschaft grant GR 1717/5. O.H.R.-A. acknowledges funding from the European Union’s Horizon 2020 research and innovation program under a Marie Skłodowska-Curie grant agreement (665593) awarded to the Science and Technology Facilities Council (STFC). C.S.-S. acknowledges support from CONICYT-Chile through the FONDECYT Postdoctoral Project (grant 3170778). S.S.-D. and A.H. thank the Spanish Ministry of Economy and Competitiveness (MINECO) for grants AYA2015-68012-C2-1 and SEV2015-0548. S.E.d.M. has received funding under the European Union’s Horizon 2020 research and innovation program from the European Commission under a Marie Skłodowska-Curie grant agreement (661502) and the European Research Council (ERC) (grant agreement 715063). M.Ga. and F.N. acknowledge Spanish MINECO grants FIS2012-39162-C06-01 and ESP2015-65597-C4-1-R. M.Gi. acknowledges financial support from the Royal Society (University Research Fellowship) and the European Research Council (ERC StG-335936, CLUSTERS). R.G.I. thanks the STFC for funding his Rutherford fellowship under grant ST/L003910/1 and Churchill College, Cambridge, for his fellowship and access to their library. V.K. acknowledges funding from FONDECYT-Chile fellowship grant 3160117. J.M.A. acknowledges support from the Spanish government Ministerio de Economía y Competitividad (MINECO) through grant AYA2016-75 931-C2-2-P. N.M. acknowledges the financial support of the Bulgarian National Science Fund under grant DN08/1/13.12.2016. The Space Telescope Science Institute is operated by the Association of Universities for Research in Astronomy under NASA contract NAS5-26555. C.J.E. is also a visiting professor at the University of Edinburgh. The raw VFTS observations are available from ESO’s Science Archive Facility at under project ID 182.D-0222. A web interface for the Bonnsai software is available at Input and derived stellar parameters for all stars used in this study are provided in table S3 and in machine-readable form in data S1, the best-fitting SFH can be found in data S2, and our Python code for determining the stellar maximum birth masses is available in data S3.
View Abstract

Navigate This Article