Open Access
Issue
A&A
Volume 659, March 2022
Article Number A163
Number of page(s) 16
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/202243048
Published online 07 April 2022

© F. Martins and A. Palacios 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Very massive stars (VMSs) have masses in excess of 100 M (Vink 2015). It is an open question how far above this threshold the mass of stars extends. Analyses of the stellar content of young massive clusters in the Galaxy and the Magellanic Clouds have led to the prediction of an upper mass cut-off of 150–200 M (Weidner & Kroupa 2004; Oey & Clarke 2005; Figer 2005). Dynamical measurements of stellar masses in binary systems are so far consistent with this limit: none of the most massive stars observed in spectroscopic binaries has a mass higher than 150 M (Schnurr et al. 2008; Tehrani et al. 2019; Shenar et al. 2021). The current record holder for an eclipsing system, in which the dynamical mass can be probed directly, without an assumption about the inclination, is NGC3603-A1, in which the primary component has a present-day mass of 116 ± 31 M (Schnurr et al. 2008). Masses can also be obtained from stellar luminosity and comparison to theoretical predictions of evolutionary models. This can be done for binary components and for single stars (Martins 2015). Using this method, Crowther et al. (2010) reported masses up to 320 M in the Large Magellanic Cloud (LMC) cluster R136. Although these mass estimates have been revised since then, they are still higher than 200 M for a few objects (Bestenlehner et al. 2020).

Since VMSs evolve rapidly and because of the shape of the initial mass function (IMF), young massive star clusters are the preferred places to search for VMSs. Wofford et al. (2014) presented HST/COS ultraviolet (UV) spectroscopy of these clusters in the blue compact dwarf galaxy NGC 3125, at a metallicity close to that of the LMC. Standard population synthesis models with a Salpeter IMF extending up to 100 M were unable to reproduce all lines, in particular the broad He II 1640 emission. Because the spectrum of cluster NGC3215-A1 is similar to that of R136 in the LMC, where VMSs are observed, the authors concluded that stars with masses in excess of 100 M are a necessary ingredient of population synthesis models (see also Smith et al. 2016; Senchyna et al. 2017; Leitherer et al. 2018). Incorporating updated evolutionary tracks of VMSs in these models, Senchyna et al. (2021) confirmed that better agreement between observed and predicted equivalent width (EW) of UV lines is obtained for extreme star-forming regions. However, they relied on atmosphere models that are not specifically designed for VMSs. On even broader spatial scales, Goswami et al. (2021) studied the [O/Fe] ratio of Galactic thin- and thick-disk stars and concluded that very massive stars up to 350 M are required to explain the properties of the latter group of stars. It is thus necessary to develop tailored evolutionary and atmosphere models of VMSs to study young stellar populations in the Local Universe.

A key ingredient of modelling VMSs is their strong stellar winds. Martins et al. (2008) showed that the mass-loss rates of the most massive stars in the Galactic Arches cluster scaled with luminosity, albeit with a significant offset compared to the relation followed by O-type stars (see also Hainich et al. 2014). In parallel, Gräfener & Hamann (2008) investigated the driving of stellar winds in late-type WN stars. They concluded that very luminous, H-rich WNh stars developed optically thick radiatively driven winds mainly because of their proximity to the Eddington limit. Subsequent investigations confirmed that VMSs have higher mass-loss rates than more classical and less massive O-type stars (Vink et al. 2011; Gräfener et al. 2011). Using observations of young massive stars in 30 Dor, Bestenlehner et al. (2014) derived a mass-loss rate recipe for VMSs in the LMC (see also Bestenlehner 2020). Gräfener (2021) used a combination of this recipe and that of Vink et al. (2001) to study the impact of mass loss on the evolution of VMSs in the LMC. However, the spectral appearance was not studied. Conversely, Gräfener & Vink (2015) presented the first synthetic spectra of VMSs at different metallicities, but provided only snapshot spectra for a specific set of stellar parameters (only one value of effective temperature and luminosity).

In the present paper, we describe a pilot study aiming at providing synthetic spectra of VMSs along their evolutionary paths. We concentrate on an LMC-like metallicity to benefit from the mass-loss recipe of Bestenlehner et al. (2014) and Gräfener (2021). We first describe our evolutionary models and how they compare to previous studies (Sect. 2). After this, we present the calculation of synthetic spectra and describe the spectroscopic evolutionary sequences in the optical and UV spectral ranges (Sect. 3). We explore the impact of our models on population synthesis in Sect. 4, and we conclude in Sect. 5.

2 Evolutionary models

In this section, we present the main ingredients we used to calculate evolutionary tracks of VMSs. We describe in detail our implementation of the mass-loss recipe of Gräfener (2021) and then compare our results to grids of models that are available in the literature.

2.1 Physical ingredients

As in Martins & Palacios (2021), we computed non-rotating evolutionary models with the code STAREVOL (Siess et al. 2000; Lagarde et al. 2012; Amard et al. 2019). We modelled stars with initial masses of 150, 200, 250, 300, and 400 M. We recall here the main physical ingredients relevant for this study. We assumed an Eddington grey atmosphere as outer boundary condition to the stellar structure equations. We used the Asplund et al. (2009) solar chemical composition as a reference, with Z = 0.0134. A calibration of the solar model with these input physics leads to an initial helium mass fraction Y = 0.2687 at solar metallicity. We used the corresponding constant slope ΔYZ = 1.60 (with the primordial abundance Y0 = 0.2463 based on WMAP-SBBN by Coc et al. 2004) to compute the initial helium mass fraction at Z = 5.50 × 10−3, corresponding to [Fe/H] = −0.4 and typical of the LMC stars. We then scaled all the abundances accordingly. The OPAL opacities used for these models comply to this scaled distribution of nuclides. We did not include specific α-element enhancement in our models. We described the convective instability using the mixing-length theory with αMLT = 1.6304, and we used the Schwarzschild instability criterion to define the boundaries of convective regions. We added a step overshoot at the convective core edge and adopted αov = 0.1Hp, with Hp being the pressure scale height. We used the thermonuclear reaction rates from the NACRE II compilation (Xu et al. 2013b) for mass number A < 16, and those from the NACRE compilation (Angulo et al. 1999) for more massive nuclei up to Ne. The proton captures on nuclei more massive than Ne are from Longland et al. (2010) or Iliadis et al. (2001). The network was generated via the NetGen server (Xu et al. 2013a).

2.2 Mass-loss treatment

We implemented the mass-loss treatment prescribed by Gräfener (2021) in our models. The mass-loss recipe of VMSs depends on the optical thickness of their winds. While they appear as O stars, their winds are optically thin. However, when they enter the WNh phase (see Sect. 3.2.1 for a definition) while still on the main sequence, their winds become optically thick and are not properly described by the mass loss of classical core He-burning Wolf-Rayet stars.

Gräfener (2021) thus proposed a criterion to switch between the optically thin wind regime (O-type star) and the optically thick wind regime (WNh-type star). This criterion is tailored for LMC VMSs; it was fitted on observed mass-loss rates for VMSs in 30 Dor as derived by Bestenlehner et al. (2014). It is based on the sonic-point optical depth τs in the wind, which for a stellar model with surface mass M, luminosity L, radius R, and hydrogen mass fraction X writes τsM˙uL/c1+vesc2u2,${\tau_s} \approx {{\dot M{u_\infty }} \over {L/c}}\left({1 + {{v_{esc}^2} \over {u_\infty ^2}}} \right),$(1)

where ν and νesc are the terminal and escape velocity in the wind, respectively, and v/vesceff=2.51±0.27${v_\infty }/v_{esc}^{eff} = 2.51 \pm 0.27$ following the empirical relation for O stars with Teff > 25 000 K from Lamers et al. (1995). The escape velocity νesc is assumed to coincide with the effective escape velocity, uescuesceff=2GM(1Γe)/Reff,${u_{esc}} \equiv u_{esc}^{eff} = \sqrt {2GM(1 - {\Gamma_e})/{R_{eff}},}$(2)

where Reff is the effective radius such as L=4πReff2σTeff4$L = 4\pi R_{{\rm{eff}}}^2\sigma T_{eff}^4$, and Γe is the classical Eddington factor expressed in the case of fully ionized plasma, when the electron scattering opacity χe only depends on the hydrogen mass fraction X, Γe=χeL4πcGM=104.813(1+X)LLMM.${\Gamma_e} = {{{\chi_e}L} \over {4\pi cGM}} = {10^{ - 4.813}}(1 + X){L \over {{L_\odot }}}{{{M_{_\odot }}} \over M}.$(3)

For models at the LMC metallicity, the criterion to decide between the two wind regimes is the following:

if τs<23${\tau_s} < {2 \over 3}$, then the optically thin regime → thin

if τs > 1, then the optically thick regime → thick

In the optically thin regime, which is characteristic of the OB stars, the mass-loss rate is the rate derived by Vink et al. (2001) on the hot side of the bistability jump, that is, for 27 500 ≤ Teff ≤ 50000 K, log(M˙thin)=6.697(±0.061)+2.194(±0.021)×logL*1051.313(±0.046)×logM*301.226(±0.037)×logv/vesc2+0.933(±0.064)×logTeff4000010.92(±0.90)×logTeff400002+70.85(±0.10)×logZZ$\matrix{ {\log ({{\dot M}_{{\rm{thin}}}}) = } \hfill &amp; { - \,6.697(\pm \,0.061)} \hfill \cr {} \hfill &amp; { + \,2.194(\pm \,0.021) \times \log \left({{{{L_*}} \over {{{10}^5}}}} \right)} \hfill \cr {} \hfill &amp; { - \,1.313\,(\pm \,0.046) \times \log \left({{{{M_*}} \over {30}}} \right)} \hfill \cr {} \hfill &amp; { - \,1.226\,(\pm \,0.037) \times \log \,\left({{{{v_\infty }/{v_{esc}}} \over 2}} \right)} \hfill \cr {} \hfill &amp; { + \,0.933\,(\pm \,0.064) \times \,\log \,\left({{{{T_{eff}}} \over {40000}}} \right)} \hfill \cr {} \hfill &amp; { - 10.92\,(\pm \,0.90) \times \log {{\left({{{{T_{eff}}} \over {40000}}} \right)}^2}} \hfill \cr {} \hfill &amp; { + \,70.85\,(\pm \,0.10) \times \log \,\left({{Z \over {{Z_\circ }}}} \right)} \hfill \cr }$(4)

where is in M yr−1, L*, and M*. are the luminosity and mass of the star in solar units, Teff is the effective temperature in Kelvin, the ratio between terminal and escape velocities ν/νesc is taken equal to 2.6 as in Eq. (24) of Vink et al. (2001) (which is compatible within the error bars with the value of 2.51 used in Eq. (1)), and Z = 0.019.

This expression differs from that adopted by Gräfener (2021), who advocated two modifications: a global decrease in mass loss by 20%, and the use of a shallower slope of the mass loss as a function of L, with the replacement of the factor 2.194 in Eq. (4) by a factor 1.45. We decided not to apply the first modification because it was added in order to reproduce observed LMC trends, which is not our main goal here. We conserved the original factor by Vink et al. (2001) in the second modification because the modified factor of 1.45 proposed by Gräfener (2021) comes from the work of Bestenlehner et al. (2014), where it characterized the wind-momentum1 luminosity relation (see their Eq. (7) and Table 5) and not the mass-loss luminosity relation. In practice, only the earliest phase (above 53 000 K) of the 150 M track is in the optically thin regime: the remainder of the 150 M track as well as the entire track for masses above 200 M are spent in the optically thick regime. Hence the choice of the mass-loss recipe in the optically thin regime has little effect on our results.

In the optically thick regime, the mass-loss rate comes from Bestenlehner et al. (2014), log(M˙thick)=5.22×log(Γeff)0.5×log(D)2.6.$\log ({\dot M_{{\rm{thick}}}}) = 5.22\, \times \,\log ({\Gamma_{eff}}) - 0.5 \times \log (D) - 2.6.$(5)

with D = 10 the wind clumping factor, and Γeff is the effective Eddington factor for a rotating star taken at an intermediate latitude, Γeff=Γe+0.5×Γrot=Γe+Ω2R32GM${\Gamma_{{\rm{eff}}}} = {\Gamma_{\rm{e}}} + 0.5\,\, \times \,\,{\Gamma_{{\rm{rot}}}} = {\Gamma_{\rm{e}}} + {{{\Omega ^2}{R^3}} \over {2GM'}}$(6)

where Ω is the surface angular velocity. We focus on non-rotating models here so that Γeff = Γe. Finally, in the intermediate regime, when 2/3 ≤ τs ≤ 1, we assumed a smooth transition and determined the mass-loss rate by performing a linear interpolation between thin and thick.

2.3 Comparison to published grids and surface chemical evolution

In Fig. 1, we present the evolutionary paths of the non-rotating STAREVOL models as described in the sections above. We also show three models, for 200 M, 300 M, and 400 M computed with a classical mass-loss recipe2 combining Vink et al. (2001) and Nugis & Lamers (2000) when the models enter a Wolf-Rayet-like phase that we used in Martins & Palacios (2017). In these latter models, we did not include wind clumping in order to be consistent with the models computed by Köhler et al. (2015) with the BEC code (Heger et al. 2000; Brott et al. 2011). We show the latter models in Fig. 1 for comparison. Finally, we overplot the 200 M, 300 M, and 400 M non-rotating models computed by Gräfener (2021) with the MESA code (Paxton et al. 2019, and references therein). In order to alleviate the figure, we limited the tracks to parts in which the central hydrogen mass fractions are larger than about 0.01.

The tracks that include the classical mass-loss recipe shown here as dotted lines can be directly compared to those of Köhler et al. (2015) shown as dotted lines. The agreement is fairly good considering the different assumptions made in this study for the input physics. The STAREVOL models exhibit a redder extension of the main sequence than the BEC models, but this turning point is very sensitive to both mass loss and convection treatment, so that the observed differences are not a concern.

Our new models, including the mass-loss recipe described in Sect. 2.2, are in very good qualitative agreement with those of Gräfener (2021): the red bending of the upper the zero-age main sequence (ZAMS), the luminosity decrease along the main-sequence evolution as a result of the huge mass loss they experience, and the envelope inflation is recovered (see also Sect. 3.1). The ZAMS of our models is more blueshifted than Graefener’s due to different assumptions for the initial chemical composition. The models with M ≥ 200 M experience optically-thick winds already on the ZAMS, which explains the direct decrease in their luminosity beyond this evolutionary point. The 150 M model experiences a longer phase with optically thin winds, and thus the slope of its path in the HRD is much more similar to what is obtained when the classical mass-loss recipe is used.

Because all models above 200 M experience only optically thick winds3, the mass-loss rates predicted by the new mass-loss scheme are high from the ZAMS on, typically a factor 3 to 10 higher than for the optically thin recipe of Vink et al. (2001) (see e.g. Fig. 8 of Bestenlehner et al. 2014). This leads to an efficient stripping of the stellar envelope, as can be seen from Table 1. This leads to the early exposure of the products of hydrogen nucleosynthesis through the CNO cycle at the surface. Helium is strongly enhanced, as shown in Table 1 and Fig. 2. In this figure, we show the evolutionary tracks in the HRD for the models including the Gräfener (2021) (left panel) and classical mass loss (right panel). The surface 4He mass fraction is colour-coded along the tracks. Comparison between the left and right panels of this figure and with the values reported for the 200 M models in Tables 1 and 2 indicates that the helium enrichment is much higher than predicted by the models using the classical mass-loss recipe, with important consequences for the spectral appearance of the low-metallicity VMSs, as we discuss in the following sections.

thumbnail Fig. 1

Hertzsprung-Russell diagram (HRD) for our models with two different mass-loss treatments compared to those published in Köhler et al. (2015) and Gräfener (2021). See text for more details. Only the main-sequence part of the tracks up to Xc ≈ 0.01 is shown for all models. The long dashed gray line represents the ZAMS of the STAREVOL models..

Table 1

Parameters of the atmosphere models.

3 Spectroscopic properties

In this section, we first describe the computation of atmosphere models and synthetic spectra. We then discuss the general morphology of our theoretical spectra and compare them to spectra of real stars in the LMC. We also discuss the effect of the new mass-loss rate prescription on synthetic spectra by comparing our results to models using the recipe of Vink et al. (2001).

3.1 Atmosphere models and synthetic spectra

Along each evolutionary sequence, we selected points for which we computed an atmosphere model and the associated synthetic spectrum with the code CMFGEN (Hillier & Miller 1998). With applications to population synthesis in mind, we chose points at 0, 0.5, 1.0, 1.5, and 2.0 Myr along each track (when possible, the most massive stars living less than 2.0 Myr). A final point was taken at the redmost extension of the track. Fig. 2 shows the positions of the selected points.

We took the surface properties of the stars predicted by STAREVOL as input for the atmosphere and synthetic spectra calculations. These parameters are gathered in Table 1. The terminal velocity was set to 2.51 × υesc as assumed in Sect.2.2. Because the mass-loss recipe of Gräfener (2021) takes wind clumping into account, we also added inhomogeneities in the computation of the atmosphere models. In CMFGEN, this is done by means of a volume filling factor f that is assumed to vary as f=f+(1f)ev/vcl,$f = {f_\infty } + (1 - {f_\infty }){e^{ - v/{v_{{\rm{cl}}}}}},$(7)

where υ is the wind velocity, υcl is a characteristic velocity at which clumping becomes non-negligible, and f is the maximum clumping factor. We adopted standard values for the last two parameters: υcl = 100 km s−1 and f = 0.1 corresponding to D = 10 in Eq. (5).

The following elements were included in the calculations: H, He, C, N, O, Ne, Mg, Si, P, S, Ar, Ca, Fe, and Ni. The abundances were taken from the evolutionary model output. To compute the atmosphere structure, a microturbulent velocity of 20 km s−1 was adopted. For the synthetic spectrum, the microturbulent velocity varied from 10 km s-1 at the photosphere to 0.1 × υ at the top of the atmosphere. The spectra were calculated from 50 Å to 3 μm.

Different assumptions are made to compute the temperature and atmosphere structures in the outer regions of evolutionary models and in the inner regions of atmosphere models. It is not obvious that these structures are consistent (see Groh et al. 2014; Martins & Palacios 2017). We confirmed this consistency for all our models and found that the agreement is excellent on the ZAMS and degrades in more advanced phases. However, the agreement in our most advanced models is still reasonable. This is illustrated in Fig. 3, where we took the ZAMS and most evolved models of the 200 M0 series (models 200a and 200f) as examples. For model 200a, the internal and atmosphere structures merge almost perfectly. For the 200f model, there is a small offset, but the structures still intersect. We tested that this is no longer the case in more advanced phases. We thus chose not to compute spectra beyond the last model of each track in Fig. 2. This would require a proper treatment of the atmospheric boundary in evolutionary models and is beyond the scope of this study. Similarly, we ran tests for a 500 M model and found that a strong discrepancy between evolutionary and atmospheric structures is present already on the ZAMS. We thus limited our calculations to masses ≤ 400 M.

Interestingly, we also note in Fig. 3 that the density structures of our models displays an plateau in the outer regions. This phenomenon is usually referred as “envelope inflation” (Ishii et al. 1999; Gräfener et al. 2012) and is naturally taken into account in our models. This can be explained by our treatment of subsurface convection in the evolutionary models (see also Gräfener 2021).

thumbnail Fig. 2

HRD for the STAREVOL models on the main sequence. The most evolved tracks (150 and 200 M) are shown up to a central hydrogen mass fraction equal to 0.015. The evolution of the surface helium mass fraction is colour-coded on the tracks. The left panel shows the tracks obtained using the new mass-loss recipe by Gräfener (2021), and the right panel shows those using the old mass-loss recipe combining Vink et al. (2001) and Nugis & Lamers (2000), as in Martins & Palacios (2017). The evolutionary points at which atmosphere models and synthetic spectra have been computed are indicated with squares..

Table 2

Same as Table 1 for the 200 M models computed with the Vink et al. mass-loss rate recipe.

thumbnail Fig. 3

Temperature and density structure as a function of radius in models 200a (left) and 200f (right). In each panel, the red line is the structure of the evolutionary model, and the black line is the structure of the atmosphere model..

3.2 Spectral morphology

We illustrate the evolution of the spectroscopic appearance of our VMS models with the 200 M star in Fig. 4. The figures for the other masses are gathered in Appendix A. To better identify spectral lines, we show in Fig. 5 the contribution of individual ions to the total UV spectrum of model 200d and to the optical spectrum of model 200f. These two models were chosen because they allow the identification of the maximum number of individual lines in the UV and optical range, respectively. We distribute our synthetic spectra through the POLLUX database4 (Palacios et al. 2010).

In the next sections, we describe the morphology of our VMS spectra in the optical and UV range. We highlight the features that appear as typical of VMSs and are not observed in stars with M ≲ 100 M.

3.2.1 Optical range

In the optical range, it is immediately clear that the spectra are dominated by He II 4686 emission that increases as the star evolves. For the 200 M models, the EW of this emission ranges from −9 to −55 Å5. Other He II lines are predicted in absorption near the ZAMS and in emission in more advanced phases. Hydrogen lines of the Balmer series are seen either in absorption or in full emission, depending on the line and the state of evolution. In more evolved models, the contribution of H I lines weakens because the hydrogen content is reduced (see Table 1). This implies that the emission features at the position of the Balmer hydrogen lines become dominated by He II lines as evolution proceeds; see for instance the emission at the position of Hβ in model 200f (Fig. 5, right panel).

Nitrogen lines from N III, N IV, and/or N V are the next strongest features in the optical range. Their relative strength varies depending on the ionization of the atmosphere: N V is present at the highest temperatures and recombines into N IV and N III when Teff decreases and the wind density increases. For instance, N V 4605-20 and N IV 4058 are predicted in model 200d, where N III 4640 is very weak. In model 200f, N V 4605-20 has almost vanished, while N III 4640 is comparable to N IV 4058. In addition to ionization effects, the rapid increase in nitrogen content at the surface of the stars causes a global strengthening of all nitrogen lines. Wind effects also play a role (Rivero González et al. 2011, 2012): the presence of velocity fields leads to a desaturation of some key transitions that control the populations of the energy levels of optical lines. Consequently, the strength of these lines strongly depends on the wind properties. We stress that N III is the main contributor to the emission complex located near the position of Hδ in model 200f (and similarly evolved models). Hi 4101 and Si IV 4089-4116 emissions are present, but much weaker than N III 4097.

These spectral properties are typical of either O supergiants with strong winds (Olf stars), transition objects (Olf/WN), or hydrogen-rich WN (WNh) stars, as witnessed by the emission in He II 4686 (Martins 2018; Smith et al. 1996). The morphology of Hβ can be used to further distinguish between early OIf, OIf/WN, and WNh stars, as demonstrated by Crowther & Walborn (2011): Hβ is in absorption, with a P-Cygni profile, and in emission at these spectral types. Near the ZAMS, none of our models shows He I lines. Together with the Hβ morphology, this indicates a spectral type O2-3.5If/WN5-7 for the early phases of the 200 M star evolution. As Teff decreases and the wind strengthens, the spectral appearance shifts to WN6-7h. This trend is present for other masses and is modulated by the wind strength. The 150 M series has a lower wind density on average than the 200 M series, implying a start as O2-3.5If instead of O2-3.5If/WN5-7.

We predict a morphology of optical lines that is typical of VMSs observed in 30 Dor, for instance (see Bestenlehner et al. 2020 and Sect. 3.3). Less massive O stars on the main sequence do not show these emission features (Sota et al. 2011, 2014). He II 4686 is usually in absorption, and none of the hydrogen Balmer lines shows emission (except sometimes Ha). Additional examples of LMC O star optical spectra can be found in the supplementary material of Bestenlehner et al. (2020).

thumbnail Fig. 4

UV (left panel) and optical (right panel) synthetic spectra of the M = 200 M series. The spectra are normalized and shifted upward for clarity. The main lines are indicated. Evolution proceeds from top to bottom..

thumbnail Fig. 5

UV spectrum of the 200d model (left) and optical spectrum of the 200f model (right). In each panel, the black line is the total spectrum, and coloured lines indicate the contribution of individual ions..

3.2.2 UV range

The typical UV spectrum of our VMS models is made of several P-Cygni profiles originating from different ions (Fig. 5, left panel). As in the optical range, He II lines are relatively strong. He II 1640 starts as an emission line on the ZAMS (model 200a) and develops rapidly (from model 200c) into the strongest P-Cygni profile of the UV range as the star evolves (see Fig. 4). The increase in wind density and of the surface helium content causes this behavior. With the exception of the first two models of the 150 M series (Fig. A.1), He II 1640 is always at least comparable in strength to C IV 1550, and often stronger. Its full width at half maximum is about 1000 km s−1. Its EW can reach up to −30 Å. These properties make He II 1640 a unique feature of VMSs on the main sequence.

C IV 1550 is present in all spectra with a prominent P-Cygni profile. Figure 4 shows that its strength remains roughly constant along the sequence. The expected increase due to the higher wind density is counter-balanced by the reduction of the carbon surface abundance by more than a factor of 10 between models 200a and 200f. The EW of the absorption (emission) part of C IV 1550 is 14 (−10 Å) and 9 (−9 Å) in these two models, respectively.

N IV 1720 is weak near the ZAMS, but rapidly develops into a strong P-Cygni profile, mainly because of the fast nitrogen enrichment of the surface (increase by a factor of 10 between models 200a and 200c). The intensity of N IV 1720 is comparable to or even exceeds that of C IV 1550.

N IV 1486 is predicted in emission towards the end of the series. It becomes one of the strongest features of VMS spectra. It is not observed in 40–85 M O-type supergiants with strong winds (OIf stars) even at solar metallicity (Walborn et al. 1985; Bouret et al. 2012). However, it is a rather common feature of WN and WNh stars (Hamann et al. 2006; Hainich et al. 2014).

At the bottom of Fig. 5, the Fe V (and to a lesser extent, Fe VI) line forests shape the morphology of the spectrum around 1400 Å (1300 Å). These forests appear as broad features in the most advanced phases of our models. The Fe V forest reduces the emission peak of the Si IV 1400 doublet, which would otherwise appear as a well-developed double P-Cygni profile.

Finally, Ar V 1194 is a strong line partly blended with the weaker S V 1188 feature. Its proximity to Lyα makes it difficult to actually observe because of broad interstellar absorption (as seen in Figs. 6 and 10). Nonetheless, to the best of our knowledge, this feature is not mentioned in atlases of O and Wolf-Rayet stars and might be a metallicity indicator in low-extinction regions.

In view of this general morphology, we conclude that a strong He II 1640 emission is a clear signature of VMSs. The presence N IV 1486 is another specific feature, and N IV 1720 can also be associated with VMSs when its strength is comparable to that of other P-Cygni profiles in the UV range.

thumbnail Fig. 6

UV spectrum of the WN5h stars R136 al and RMC 146 (blue middle spectra) with our models 250c and 300b..

3.3 Comparison to LMC stars

To determine the relevance of our synthetic spectra for representing the general morphology of VMSs at a metallicity representative of the LMC, we compare them with the normalized spectrum of the WN5h star R136 a1 in Fig. 6. The data are from Crowther et al. (2016) (see their Fig. B1). According to Bestenlehner et al. (2020), R136 a1 has a luminosity of 106.79 ±0.10 L. Hence we overplot the spectra of our 300 M and 250 M series that resemble the spectrum of R136 a1 most. The morphology of the synthetic spectra is similar to that of the observed spectrum. The strong P-Cygni profiles of N V 1240 and C IV 1550 are well accounted for. He II 1640 is the strongest line in models 300b and 250c, as observed for R136 a1. S V 1500 and O V 1371 are predicted as too strong in our models. Si IV 1400 is also present in our models.

Another sanity check was performed with the WN5h star RMC 146 in the LMC6. He II 1640 is observed as the strongest emission line as in models 300b and 250c. C iv 1550 and N IV 1720 have comparable emission peaks, as in the models, although the absolute emission level is slightly lower than predicted. For N IV 1720, this can partly be due to a different initial nitrogen content: we used a scaled solar value, which is probably higher than the baseline LMC nitrogen abundance (Dopita et al. 2019). S V 1500 and Si Iv 1400 are not observed in RMC 146. Si IV 1400 is sensitive to Teff (and thus age), and so is O V 1371, which is observed to be weaker than predicted in model 300b.

Since we did not aim at a detailed fit of these two stars but at a morphological comparison, we conclude that our models are satisfactory to account for the general UV features seen in R136 a1 and RMC 146. When optical data for these objects are publicly available, we will perform a similar comparison with our models. A preliminary qualitative inspection of Fig. S12 of Bestenlehner et al. (2020) indicates that our synthetic spectra look quite similar to the observed spectrum of R136 a1.

3.4 Comparison to Vink et al.’s mass-loss rates

The effect of the new mass-loss rate recipe on the spectra of VMSs can be assessed from Table 2 and Fig. 7. We gathered the parameters of models computed with the mass loss prescription of Vink et al. (2001) for the 200 M series (200V series in the following). For these calculations, we used the 200 M track presented in the right panel of Fig. 2. We stress that these computations do not include clumping. In the computations using the mass-loss recipe of Gräfener (2021), clumping is taken into account for optically thick winds, but not for optically thin winds, see Sect. 2.3. Clumping would have affected mostly the red wing of some emission lines (e.g. He II 4686 and Balmer lines; see Hillier 1991; Martins et al. 2009). It would also have improved the predicted shape of some UV features, in particular, P V 1118-28, O V 1371, and N IV 1720 (Bouret et al. 2005; Fullerton et al. 2006). For our computations, we chose the same ages as in the original series to compute the synthetic spectra. Due to the reduced mass-loss rates, the 200 M track extends towards lower Teff (see Fig. 2).

With the Vink et al. mass-loss rates, He II 1640 does not develop into a strong P-Cygni profile and remains relatively weak compared to other UV lines (Fig. 7, left panel). Its EW (emission part) barely reaches −1.5 Å in model 200Vd. The same conclusion applies to He II 4686 (emission EW smaller than 1.2 Å). This is due to a combination of a lower mass-loss rate and a lower helium content. Table 2 shows that is between 0.1 and 0.6 dex lower in the 200V series than in the initial 200 models. The helium mass fraction never exceeds 0.3 in the 200V models, while it reaches 0.85 in the initial models, as also illustrated in Fig. 2.

In contrast to the initial 200 M sequence, C iv 1550 increases along the 200V sequence, except for the last model, in which the low Teff implies a low C iv fraction. This trend is explained by the increase in wind density and a relatively constant surface carbon abundance, at least in models 220Va, 200Vb, and 200Vc. Only in model 200Ve, at 2 Myr, is the drop in abundance sufficient to counterbalance the wind effect. In the initial series, the C abundance drop takes place much earlier, before 1 Myr.

With the Vink et al. mass-loss rates, the 200 M track extends towards lower Teff where the Si Iv 1400 doublet can develop into strong P-Cygni profiles. In the initial series, Teff is higher, and as described previously, the Fe V line forest develops into relatively broad features that reduce the Si IV 1400 lines.

Another key difference is the shape of the iron line forests. With the weaker mass-loss rates implied by the Vink et al. recipe, these lines are formed closer to the photosphere where wind velocity is low. Hence, lines do not overlap because of Doppler shifts. This is different when the stronger mass-loss rates of Graefener are used. In this case, lines are formed higher in the atmosphere where velocities imply Doppler shifts comparable to the intrinsic line width. Hence blending creates the broad absorption features seen in the left panel of Fig. 4, for example.

The spectrum of VMSs is thus very different depending on the mass-loss rate prescription. Comparisons to LMC stars (see previous section) revealed that models need to produce He II 1640 at least as strong as C IV 1550 to correctly account for the observed UV morphology. Inspection of Fig. 7 shows that this is not achieved with the Vink et al. recipe. Our study thus favours the recipe of Gräfener (2021) to model VMSs.

thumbnail Fig. 7

Same as Fig. 4, but for models with the Vink et al. (2001) mass-loss rate recipe..

4 Effects on population synthesis

In this section, we investigate the impact of VMS spectra on population synthesis models incorporating stars up to 300 M. We show that a proper treatment of stellar evolution and stellar atmospheres with the mass-loss recipe relevant for VMSs has a strong impact on the appearance of young stellar populations. We first discuss theoretical models and how we add VMSs, before making comparisons to young massive clusters.

4.1 Theoretical predictions

To build population synthesis models incorporating VMSs, we retrieved the spectral energy distribution (SED) of the BPASS single-star population synthesis models7 for a metallicity Z = 0.006, a Salpeter IMF above 0.5 M, and an instantaneous burst of star formation. The metallicity is the closest to Z = 1/2.5 Z in the BPASS models. We refer to Eldridge et al. (2017) for a description of the BPASS models.

We first used BPASS models that extend up to 100 M. To assess the effect of VMSs on population synthesis models, we added to these BPASS models the contribution of stars in the mass range 100–300 M. For this, we used the same mass function as in the BPASS calculations and added the contribution of stars in the mass bins 100–175, 175–225, 225–275, and 275– 300 M0 as detailed in Appendix B. For each mass bin, we used the SEDs of the 150, 200, 250, and 300 M series. We focused on two ages − 1 and 2 Myr - because VMSs rapidly disappear after this.

Figure 8 shows the results of this exercise in the UV wavelength range. VMSs contribute almost the same amount of light as all 0.1−100 M stars. Some lines are almost exclusively produced by VMSs. This includes He II 1640 and N IV 1720. P V 1117-28, C iv 1169, N v 1240, and C iv 1550 are also affected, they are stronger when VMSs are included, especially at 2 Myr. At this age, N IV 1486 appears in the integrated spectrum because of VMSs.

The BPASS database provides models that include stars up to 300 M. For this purpose, both the evolutionary tracks and the stellar atmospheres for 100–300 M stars were computed with the mass-loss rates of Vink et al. (2001). This causes differences compared to population synthesis models based on the appropriate treatment of VMS winds. We illustrate this in Fig. 9. All wind lines are stronger when proper VMS mass-loss rates are included. This is again striking for He II 1640 and N IV 1720. The former is produced exclusively when our VMS models are taken into account. C IV 1550 is also affected, being stronger in our predictions. At 2 Myr N IV 1486 is present only in our models. The strength of all lines affected by VMSs depend on the number of such objects. We recall that in our tests we extended the Salpeter mass function up to 300 M. A different IMF slope and/or a different upper mass cut-off would modulate the shape of all strong wind-sensitive lines. Our tests are also made for a burst of star formation, and further studies should investigate various star formation histories.

Senchyna et al. (2021) produced improved population synthesis models with the Bruzual & Charlot code (Vidal-García et al. 2017). Unlike the BPASS models, they include new evolutionary tracks from Chen et al. (2015) that take the specific mass-loss rates of VMSs into account. The mass-loss recipe used by Chen et al. is that of Vink et al. (2011) and is slightly different from that of Grafener (2021). In particular, the scaling of the mass-loss rate with the Eddington factor is steeper in the Grafener (2021) recipe. When comparing our 250 M track with the 250 M track of Chen et al. (2015) at Z = 0.0068, we find that on the main sequence, our mass-loss rates are 0.3 dex higher on average than those of the Chen et al. track.

According to Senchyna et al. (2021), the updated Bruzual & Charlot models rely on the PoWR9 grids of Wolf-Rayet models for synthetic spectra (Todt et al. 2015). Hence dedicated VMS synthetic spectra are not included. In particular, this implies that the mass-loss rates and helium mass fraction used for the synthetic spectra are not consistent with those of the evolutionary tracks of Chen et al. (2015). Nonetheless, Senchyna et al. (2021) argued that VMSs are necessary to properly reproduce the UV spectrum of extreme star-forming galaxies. In their models, the EW of He II 1640 and C IV 1550 can only be matched with VMSs. The absorption (emission) EW of line C IV 1550 (He II 1640) in our population synthesis models is 8.6 (−4.7) Å at 1 Myr and 8.3 (−9.0 Å) at 2 Myr. From Fig. 4 of Senchyna et al. (2021), this is even larger than what the Bruzual & Charlot models predict, although they use a constant formation history rather than a burst of star formation, which partly explains the difference. The differences in mass-loss rates may also partly explain that we predict stronger lines.

thumbnail Fig. 8

UV spectrum of single-star BPASS population synthesis models at 1 Myr (left panel) and 2 Myr (right panel) shown in black. The blue spectrum is the contribution of stars with masses in the range 100–300 M. The red line is the spectrum of the entire population. In each panel the comparison between the normalized spectrum of the population with (red) and without (black) VMS is shown at the bottom..

thumbnail Fig. 9

Normalized UV spectrum of population synthesis models. Black lines are the BPASS model extending to 300 M. The red line is the combination of the BPASS model up to 100 M and the additional contribution of VMSs from our models and a Salpeter IMF slope for 100-300 M. The top (bottom) panel is for a burst of star formation after 1 (2) Myr..

thumbnail Fig. 10

Normalized UV spectrum of the R136 cluster (black) with our population synthesis model including VMS at 1 Myr (top, red) and the BPASS model for an upper mass cut-off of 300 M (bottom, green). The synthetic spectra have been normalized and degraded to a spectral resolution of 1500..

4.2 Young massive clusters

In this section we now investigate whether our population synthesis models allow a better fit to the observed spectrum of young massive clusters.

The first direct test is shown in Fig. 10. We compare the integrated UV spectrum of the R136 cluster in the LMC (Crowther et al. 2016) to our simple population synthesis. We also show the BPASS predictions with an upper mass cut-off of 300 M. Our test model matches the observed He II 1640 better than the BPASS spectrum. C IV 1550 is also stronger in our model and agrees better with the observed profile. O IV 1340 and O V 1371 are slightly too strong in our model compared to the spectrum of R136. The comparison of Si IV 1400 is not straightforward without an assessment of the interstellar absorption. We note that the IMF of the 30 Dor region, of which R136 is part, may be top heavy (Schneider et al. 2018), while we used a Salpeter mass function. Changing the IMF would increase the contribution of VMS and strengthen some of the emission lines. In spite of these shortcomings, we conclude that the morphology of the spectrum including VMS is closer overall to what is observed than the pure BPASS spectrum.

We further investigated whether our improved population synthesis model might reproduce the integrated spectrum of the super star cluster A1 in NGC3125 better. Wofford et al. (2014) showed that some of the strongest UV features, including He II 1640, could not be reproduced by classical STARBURST99 models (Leitherer et al. 2014). A comparison with our models is relevant since the metallicity of NGC3125 is estimated to be close to that of the LMC (Hadfield & Crowther 2006). We retrieved the HST/STIS G140L spectrum of the cluster from the MAST archive10. The spectrum has a resolution of about 1000 and covers the range 1150–1700 Å. We normalized the spectrum and cut the spectral region around Lyα because it is strongly absorbed. We show how our population synthesis models including VMS at 2 Myr compares to the spectrum of NGC3125 A1 in Fig. 11. The agreement is good, with Si IV 1400 and He II 1640 being relatively well accounted for given that we do not explore the full parameter space (age, IMF shape, star formation history). In particular, we are able to reproduce the broad He II 1640 emission, which was the main problem encountered by Wofford et al. (2014). N V 1240 and N IV 1486 are present in both our model and the observation, but we predict a slightly too strong emission. As already discussed, this may be due to the adopted initial nitrogen content of our models. Unfortunately, N IV 1720 is not observed and cannot be used to investigate this issue. We also note that N IV 1486 is very sensitive to age, appearing and increasing after ~1.5 Myr (see Fig. 4). C iv 1550 is also slightly too strong in our model. Future investigation with proper population synthesis models will help understand if this can be solved by tuning model parameters. The key points of our comparison are that we can reproduce He II 1640 without invoking an extremely top-heavy mass function, as in Wofford et al. (2014), and that we predict N IV 1486 emission (although it is too strong).

Hadfield & Crowther (2006) studied the optical spectrum of the clusters in NGC3125. They used the classical blue and red bumps near 4650 and 5800 Å, respectively, to infer the ratio of WN to WC stars. They relied on template spectra of classical Wolf-Rayet stars. Here we caution that the presence of VMSs may alter the WN/WC ratios determined by Hadfield & Crowther (2006) since both bumps are present in the spectra of VMSs. VMSs could thus contribute a significant fraction of the flux at the corresponding wavelength, reducing the contribution of classical WR stars. Here again, future studies with population synthesis models including VMSs are necessary. From our population synthesis experiments and comparisons to the spectra of young massive clusters, we conclude that it is important to include a proper treatment of VMS evolution and atmospheres to correctly predict the integrated light of young starburst, especially in the UV spectral range.

thumbnail Fig. 11

Normalized UV spectrum of the NGC3125 A1 cluster (black) with our population synthesis model including VMS at 2 Myr (red). The synthetic spectrum has been normalized and degraded to a spectral resolution of 1000 typical of the STIS/G140L grating..

5 Concluding remarks

We have presented evolutionary tracks for VMSs including the mass-loss recipe of Grafener (2021). We used the code STAREVOL. We focused on the main sequence for initial masses of 150, 200, 250, 300, and 400 M. These tracks were compared to previous calculations and agree well. VMSs loose a significant amount of mass during their main-sequence evolution, so that they exhibit CNO-cycle processed material, including large amounts of helium, very early on during this evolutionary phase.

We subsequently computed atmosphere models with the code CMFGEN. We selected points equally spaced in age along each evolutionary track. We showed that VMSs have an optical and UV spectrum that is dominated by He II lines for almost the entire main sequence. He II 4686 is the strongest optical line in all phases, and He II 1640 is stronger than any other UV line shortly after the star leaves the ZAMS. In the optical range, nitrogen lines are the second characteristic features. Balmer hydrogen lines are predicted in emission for most of the evolution, making the spectra of our VMS models typical of WNh stars. In the UV, a well-developed N IV 1720 P-Cygni profile, the presence of broad features caused by the Fe V forest, and emission in N IV 1486 are features specific to VMSs. Our synthetic spectra are morphologically similar to stars RMC 146 and R136 a1 in the LMC.

We tested the impact of VMS on the integrated spectrum of a young starburst. We added the contribution of VMSs with masses between 100 and 300 M0 to a BPASS population synthesis model at Z = 0.006. We showed that taking the spectra of VMSs properly into account is crucial, especially in the UV range, where VMSs dominate the integrated spectrum: He II 1640 emission is entirely caused by the presence of VMSs, N IV 1720 is much stronger when VMSs are included; C IV 1550 is also affected; and VMSs produce N IV 1486 emission after ~ 1–1.5 Myr. VMSs contribute almost as much light as the rest of the population in that wavelength range. Our population synthesis test model including VMSs reproduces the integrated spectrum of R136 better than the BPASS model that extends to 300 M. The reason is the proper treatment of mass loss in both the evolutionary models and synthetic spectra of VMSs. Including VMSs also helps to reproduce the strong UV lines of the super star cluster A1 in NGC3125.

In our study, we have worked with single stars. Some of the most massive stars known to date are members of binary |systems: NGC3602-A1 (Schnurr et al. 2008), R145 (Schnurr et al. 2009), Arches-F2 (Lohr et al. 2018), Melnick 34 (Tehrani et al. 2019), and R144 (Shenar et al. 2021). At least one of the components usually shows a spectrum typical of a WNh star. Including binary stars in population synthesis models extending up to 300 M is thus necessary. However, we anticipate that our conclusion regarding the need for a proper treatment of the wind physics of VMSs remains at least qualitatively valid. Integrated spectra of young star-forming galaxies are usually obtained at relatively low spectral resolution where the component of the massive binary systems listed above are blended. The light is thus dominated by the most massive component, or in case of nearly equal-mass systems, is the sum of the two VMS components.

In this study, we have focused on a metallicity Z = 1/2.5 Z that can be viewed as typical of the LMC. We chose this metallicity because the mass-loss recipe for VMSs was mainly established from the analysis of LMC stars (Bestenlehner et al. 2020; Gräfener 2021). The LMC also hosts the most massive stars that can be observed individually and to which we can compare our predictions. In future studies, we will extend our predictions to a wider range of metallicities. We will also test the ability of our population synthesis models to reproduce the UV spectrum of young star-forming galaxies, in particular, their stellar He II 1640 emission (e.g. Erb et al. 2010; Saxena et al. 2020; Marques-Chaves et al. 2020, 2021).

Acknowledgements

We thank an anonymous referee for a positive report and suggestions of clarification. We warmly thank John Hillier for making the code CMFGEN available to the community and for constant help with it. This work made use of v2.2.1 of the Binary Population and Spectral Synthesis (BPASS) models as last described in Stanway & Eldridge (2018). We thank Paul Crowther and Joachim Bestenlehner for sharing the UV spectra of R136 a1 and of the entire R136 population. We acknowledge interesting discussions with Daniel Schaerer and Rui Marques-Chaves about the UV spectra of young starburst galaxies. Some of the data presented in this paper were based on observations made with the NASA/ESA HST, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). The STScI is operated by the association of Universities for Research in Astronomy, Inc. under the NASA contract NAS 5-26555. These observations are associated with program 12172. This research has made use of NASA’s Astrophysics Data System.

Appendix A Synthetic spectra throughout all mass sequences

In Figs. A.1 to A.4 we show the UV and optical spectra along the 150, 250, 300, and 400 M evolutionary tracks.

thumbnail Fig. A.1

Synthetic UV (left) and optical (right) spectra of the 150 M models.

thumbnail Fig. A.2

Synthetic UV (left) and optical (right) spectra of the 250 M models.

thumbnail Fig. A.3

Synthetic UV (left) and optical (right) spectra of the 300 M models.

thumbnail Fig. A.4

Synthetic UV (left) and optical (right) spectra of the 400 M models.

Appendix B Population synthesis

The BPASS SEDs are given for a cluster containing 106 M. According to the BPASS user manual11, the number of stars between 0.1 M0 and a maximum mass M2 is N(0.1,M2)=C×(0.1M1Mα1dM+M1α1M1M2Mα2dM).$N(0.1,{M_2}) = C \times \left({\int_{0.1}^{{M_1}} {{M^{{\alpha_1}}}} dM + M_1^{{\alpha_1}}\int_{{M_1}}^{{M_2}} {{M^{{\alpha_2}}}} dM} \right).$(B.1)

For our comparisons, we selected the BPASS models with M1 = 0.5 M, α1 = −1.3, α2 = −2.35, and M2 = 100 M. The total mass of the cluster is given by M(0.1,M2)=C×(0.1M1M1+α1dM+M1α1M1M2M1+α2dM).$M(0.1,{M_2}) = C \times \left({\int_{0.1}^{{M_1}} {{M^{1 + {\alpha_1}}}} dM + M_1^{{\alpha_1}}\int_{{M_1}}^{{M_2}} {{M^{1 + {\alpha_2}}}} dM} \right).$(B.2)

The normalization constant C has a value of 1.23 × 105 for a total mass equal to 106. For this set up, the number of stars with masses above 100 M and in the mass bin [Mα;Mb] is N(Ma,Mb)=C×M1α1MaMbMσ2dM,$N({M_a},{M_b}) = C \times M_1^{{\alpha_1}}\int_{{M_a}}^{{M_b}} {{M^{{\sigma_2}}}} dM,$(B.3)

which gives 236.56, 60.30, 35.42, and 12.16 stars for the mass bins 100–175, 175–225, 225–275, and 275–300 M0, respectively. To take VMSs (M> 100 M) in population synthesis into account, we thus added the spectra of the 150, 200, 250, and 300 M multiplied by these numbers to the SED provided by BPASS (for a cluster hosting stars with masses between 0.1 and 100 M). We did this at two ages: 1 and 2 Myr.

References

  1. Amard, L., Palacios, A., Charbonnel, C., et al. 2019, A&A, 631, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [Google Scholar]
  3. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bestenlehner, J. M. 2020, MNRAS, 493, 3938 [Google Scholar]
  5. Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918 [Google Scholar]
  7. Bouret, J. C., Lanz, T., & Hillier, D. J. 2005, A&A, 438, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bouret, J. C., Hillier, D. J., Lanz, T., & Fullerton, A. W. 2012, A&A, 544, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068 [Google Scholar]
  11. Coc, A., Vangioni-Flam, E., Descouvemont, P., Adahchour, A., & Angulo, C. 2004, ApJ, 600, 544 [NASA ADS] [CrossRef] [Google Scholar]
  12. Crowther, P. A., & Walborn, N. R. 2011, MNRAS, 416, 1311 [Google Scholar]
  13. Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731 [Google Scholar]
  14. Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624 [Google Scholar]
  15. Dopita, M. A., Seitenzahl, I. R., Sutherland, R. S., et al. 2019, AJ, 157, 50 [Google Scholar]
  16. Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058 [Google Scholar]
  17. Erb, D. K., Pettini, M., Shapley, A. E., et al. 2010, ApJ, 719, 1168 [Google Scholar]
  18. Figer, D. F. 2005, Nature, 434, 192 [Google Scholar]
  19. Fullerton, A. W., Massa, D. L., & Prinja, R. K. 2006, ApJ, 637, 1025 [Google Scholar]
  20. Goswami, S., Slemer, A., Marigo, P., et al. 2021, A&A, 650, A203 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gräfener, G. 2021, A&A, 647, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Gräfener, G., & Hamann, W. R. 2008, A&A, 482, 945 [Google Scholar]
  23. Gräfener, G., & Vink, J. S. 2015, A&A, 578, L2 [Google Scholar]
  24. Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Groh, J. H., Meynet, G., Ekström, S., & Georgy, C. 2014, A&A, 564, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hadfield, L. J., & Crowther, P. A. 2006, MNRAS, 368, 1822 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hainich, R., Ruhling, U., Todt, H., et al. 2014, A&A, 565, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Hamann, W. R., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hillier, D. J. 1991, A&A, 247, 455 [NASA ADS] [Google Scholar]
  32. Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
  33. Iliadis, C., D’Auria, J. M., Starrfield, S., Thompson, W. J., & Wiescher, M. 2001, ApJS, 134, 151 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ishii, M., Ueno, M., & Kato, M. 1999, PASJ, 51, 417 [NASA ADS] [Google Scholar]
  35. Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Kudritzki, R. P., Puls, J., Lennon, D. J., et al. 1999, A&A, 350, 970 [Google Scholar]
  37. Lagarde, N., Decressin, T., Charbonnel, C., et al. 2012, A&A, 543, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269 [Google Scholar]
  39. Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14 [Google Scholar]
  40. Leitherer, C., Byler, N., Lee, J. C., & Levesque, E. M. 2018, ApJ, 865, 55 [NASA ADS] [CrossRef] [Google Scholar]
  41. Lohr, M. E., Clark, J. S., Najarro, F., et al. 2018, A&A, 617, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Longland, R., Iliadis, C., Champagne, A. E., et al. 2010, Nucl. Phys. A, 841, 1 [NASA ADS] [CrossRef] [Google Scholar]
  43. Marques-Chaves, R., Álvarez-Marquez, J., Colina, L., et al. 2020, MNRAS, 499, L105 [NASA ADS] [CrossRef] [Google Scholar]
  44. Marques-Chaves, R., Schaerer, D., Álvarez-Márquez, J., et al. 2021, MNRAS, 507, 524 [NASA ADS] [CrossRef] [Google Scholar]
  45. Martins, F. 2015, in Astrophysics and Space Science Library, 412, Very Massive Stars in the Local Universe, ed. J. S. Vink, 9 [NASA ADS] [CrossRef] [Google Scholar]
  46. Martins, F. 2018, A&A, 616, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Martins, F., & Palacios, A. 2017, A&A, 598, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Martins, F., & Palacios, A. 2021, A&A, 645, A67 [EDP Sciences] [Google Scholar]
  49. Martins, F., Hillier, D. J., Paumard, T., et al. 2008, A&A, 478, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Martins, F., Hillier, D. J., Bouret, J. C., et al. 2009, A&A, 495, 257 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227 [NASA ADS] [Google Scholar]
  52. Oey, M. S., & Clarke, C. J. 2005, ApJ, 620, L43 [NASA ADS] [CrossRef] [Google Scholar]
  53. Palacios, A., Gebran, M., Josselin, E., et al. 2010, A&A, 516, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  55. Rivero González, J. G., Puls, J., & Najarro, F. 2011, A&A, 536, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Rivero González, J. G., Puls, J., Massey, P., & Najarro, F. 2012, A&A, 543, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Saxena, A., Pentericci, L., Mirabelli, M., et al. 2020, A&A, 636, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69 [NASA ADS] [CrossRef] [Google Scholar]
  59. Schnurr, O., Casoli, J., Chené, A. N., Moffat, A. F. J., & St-Louis, N. 2008, MNRAS, 389, L38 [NASA ADS] [CrossRef] [Google Scholar]
  60. Schnurr, O., Moffat, A. F. J., Villar-Sbaffi, A., St-Louis, N., & Morrell, N. I. 2009, MNRAS, 395, 823 [Google Scholar]
  61. Senchyna, P., Stark, D. P., Vidal-García, A., et al. 2017, MNRAS, 472, 2608 [NASA ADS] [CrossRef] [Google Scholar]
  62. Senchyna, P., Stark, D. P., Charlot, S., et al. 2021, MNRAS, 503, 6112 [NASA ADS] [CrossRef] [Google Scholar]
  63. Shenar, T., Sana, H., Marchant, P., et al. 2021, A&A, 650, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593 [Google Scholar]
  65. Smith, L. F., Shara, M. M., & Moffat, A. F. J. 1996, MNRAS, 281, 163 [NASA ADS] [CrossRef] [Google Scholar]
  66. Smith, L. J., Crowther, P. A., Calzetti, D., & Sidoli, F. 2016, ApJ, 823, 38 [NASA ADS] [CrossRef] [Google Scholar]
  67. Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011, ApJS, 193, 24 [Google Scholar]
  68. Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014, ApJS, 211, 10 [Google Scholar]
  69. Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75 [NASA ADS] [CrossRef] [Google Scholar]
  70. Tehrani, K. A., Crowther, P. A., Bestenlehner, J. M., et al. 2019, MNRAS, 484, 2692 [NASA ADS] [CrossRef] [Google Scholar]
  71. Todt, H., Sander, A., Hainich, R., et al. 2015, A&A, 579, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Vidal-García, A., Charlot, S., Bruzual, G., & Hubeny, I. 2017, MNRAS, 470, 3532 [Google Scholar]
  73. Vink, J. S. 2015, in Astrophysics and Space Science Library, Vol. 412, Very Massive Stars in the Local Universe [CrossRef] [Google Scholar]
  74. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Vink, J. S., Muijres, L. E., Anthonisse, B., et al. 2011, A&A, 531, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Walborn, N. R., Nichols-Bohlin, J., & Panek, R. J. 1985, NASA Reference Publication, 1155 [Google Scholar]
  77. Weidner, C., & Kroupa, P. 2004, MNRAS, 348, 187 [NASA ADS] [CrossRef] [Google Scholar]
  78. Wofford, A., Leitherer, C., Chandar, R., & Bouret, J.-C. 2014, ApJ, 781, 122 [Google Scholar]
  79. Xu, Y., Goriely, S., Jorissen, A., Chen, G. L., & Arnould, M. 2013a, A&A, 549, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Xu, Y., Takahashi, K., Goriely, S., et al. 2013b, Nucl. Phys. A, 918, 61 [Google Scholar]

1

The modified wind momentum is defined by M˙υR*/R${\rm{\dot M}}\upsilon \infty \sqrt {{R_*}/{R_\odot }}$ following Kudritzki et al. (1999).

2

This mass-loss treatment is similar to that referred to as the Dutch scheme in the MESA code.

3

For the 150 M model, the winds are optically thin down to 53 000 K and optically thick after that.

5

We adopt the usual convention of negative EW for emission.

6

The observed spectrum of RMC 146 is SWP04140 in the IUE database of the MAST archive (https://archive.stsci.edu/iue/search.php).

8

This value of Z is the closest in the Chen et al. grid to our adopted value, Z = 0.0055.

All Tables

Table 1

Parameters of the atmosphere models.

Table 2

Same as Table 1 for the 200 M models computed with the Vink et al. mass-loss rate recipe.

All Figures

thumbnail Fig. 1

Hertzsprung-Russell diagram (HRD) for our models with two different mass-loss treatments compared to those published in Köhler et al. (2015) and Gräfener (2021). See text for more details. Only the main-sequence part of the tracks up to Xc ≈ 0.01 is shown for all models. The long dashed gray line represents the ZAMS of the STAREVOL models..

In the text
thumbnail Fig. 2

HRD for the STAREVOL models on the main sequence. The most evolved tracks (150 and 200 M) are shown up to a central hydrogen mass fraction equal to 0.015. The evolution of the surface helium mass fraction is colour-coded on the tracks. The left panel shows the tracks obtained using the new mass-loss recipe by Gräfener (2021), and the right panel shows those using the old mass-loss recipe combining Vink et al. (2001) and Nugis & Lamers (2000), as in Martins & Palacios (2017). The evolutionary points at which atmosphere models and synthetic spectra have been computed are indicated with squares..

In the text
thumbnail Fig. 3

Temperature and density structure as a function of radius in models 200a (left) and 200f (right). In each panel, the red line is the structure of the evolutionary model, and the black line is the structure of the atmosphere model..

In the text
thumbnail Fig. 4

UV (left panel) and optical (right panel) synthetic spectra of the M = 200 M series. The spectra are normalized and shifted upward for clarity. The main lines are indicated. Evolution proceeds from top to bottom..

In the text
thumbnail Fig. 5

UV spectrum of the 200d model (left) and optical spectrum of the 200f model (right). In each panel, the black line is the total spectrum, and coloured lines indicate the contribution of individual ions..

In the text
thumbnail Fig. 6

UV spectrum of the WN5h stars R136 al and RMC 146 (blue middle spectra) with our models 250c and 300b..

In the text
thumbnail Fig. 7

Same as Fig. 4, but for models with the Vink et al. (2001) mass-loss rate recipe..

In the text
thumbnail Fig. 8

UV spectrum of single-star BPASS population synthesis models at 1 Myr (left panel) and 2 Myr (right panel) shown in black. The blue spectrum is the contribution of stars with masses in the range 100–300 M. The red line is the spectrum of the entire population. In each panel the comparison between the normalized spectrum of the population with (red) and without (black) VMS is shown at the bottom..

In the text
thumbnail Fig. 9

Normalized UV spectrum of population synthesis models. Black lines are the BPASS model extending to 300 M. The red line is the combination of the BPASS model up to 100 M and the additional contribution of VMSs from our models and a Salpeter IMF slope for 100-300 M. The top (bottom) panel is for a burst of star formation after 1 (2) Myr..

In the text
thumbnail Fig. 10

Normalized UV spectrum of the R136 cluster (black) with our population synthesis model including VMS at 1 Myr (top, red) and the BPASS model for an upper mass cut-off of 300 M (bottom, green). The synthetic spectra have been normalized and degraded to a spectral resolution of 1500..

In the text
thumbnail Fig. 11

Normalized UV spectrum of the NGC3125 A1 cluster (black) with our population synthesis model including VMS at 2 Myr (red). The synthetic spectrum has been normalized and degraded to a spectral resolution of 1000 typical of the STIS/G140L grating..

In the text
thumbnail Fig. A.1

Synthetic UV (left) and optical (right) spectra of the 150 M models.

In the text
thumbnail Fig. A.2

Synthetic UV (left) and optical (right) spectra of the 250 M models.

In the text
thumbnail Fig. A.3

Synthetic UV (left) and optical (right) spectra of the 300 M models.

In the text
thumbnail Fig. A.4

Synthetic UV (left) and optical (right) spectra of the 400 M models.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.